Generic equilibria are derived for turbulent relaxing plasmas via an entropymaximization procedure that accounts for the shorttime conservation of certain collisionless invariants. The conservation of these collisionless invariants endows the system with a partial memory of its prior conditions but is imperfect on long time scales due to the development of a turbulent cascade to small scales which breaks the precise conservation of phase volume making this memory imprecise. The equilibria are still determined by the shorttime collisionless invariants but the invariants themselves are driven to a universal form by the nature of the turbulence. This is numerically confirmed for the case of beam instabilities in onedimensional electrostatic plasmas where sufficiently strong turbulence appears to cause the distribution function of particle energies to develop a universal powerlaw tail with exponent 2.
['Physics - Plasma Physics', 'Astrophysics - High Energy Astrophysical Phenomena', 'Condensed Matter - Statistical Mechanics']
Relaxation to universal nonMaxwellian equilibria in a collisionless plasma
{'Relaxation to universal non-Maxwellian equilibria in a collisionless plasma': "R. J. Ewart a,b,1 , M. L. Nastac a,c , P. J. Bilbao d , T. Silva d , L. O. Silva d , and A. A. Schekochihin a,e  \nThis manuscript was compiled on June 31, 2024  \nGeneric equilibria are derived for turbulent relaxing plasmas via an entropy-maximization procedure that accounts for the short-time conservation of certain collisionless invariants. The conservation of these collisionless invariants endows the system with a partial 'memory' of its prior conditions, but is imperfect on long time scales due to the development of a turbulent cascade to small scales, which breaks the precise conservation of phase volume, making this memory imprecise. The equilibria are still determined by the short-time collisionless invariants, but the invariants themselves are driven to a universal form by the nature of the turbulence. This is numerically confirmed for the case of beam instabilities in one-dimensional electrostatic plasmas, where sufficiently strong turbulence appears to cause the distribution function of particle energies to develop a universal power-law tail, with exponent -2 .  \nplasma | statistical mechanics | non-thermal particle acceleration  \nT he naive application of statistical mechanics would imply that one should find a universe full to the brim with plasma in local Maxwell-Boltzmann equilibrium. In fact, the plasmas in our Universe are not, by and large, in such an equilibrium. This is true in many settings, from cosmic rays populating our Galaxy and beyond (1), down through the solar neighbourhood (2, 3), and all the way to Earth-based experiments (4-6). It is not at all difficult, of course, to motivate why this should be the case. The process by which plasmas are driven towards Maxwell-Boltzmann equilibrium-two-body interparticle collisions (7-10)typically takes place over a time scale much longer than that associated with the evolution set by the mean (averaged over length scales far larger than the interparticle separation) electric and magnetic fields that the plasma inherently generates. In systems with large-scale (system-size) gradients, these fields may be driven unstable, triggering not just a zoo but an entire ecosystem of plasma instabilities (11, 12), with fluctuations growing and reaching amplitudes at which they can react back upon their progenitors, potentially altering the mean distribution and relaxing the system towards some semblance of stability, but by no means Maxwellianity. Despite all this, the distributions that we observe do possess some universal features (or fall into a finite number of universality classes (13-16)).  \nIt has therefore been a question of interest since the earliest observations of non-thermal distributions whether they can be explained by an overarching theory of relaxation. One such theory, proposed by Lynden-Bell (17) in the context of selfgravitating systems, gave hope of an affirmative answer to this question by appealing to the statistical-mechanical principle of maximum entropy even in the absence of particle collisions to enforce it. Unlike the collisional approach, Lynden-Bell's theory made use of the fact that the distribution function f ( x , v ) of a collisionless system naturally obeys a Liouville equation and, therefore, conserves an infinite family of invariants: the conservation of 'phase volume' in position and velocity space ( x , v ) implies that, for any function G ( f ),  \n∫ d x d v G ( f ) = const . [1]  \nTherefore, in line with the dictates of statistical mechanics, the conservation of this continuum of quantities, known as 'Casimir invariants' (18), must be respected when the entropy of the system (appropriately defined) is maximized-effectively giving the system an enhanced memory of its initial conditions. Despite the promise of this theory, its actual application has been limited due to the inherent complexity and apparent non-universality of the resulting equilibria (19-23). A conceptually graver concern is the reliance on the precise conservation of the Casimir invariants. In any real system, their conservation was only ever going to be approximate, with molecular chaos eventually reasserting itself. Aggravating this concern further,  \n|  \n|  \n|  \nno. XX  \n|", 'Significance Statement': "A fundamental tenet of thermodynamics is that chaotic systems will relax to maximum-entropy states. In plasmas, the chaos is conventionally provided by interparticle collisions and the universal maximumentropy equilibrium is a Maxwellian distribution. However, in collisionless plasmas, the chaotic state is due to collective turbulent dynamics. In this paper, we argue theoretically, and show numerically, that such plasmas still relax towards universal equilibria, but are non-Maxwellian, featuring a power-law distribution of particle energies. These equilibria are still obtained via a maximumentropy principle, but taking into account the conservation, on short time scales, of phase volume in a collisionless system, and a subtle form of breaking of this conservation on longer time scales, due to the nature of collisionless plasma turbulence.  \nAuthor affiliations: a Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford, OX1 3PU, UK; b Balliol College, Oxford, OX1 3BJ; c St John's College, Oxford, OX1 3JP; d GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisbon, Portugal; e Merton College, Oxford, OX1 4JD  \nR.J.E., L.O.S., and A.A.S. designed research; R.J.E., M.L.N., P.J.B., T.S., L.O.S, and A.A.S. contributed new reagents/analytic tools; R.J.E., M.L.N., P.J.B, and T.S. analyzed data; R.J.E., M.L.N., P.J.B., T.S., L.O.S., and A.A.S. wrote the paper.  \nThe authors have no competing interests to declare.  \nFig. 1. Self-consistent particle-in-cell simulation of the temporal evolution (from top to bottom) of the phase-space density f ( x, v ) for the electron-only two-stream instability, visualised across three time snapshots: 50 ω -1 pe , 500 ω -1 pe , and 2500 ω -1 pe (where ω -1 pe is the electron plasma frequency). The two counter-propagating beams, of speed v b , (top in red) rapidly go unstable generating the phase-space structure shown at 50 ω -1 pe . As the system evolves (middle and bottom snapshots), the formation and merger of turbulent electron holes is clearly evident, which drives the mixing of the phase space and relaxation towards a collisionless equilibrium. The color scale is normalised to the peak value of f ( x, v ) in each snapshot.  \n<!-- image -->  \nx  \npe/  \nv  \nb  \nit has been demonstrated both theoretically (24-26) and numerically (27, 28) that even for systems with nominally weak interparticle collisions-such that the time scale for relaxation to a Maxwellian equilibrium is much longer than dynamical times-the Casimir invariants are broken on time scales competitive with the evolution of the system (provided there is a sufficient level of turbulence to stir it), rendering their status as invariants questionable at best.  \nIn this paper, we show that turbulent plasmas do achieve Lynden-Bell equilibria, and that, rather than doing this in spite of the breaking of Casimir invariants, they manage it in tandem with this. We further show that the LyndenBell equilibrium that is achieved has a universal high-energy asymptotic. We carry out this study the aid of the particlein-cell (PIC) code OSIRIS (29), applied to one of the most classically studied turbulent collisionless systems: the twostream instability. The very earliest simulations (30, 31) showed vividly that this instability leads a system composed of long-lived structures in phase space (cf. 'BGK modes' (32)), known as electron 'holes' (see figure 1). These holes, since studied analytically, numerically, and observationally (see, e.g., (33, 34) and references therein) move around, merge, and generically represent a relaxing turbulent state. It is precisely this turbulence that drives the system's relaxation towards the Lynden-Bell equilibrium. However, the same turbulence generates small velocity- and position-space structures (in the manner of a turbulent cascade predicted by (35)), which cause the Casimir invariants to evolve with time. In the course  \n|  \nof this evolution, the underlying Lynden-Bell equilibrium to which the system wants to relax is gradually changed and the system adjusts to reach this evolving target equilibrium. Thus, the system's precise memory of initial conditions is replaced with a 'turbulent amnesia': turbulent fluctuations are perpetually trying to push the system towards its LyndenBell equilibrium, but the goalposts are continually moved by the breaking of Casimir invariants driven by those same fluctuations. Eventually the Casimir invariants themselves reach a steady state, previously conjectured by (36), which causes the system to converge to a final, universal LyndenBell equilibrium. This equilibrium exhibits a particle-energy distribution possessing a power-law tail with exponent -2, confirming the existence of universal equilibria in strongly turbulent relaxing plasmas.", 'Lynden-Bell equilibria': "The formalism of Lynden-Bell statistical mechanics is built around the argument that, during relaxation, the plasma becomes highly chaotic and disordered in phase space. A prototypical example of this, a plasma experiencing a twostream instability (the simulation and statistical mechanics of which will be discussed below), is shown in figure 1: the plasma starts as two cold counter-propagating electron beams in phase space and is mixed until a myopic eye cannot distinguish the forest (mean phase-space density) from the trees (exact phase-space density). Because of this disorder, and in light of the conservation of Casimir invariants [1], it is  \n<!-- image -->  \nFig. 2. Left panel: the evolution of the mean distribution function N ( ε ) of particle energies ε = mv 2 / 2 for the electron two-stream instability. The energies are normalized to the initial beam-particle energy ε b = mv 2 b / 2 ; n 0 is the mean particle number density. The solid black line shows a Maxwellian of the same energy as the initial condition. The lines shown in shades of the blue-to-red color palette are the distributions from the initial condition (dark blue) to t = 3200 ω -1 pe (dark red) at successive time intervals of 50 ω -1 pe . Right panel: the evolution of the 'waterbag content' ρ ( η ) , over time, computed from [7] for the fine-grained phase-space density f ( x , v ) , as detailed in the Materials and Methods. The color scheme for time evolution is the same as in the left panel. The universal asymptotics of N ( ε ) and ρ ( η ) are plotted as dashed lines showing the development of a power-law tail [17] and the agreement with the predicted scaling of the waterbag content [15].  \n<!-- image -->  \nhelpful to consider the exact phase-space density f ( x , v ) to be a random quantity and to describe the statistical state of the plasma by the probability density P ( v , η ) (which can be defined over a large spatial average) that, at a given velocity v , this phase-space density f ( x , v ) takes the value η (37, 38). Once P ( v , η ) is known, the mean phase-space density 〈 f 〉 can always be computed as the first moment in η of P ( v , η ):  \n〈 f 〉 ( v ) = ∫ d η ηP ( v , η ) . [2]  \nTo determine P ( v , η ), in equilibrium, one assumes that the system maximizes the Lynden-Bell entropy  \nS LB = -∫ d v d η P ( v , η ) ln P ( v , η ) [3]  \nsubject to reasonable constraints ∗ . The most obvious such constraint is the conservation of probability:  \n∫ d η P ( v , η ) = 1 . [4]  \nThe next most unobjectionable constraint is that the total energy of the system is fixed:  \n∫ d v 1 2 m | v | 2 ∫ d η ηP ( v , η ) = E tot , [5]  \nwhere we have neglected contributions from the potential energy (which we will find to be smaller by a factor ∼ 100 in our simulations). The final constraint that must be enforced is the conservation of phase volume, mathematically stated as  \n∫ d v P ( v , η ) = ρ ( η ) . [6]  \nThis simply reads that the volume-integrated probability of the exact phase-space density f taking the value η is equal to some prescribed function of η . This function, referred to by us as the 'waterbag content', can be computed from the exact phase-space density of the plasma:  \nρ ( η ) = 1 V ∫∫ d x d v δ ( f ( x , v ) -η ) , [7]  \nwhere V is the spatial volume of the system. In effect, ρ ( η ) measures the phase-space volume occupied by each level set η of f . This indeed captures all Casimir invariants [1], which can be recovered as weighted integrals of ρ ( η ) † :  \n∫ d x d v G ( f ) = V ∫ d η G ( η ) ρ ( η ) . [8]  \nThus, P ( v , η ) naturally enables the book-keeping of the level sets of the phase-space density. There is a subtle complication (discussed in greater detail in (36, 39)) as to how η = 0 (the empty level set) is handled, since ρ ( η ) cannot be finite at η = 0 for systems with unbounded velocity  \ndomains. Treating this subtlety carefully and maximising the entropy [3] subject to the constraints [4-6] gives the LyndenBell equilibria  \nP ( v , η ) = δ ( η ) + e -βηε ( v ) F ( η ) 1 + ∫ η> 0 d η ' e -βη ' ε ( v ) F ( η ' ) , [9]  \nwhere ε ( v ) = m | v | 2 / 2 and β and F ( η ) are Lagrange multipliers, known, in analogy with conventional thermodynamics, as the inverse 'thermodynamic temperature' and the 'fugacity', respectively. They are specified by enforcing the constraints [5] and [6] for the function [9] (which has already been normalized to satisfy [4]).  \n'Turbulent amnesia'. The formalism of Lynden-Bell is, therefore, now completely prescriptive: in specifying the initial condition, one chose the energy of the system E tot and the waterbag content ρ ( η ) through [7]. Strict adherence to Lynden-Bell's theory would then imply that the system should relax towards the appropriately solved equilibrium [9]. This is however, not the case. In any real situation, including any numerical simulation, one will find that the mean phase-space density will continue evolving even after the initial instability is quenched-and the putative collisionless invariants ρ ( η ) will evolve as well. This is manifest in figure 2, showing the energy distribution N ( ε ) of particles, and their waterbag content ρ ( η ) for a plasma undergoing the nearly collisionless two-stream instability depicted in figure 1. It is, of course, clear why ρ ( η ) should change: no truly collisionless plasma (or simulation) can exist, the Liouville equation is never perfectly satisfied, and so the memory of the initial conditions encoded by ρ ( η ) cannot be preserved forever. Why ρ ( η ) evolves relatively fast, even for nearly collisionless systems, will be discussed below. For the purposes of formulating a theoretical scenario and verifying it numerically, it suffices to know that ρ ( η ) does evolve. We therefore propose an amendment to the LyndenBell theory: the principle of 'turbulent amnesia'. Under this scheme, the collisionless dynamics of the plasma push the system towards the Lynden-Bell equilibrium computed using the time-evolving ρ ( η ), which is changed as the turbulence scrambles the system's long-term memory of its prior conditions. This is somewhat analogous to the way in which collisional plasmas undergoing heating will pass through a sequence of Maxwellian distributions with distinct temperatures even though the energy is not a conserved quantity: systems still strive to maximize entropy rapidly even when their invariants are imperfect. We shall return to the theoretical justification of this scenario once we have established its validity numerically.", 'Numerical verification': "Since plasmas expected to relax to Lynden-Bell equilibria should be as collisionless as possible, the best testbed currently feasible for such relaxation is 1D-1V simulations, which afford the largest number of particles per cell, and, therefore, the lowest collisionality due to the reduced noise levels (40, 41). As well as this, one can only expect the Lynden-Bell relaxation to be relevant when the plasma is sufficiently turbulent because inherent in the expectation of rapid maximization of entropy is the assumption of nearperfect mixing in phase space (we shall return below to the  \n|  \nFig. 3. Comparison of the directly measured mean phase-space density of the system (solid lines) and the Lynden-Bell equilibrium obtained by taking the first η moment [2] of [9] (dashed black lines) at two different times. At early times (blue), the mean phase-space density closely follows the Maxwellian associated with ρ ( η ) computed from the initial condition (the two beams, which have a single value of η ), while at late times (red), it is better fit by the Lynden-Bell equilibrium associated with the evolved ρ ( η ) shown in figure 2 (right panel). The difference between the Lynden-Bell equilibrium and the simulated distribution is highlighted.  \n<!-- image -->  \nquestion of how perfect it really is). We therefore test our relaxation scenario for two-beam plasmas in 1D-violently unstable situations for which the saturated state is not likely to be obtained via quasilinear theory (cf. (42)) ‡ . The beam instabilities that we study are described in detail in the Materials and Methods: the electron-only two-stream instability and the electron-positron two-stream instability. In the interest of simplicity, we have neglected the ion dynamics in the electron-only two-stream instability. Undoubtedly, this will prove to be an oversimplification at late times and on long length scales (missing the possibility of a plethora of ionscale physics: see, e.g., (34)), but it is the simplest possible framework to demonstrate our theory.  \nIn figures 1 and 2, we see that, at early times, ρ ( η ) is well approximated by a single delta function, in line with the initial condition  \nρ ( η, t = 0) = 2∆ v b δ ( η -η max ) , [10]  \nwhere ∆ v b is the beam width. At these early times, the mean phase-space density, also shown in figure 2, becomes very nearly Maxwellian. As time progresses, however, the mean phase-space density deviates from the Maxwellian equilibrium, forming an energy distribution with a power-law tail N ( ε ) ∼ ε -2 , while its waterbag content develops a lowη asymptotic ρ ( η ) ∼ η -1 . For the measured values of ρ ( η ) and total kinetic energy, we can now treat equation [9] as a black box (the numerical method for this is detailed in  \nthe Materials and Methods) and compare its output to the measured mean phase-space density. Such a comparison for early and late times is shown in figure 3 for the the electrononly two-stream instability. We have done this for initial beam widths ∆ v b between 2% and 50% of the beam speed v b and found good agreement with [9] up to around 25% of v b , beyond which point the instability is insufficiently violent to drive strong relaxation over the times and domain sizes that were simulated.  \nQuantitatively, this is the fundamental result of this paper. Viewed as a purely thermodynamic tool, the Lynden-Bell statistical mechanics correctly predicts the relaxation of the mean phase-space density in a two-stream unstable plasma. Qualitatively, however, it is possible to deduce how this thermodynamic tool operates, and why its output may be generic to many turbulent relaxing systems.  \nUniversal Lynden-Bell equilibria. We first note that the mean distribution [9] has a form very similar to the Fermi-Dirac distribution, owing to the analogy between the Pauli exclusion principle and the incompressibility of phase space: two fermions cannot be in the same quantum state and two level sets of phase space cannot be forced together. In a further analogy, the Lynden-Bell equilibrium can be roughly construed as the competition between two pieces of physics: the tendency of each level set to form a Maxwellian equilibrium η by η (as in the numerator of [9]) and the incompressibility of phase space (expressed by the denominator of [9] being greater than unity).  \nWhen the incompressibility condition wins this competition (which is the classical analogue of the system being nearly degenerate), the system is very close to its ground state, known, in plasma physics, as the Gardner state (43). The existence of such ground states is the inevitable consequence of an exclusion principle. The 'Gardner-restacked' minimum-energy counterpart to any given f is a distribution function f G ( ε ) that has ρ ( η ) identical to that of f , but is a monotonically decreasing function of solely the particle energy ε = m | v | 2 / 2. In this way, any further reduction of energy would require the distribution function to be larger at lower velocities, impossible without the compression of phase space, which is forbidden. Therefore, the ground states f G can be computed implicitly from-and are in one-to-one correspondence with- the waterbag content ρ ( η ):  \nρ ( η ) = ∫ d v δ ( f G -η ) 1D = -2 d f -1 G d η . [11]  \nThe last equality is correct in 1D; in higher dimensions, it would involve the density of states.  \nThis ground state, once computed, sets an important energy scale of the system: the energy of the Gardner distribution E G . At energies much larger than this (such as in strongly turbulent systems), one should expect the effect of phase-volume exclusion to become sub-dominant to the effect that pushes each level set towards a Maxwellianequivalent to the denominator of [9] being approximately unity. Then one should anticipate the solution to have the approximate form  \nP ( v , η > 0) ≈ e -βηε ( v ) F ( η ) . [12]  \n|  \nThe fugacity F ( η ) can then immediately be deduced from [7]:  \nF ( η ) ≈ ρ ( η ) [∫ d v e -βηε ( v ) ] -1 1D = √ βηm 2 π ρ ( η ) . [13]  \nThe mean phase-space density can now be computed from [12] and [13] via [2]: in 1D-1V,  \n〈 f 〉 ( v ) ≈ ∫ η max η min d η √ βηm 2 π ηρ ( η ) e -βηε ( v ) . [14]  \nThe meaning of equation [14] is physically transparent. Each level set relaxes to a Maxwellian distribution whose temperature is inversely proportional to its phase-space density: less dense portions behave as though they were lighter particles and more dense portions as though they were heavier. The resulting mean phase-space density is, therefore, a superposition of many Maxwellians with relative abundances set by ρ ( η ). At early times in our simulation, ρ ( η ) is approximately a delta function in η , as in [10] (by design). This delta function selects from the integral in [14] a single η -so the mean phase-space density becomes an actual Maxwellian. This is precisely what is seen in figures 2 and 3. We further note that this Maxwellianization is a fundamentally collisionless effect : it is an entropic property of collisionless plasma that turbulence should initially want to push mono-energetic beams towards a Maxwellian equilibrium. That this is a collisionless effect is also obvious as the mean phase-space density is inhomogeneous and continues to evolve after reaching a Maxwellian (the behaviour that would be forbidden by collisional dynamics).  \nAs the collisionless system evolves, ρ ( η ) changes form, as seen in figure 2, developing a lowη asymptotic for which ρ ( η ) ∼ η -1 . We may therefore assume that ρ ( η ) has a time-asymptotic limit  \nρ ( η ) → η -1 v b G ( η η max ) , [15]  \nwhere G is a dimensionless function that only has strong dependence on η near the maximum phase-space density η = η max . Then [14] can be written as  \n〈 f 〉 ( v ) = v b ε ( v ) 3 / 2 ∫ βη max ε ( v ) βη min ε ( v ) d¯ η √ m ¯ η 2 βπ G ( ¯ η βη max ε ( v ) ) e -¯ η , [16]  \nwhere we have changed the integration variable to ¯ η = βηε ( v ). Therefore, for velocities such that βη max ε ( v ) /greatermuch 1 and βη min ε ( v ) /lessmuch 1, the ¯ η integral will be a weak function of v , giving 〈 f 〉 a universal power-law tail. Recasting this asymptotic as a particle-energy distribution, we get  \nN ( ε ) ∝ ε -1 / 2 〈 f 〉 ∼ ε -2 . [17]  \nIt is easy to show that, for ρ ( η ) ∼ η -1 , this highenergy scaling of the particle-energy distribution is true in any number of dimensions (36). There is a subtlety associated with what the value of β turns out to be and, therefore, whether neglecting degeneracy effects was ever a good approximation. This is treated in detail in (36). Namely, it turns out that, while neglecting phase-space degeneracy is quantitatively incorrect, this does not invalidate the headline  \nFig. 4. Same as figure 1 but for the electron-positron two-stream instability, showing the distribution of electrons (which is statistically identical to that of the positrons). While the mixing of the distribution still occurs, the formation of large-scale turbulent structures present in figure 1 is entirely absent here, so the resulting distribution is much closer to its ground state. Details of the simulation setup and parameters are documented in the Materials and Methods.  \n<!-- image -->  \nx  \npe/  \nv  \nb  \nresult [17]. This occurs because the denominator of [9] is substantially different from unity only at low velocities where high-density level sets crowd each other out in the phase space, whereas at larger velocities, this effect is unimportantbut this is precisely where we expect the ε -2 power law to emerge.  \nUniversality of waterbag content. Thus, the universal highenergy tail N ( ε ) ∝ ε -2 is a direct consequence of the waterbag content of the system tending towards a ρ ( η ) ∼ η -1 asymptotic at low η -as indeed it does in our numerical experiment. While we do not know how to prove formally that this must happen, it can be justified qualitatively in the following way.  \nAs follows from [8], all Casimir invariants of the system can be recovered as moments of ρ ( η ). In particular, the zeroth moment of ρ ( η ) is the volume of phase space in which the phase-space density f takes a non-zero value. This will clearly be finite for our initial setup with two beams. As the turbulence stirs the plasma, f becomes highly filamented, developing ever sharper gradients in phase space. As collisions (equivalently, particle noise (41)) smooth out these sharp gradients of f , the volume in which f is nonzero should grow (this is manifest in figure 1). However, the first moment of ρ ( η ) (which is the particle number) must stay fixed, implying that ρ ( η ) must decrease at large η and increase at small η : most of the non-empty phase space must be occupied by relatively low phase-space densities . This is manifestly (although not uniquely) satisfied by ρ ( η ) ∼ η -1 .  \n|  \nAnother argument to the same effect is as follows. The energy of the Gardner distribution corresponding to ρ ( η ) (the ground state defined by [11]), can only increase under the action of collisional phase-space diffusion (see, e.g., (44)). It seems reasonable to conjecture, although not easy to prove rigorously, that f G should become more generic as it is thus heated. As was shown by (36), for a wide class of f G -all functions with any form of exponential decay at high energies § -have waterbag content with the asymptotic ρ ( η ) ∼ η -1 as η → 0.  \nDegenerate equilibria. While the formation of the ρ ( η ) ∼ η -1 asymptotic may be generic, it is not the only ingredient necessary to achieve the ε -2 power law. This is evidenced in our second numerical experiment: the (non-relativistic) electron-positron two-stream instability-a purely numerical invention (although a reality in the relativistic setup (45, 46)), but an interesting case study because it has identical linear physics to the electron-only two-stream instability (up to a rescaling of time) but exhibits a vastly different saturation scenario. In the lower right panel of figure 5, we see that the turbulence again pushes the waterbag content towards the asymptotic ρ ( η ) ∼ η -1 . However, the mean phase-space density does not form an ε -2 tail because the assumption of non-degeneracy [12] is completely violated. This can be seen in the upper right panel of figure 5. Initially, as the beams are two thin slivers in phase space, the Gardner restacking would amount to just placing the beams  \nFig. 5. Left panel: comparison of the Gardner distribution function of electrons [11] (dashed lines) computed from the measured waterbag content ρ ( η, t ) and the mean phase-space density of electrons (solid lines) at a range of times during the evolution of the electron-positron two-stream instability. Upper right panel: ratios of the Gardner energy E G to the total system energy E tot during the evolution of the electron-positron (dot dashed) and electron-only (dotted) two-stream instabilities. The degenerate and non-degenerate limits are highlighted in red ( E G = E tot ) and blue ( E G = 0 ), respectively. Lower right panel: evolution of the waterbag content ρ ( η ) (of electrons) over time, as in figure 2, but this time for the electron-positron two-stream instability. Thus, while the electron-positron system still achieves a universal waterbag content, its Gardner energy has grown so close to its actual energy that it has been frozen into its ground state.  \n<!-- image -->  \n<!-- image -->  \n3000  \nat v = 0. The resulting ground state would have an energy E G much smaller than the total kinetic energy E tot , making the system highly non-degenerate. However, as the turbulence breaks ρ ( η ) conservation, the energy of the underlying Gardner distribution grows, causing the system to become more degenerate. In the case of the electrononly two-stream instability, the formation and persistence of coherent structures (phase-space holes) seen in figure 1 causes the growth of the Gardner energy to saturate. In contrast, for the electron-positron instability, the holes that do form fail to merge into large-scale, large-amplitude structures (see figure 4). As a result, the Gardner energy increases to meet the system's total energy and the evolution freezes in a degenerate state (a phenomenon somewhat similar to 'incomplete relaxation' (47)). This again amounts, of course, to the system reaching its Lynden-Bell equilibrium, but in a much more trivial manner: the distribution simply becomes the Gardner distribution, as seen in the left panel of figure 5.", 'Breaking of Casimir invariants and the phase-space cascade': "Let us now turn our attention to the question of why ρ ( η ) is able to evolve relatively quickly (and how quickly) even in low-collisionality systems. This is because the thorough mixing in phase space, which is required for the system to reach the Lynden-Bell maximum-entropy state, generates progressively smaller scales in velocity and position space (as seen in figure 1) until, at sufficiently small scales, collisions, however small is their rate, act to smooth the distribution function, altering ρ ( η ). For electrostatic plasmas, this process  \nhas been shown theoretically (26) and numerically (35) to be able to be described as a turbulent cascade through phase space of a representative Casimir invariant (cf. (25, 48-50))  \nC 2 = ∫ d η η 2 ρ ( η ) = 1 L ∫ d x d v f 2 , [18]  \nwhere L is the domain size. The distribution of this invariant across spatial and velocity scales can be quantified by its spectrum C ( k, s ) = | f k,s | 2 , where  \nf k,s = 1 L ∫ d x d v e -ikx + isv f ( x, v ) [19]  \nis the Fourier transform of the phase-space density in position and velocity space. Since C 2 is a quadratic norm of the phase-space density, the contributions to it from the (spatial) mean [2] and perturbed δf = f -〈 f 〉 parts of f add:  \n/negationslash  \nC 2 = C 2 , 0 + δC 2 = ∫ d v 〈 f 〉 2 + 1 2 π ∑ k =0 ∫ d s C ( k, s ) . [20]  \nIt is not hard to see that the relaxation of the initially unstable state will lead to C 2 , 0 decreasing and, therefore, to δC 2 receiving the balance of the C 2 density. This results in an approximately constant flux of C 2 towards small scales (large k and s ). A cascade theory in the spirit of Kolmogorov (51) leads to the following asymptotic form of the spectrum (35):  \nC k,s ∝ { s -2 , k /lessmuch γs, k -2 , k /greatermuch γs, [21]  \nFig. 6. Left panel: cuts of C ( k, s ) as a function of the velocity wave number s (see [19]). This spectrum was taken at time t = 250 ω -1 pe during the evolution of the electron-only two-stream instability. The lines shown in the blue-to-red color palette are C ( k, s ) = | f k,s | 2 for different values of k (plotted in increments of 0 . 8 ω pe /v b ), representing the velocity-space structure of fluctuations of different scales (large scales: dark blue; small scales: dark red). The theoretically expected asymptotic scalings [21] are shown in dashed black lines. Right panel: the same spectrum C ( k, s ) , plotted in 2D to show the so-called 'critical-balance' line s ∼ γ -1 k , separating the two asymptotics in [21]. Physically, the velocity-space linear phase mixing dominates at s /lessmuch γk and the position-space nonlinear mixing of the phase-space density by the electric field dominates at s /greatermuch γk .  \n<!-- image -->  \n<!-- image -->  \nv  \ns  \nwhere γ is a typical shearing rate in phase space set by the amplitude of the electric field. Since the electric field E = -∂ϕ/∂x is determined, via Poisson's equation  \n∇ 2 ϕ = 4 πe (∫ d v f -n 0 ) , [22]  \nby the perturbed electron density (for the experiment with static ions), the spectrum of the electric fluctuations at small spatial scales can be determined from the s → 0 asymptotic of [21]:  \n| E k | 2 = 16 π 2 e 2 k 2 C k,s → 0 ∝ k -4 . [23]  \nBecause this spectrum is quite steep, the phase-space density, at any scale, is predominantly stirred by the electric field at the outer scale (53), which, for the two-stream instability, are of the order of (several times) the Debye length λ D ∼ v b /ω pe . It is because of this large-scale dominance that the shearing rate γ in [21] can be assumed to be a scale-independent constant (35). This also implies that the typical time scale at which f will be mixed all the way to the phase-space scale l c where collisions-equivalently, particle noise-start erasing the small-scale structure is (35)  \nτ c ∼ γ -1 ln λ D l c . [24]  \nEven without directly computing l c , it is clear that l c /lessmuch λ D , provided the noise floor is sufficiently low (equivalent to the plasma being weakly coupled, i.e., n 0 λ D /greatermuch 1). Thus, [24] tells us that ρ ( η ) will change on time scales that are only logarithmically longer than the dynamical time scale ( ∼ γ -1 )  \n|  \nof collisionless relaxation. This provides a modicum of justification for our scheme-confirmed by the numerical experiment with the electron two-stream instability-of evolving the mean phase-space density as an 'instantaneous' LyndenBell equilibrium coupled to a time-dependent waterbag content ρ ( η, t ).  \nThe above argument is supported by the excellent agreement between the theoretical predictions [21] and [23] and the spectra measured in our numerical simulation of the electron two-stream instability. The phase-space spectrum C ( k, s ) and the electric-energy spectrum | E k | 2 shown in figures 6 and 7, respectively, approach their theoretically predicted asymptotic forms around the time when the electron holes, vividly displayed in figure 1, begin to move around and merge. It is these dynamics that provide the vigorous mixing that leads to the phase-space cascade and ultimately pushes ρ ( η ) and, therefore, 〈 f 〉 ( v ) towards their universal forms discussed and measured above.", 'Discussion': "In this work, we have shown that statistical mechanics can be used to predict classes of universal equilibria for relaxing collisionless plasmas. We have tested this proposition on the example of an electrostatic plasma destabilized by one of the simplest, best-studied instabilities in plasma physics, the two-stream instability, and found good agreement with the theory.  \nThe theory is based on the thermodynamic approach first proposed by Lynden-Bell (17), maximising the entropy [3] subject to the conservation of phase volume [1]. The  \nFig. 7. The spectrum of the electric field for a range of times during the evolution of the electron-only two-stream instability. The lines shown in the blue-to-red color palette represent the times from t = 200 ω -1 pe (dark blue) to t = 3200 ω -1 pe (dark red), plotted in increments of 200 ω -1 pe . The dashed black lines represent the theoretical prediction of a phase-space cascade [23], | E k | 2 ∝ k -4 , which gives way to | E k | 2 ∝ k -2 , the floor due to the Poisson shot noise of discrete particles at large k (cf. (35, 52)).  \n<!-- image -->  \nadditional constraints arising from the latter can be captured by tracking the level sets of the phase-space density via the 'waterbag content' ρ ( η ) [7] and endowing the resulting equilibrium [9] with some memory of the plasma's earlier state. For the case of the electron two-stream instability, this theory predicts the formation of a Maxwellian distribution (despite the purely collisionless dynamics). This is indeed achieved at early times (see figures 2 and 3). However, unlike in a collisional regime, the system continues evolving after reaching a Maxwellian. This further evolution is driven by phase-space turbulence stirred up by electron holes generated in the early stages of the instability. We show that this turbulence drives a phase-space cascade of the form predicted by (26, 35), giving rise to small-scale structure in both velocity and position space (see figures 6 and 7). Since no system is truly collisionless (there is always a finite number of particles), this small-scale structure causes the collisionless invariant ρ ( η ) to be broken (see figures 2 and 5), causing the LyndenBell equilibrium of the system to change and the system to continue evolving on time scales that we estimate, in [24], to be only logarithmically longer than the dynamical relaxation times. At these later times, the waterbag content develops a lowη asymptotic ρ ( η ) ∼ η -1 , which we argue to be universal, being the natural ρ ( η ) associated with systems that have smooth ground-state phase-space densities. We show that this is indeed what happens in the numerical experiments featuring both the electron-only two-stream instability and the electron-positron two-stream instability. The details of the resulting equilibria in the presence of this universal waterbag content then depend strongly on the amount of  \n|  \nenergy the well-mixed system is able to retain relative to its ground state (defined by [11]). In the case of the electron-only two-stream instability, the saturated energy of the system is larger than the energy of the ground state by a factor of order unity and, as a result, the corresponding Lynden-Bell equilibrium [9] is approximately non-degenerate, featuring a distribution of particle energies that has a universal powerlaw tail ∝ ε -2 . In contrast, the case of electron-positron two-stream instability exemplifies systems where the nonconservation of ρ ( η ) causes the system to freeze in its ground state, achieving a fully degenerate Lynden-Bell equilibrium (cf. (54)).  \nThus, we have two examples representing what are likely to be two equivalence classes of universal collisionless equilibria: those that, given initial energy E tot , relax to (approximately) non-degenerate Lynden-Bell distributions such that the corresponding ground state's energy E G is a finite fraction of (or, better still, much smaller than) the system's energy E tot -and those for which E G ≈ E tot in the final state, which is, therefore, a fully degenerate Lynden-Bell equilibrium. It is the former class that features the universal high-energy tail ∝ ε -2 . What appears to help such a state into existence is the emergence of highly non-equilibrium structures-in the case of the electron two-stream instability, electron holesthat can store a certain amount of energy and engage in long-time nonlinear dynamics (moving around, merging) that stir up the plasma and keep it turbulent, rather than decaying into a ground state with no available energy. In a certain sense, this is similar to a non-equilibrium, driven system, where there is continuous injection of δC 2 (see [18]), and, presumably, other moments of ρ ( η ), into small scales-which is why the turbulent spectra that we have observed (figures 6 and 7) can be predicted by theories that assume continuous driving (26, 35).  \nThis argument contains an apparent internal contradiction: non-degenerate Lynden-Bell equilibria emerge thanks to dynamics that must clearly be inimical to the perfect satisfaction of the hypothesis of perfect mixing. Indeed, it seems unlikely that, with phase-space holes roaming the system, all parts of the phase space can be plausibly assumed equally accessible. It appears, however, that enough of it is accessible for the maximum-entropy principle to assert itself in a theoretically computable way, and that this partial accessibility is a compromise that allows the system to remain turbulent, and converge to a statistical state that is both interesting and has a modicum of universality.  \nGiven the propensity for distributions with power-law tails in energy to occur in a wide variety of collisionless plasmasincluding, but not limited to, the solar wind (e.g., (55-57)), the solar neighbourhood (e.g., (2)), and numerous numerical studies (e.g., (15, 58-60))-it is interesting to speculate whether the adjusted Lynden-Bell formalism proposed (and validated) here could be used to explain some of the observed equilibria (obviously, with due regard to how the plasmas 'entropic' desire to relax to Lynden-Bell states might compete with non-equilibrium processes such as particle sources and losses, energization by various structures, etc.). As has been abundantly clear since its inception, the strongest feature of the Lynden-Bell approach has also been its fundamental drawback: it is a thermodynamic theory and, while thermodynamic theories can be remarkably robust  \nand successful, it is also difficult to predict how badly their underlying assumptions must be broken for the theory to fail completely. The idea of evolving waterbag content ρ ( η ) offers some promise for understanding how the system should relax. There is, however, currently no theory of how to compute dynamically the time evolution of ρ ( η ), which is the final piece of the puzzle. Since it is this evolution that determines how far from its ground state the system saturates, it is this final piece that should teach us how to sort relaxing plasma systems between the two universality classes identified above.", 'Materials and Methods': "Here, we provide some details about the set-ups of our numerical simulations of the relaxation of two-stream instabilities to LyndenBell equilibria. The simulations were conducted using the PIC code OSIRIS (29) with an initial condition containing two thin beams:  \nf ( x , v ) = { n 0 2∆ v b , v b -∆ v b 2 < | v | < v b + ∆ v b 2 , 0 , otherwise , [25]  \nwhere v b is the beams' velocity and ∆ v b is their width. For the simulation of the two-stream electron-positron instability, the same initial condition was used for both species. As OSIRIS is a fully relativistic code, a beam velocity of v b = c/ 20 was chosen, so that the simulation was approximately classical with the relativistic gamma factor γ rel -1 ∼ 10 -3 . Both the simulation of the electron-only and electron-positron two-stream instabilities used a periodic 1D domain of size 25 . 6 d e (where d e is the electron skin depth) consisting of 2 15 cells. Both linear and quadratic interpolations for particles were checked, with no discernible difference. With these parameters, the energy and momentum were conserved up to diagnostic precision for the duration of the simulation. No smoothing was used on the fields or currents. For convergence tests, the simulation of the electron-only two-stream instability was repeated on a 1D periodic domain of size 210 d e with 268 , 800 cells, showing good convergence and much cleaner statistics. To our knowledge, this makes this the largest, longestrun, and most collisionless simulation of the electron-only twostream instability to date. In the main text, figures 1 and 6 use the results of the smaller electron-only simulation, while figures 2, 3, 5, and 7 use the larger simulation. Despite the fact that the peak growth rate of the two-stream instability sits at k ∼ √ 3 / 8 λ -1 D ∼ √ 3 / 8 d -1 e c/v b , the domain size was chosen as a trade off between the need for resolution of sub-Debye physics and the desire for the multiple hole mergers seen in figure 1 to drive vigorous relaxation of the distribution function towards equilibrium. The number of particles per cell was chosen to be 20 , 000 (in the case of the electron-positron simulation, 20 , 000 of each species), giving the effective plasma parameter of n 0 λ D ∼ 10 6 . Naively, this would imply that the time that it would take for collisions to modify the mean distribution function in the absence of turbulence should be ∼ 10 6 ω -1 pe . As a result, our simulations should be considered collisionless as far as the evolution of 〈 f 〉 is concerned, despite the fact that ρ ( η ) is modified on much shorter time scales, for the reasons explained in the main text.  \nIn order to compute the waterbag content ρ ( η ) of the exact phase-space density f , the PIC particles must be placed on a grid in ( x, v ) (since they are otherwise slivers with a finite extent in position space but zero width in velocity space). Therefore, we compute the fine-grained phase-space density f on evenly spaced rectangular bins with a width of 4 cells in position space and 0 . 0048 v b in velocity space in the interval between v = ± 0 . 25 c . The integral of the waterbag content  \nΓ( η ) = 1 L ∫ d x d v Θ( f -η ) = ∫ ∞ η d η ' ρ ( η ' ) [26]  \nis then computed assuming a piecewise-constant f on this finegrained grid. This quantity is computed on a grid of 1000 logarithmically spaced values of η between the maximum phasespace density η max (at that time step) and η min corresponding  \n|  \nto one particle per bin. We tested the results of this scheme for different bin sizes. The chosen bin size is a compromise between the need to have small bins in order to resolve fine-grained features of the phase-space density, and the need to have bins that are sufficiently large to capture many particles. Our results were insensitive to increasing or decreasing the bin dimensions by up to a factor of 16. Due to the shape of the PIC particles in position space, it is also possible to change the lower limit for the η grid (as a fine-grained cell may contain a fraction of a particle). The results are insensitive to changing this value by a factor of 10 in either direction.  \nFrom Γ( η ) and the system's energy, one knows everything required for the Lynden-Bell equilibrium to be computed numerically with correct β and F ( η ). This is done using the iterative scheme laid out in (36) ¶ . This iteration continues until a Lynden-Bell equilibrium is found for which the normalized root-mean-square error in ρ ( η ) is below 10 -4 and the error in the energy (normalized to the Gardner energy) is below 10 -3 .  \nACKNOWLEDGMENTS. It is a pleasure to acknowledge encouraging discussions with Barry Ginat, Chris Hamilton, David Hosking, and Anatoly Spitkovsky (to whom we are especially grateful for suggesting an electron-positron two-beam experiment). We thank FCCN/RNCA (Portugal) for access to MareNostrum 5 (BSC, Spain) for computational work. RJE was supported by a UK EPSRC studentship; MLN by a Clarendon Scholarship; TS by FCT, Portugal (IPFN-CEEC-INST-LA3); PJB also by FCT, Portugal (Project X-maser No.2022.02230.PTDC and Grant No. UI/BD/151559/2021); the work of AAS was supported in part by grants from STFC (ST/W000903/1) and EPSRC (EP/R034737/1), as well as by the Simons Foundation via a Simons Investigator award. This research was also supported in part by the NSF grant PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP) and benefited from many vigorous interactions with the participants of their programme on 'Interconnections between the Physics of Plasmas and Self-gravitating Systems'.  \n- 1. J Becker Tjus, L Merten, Closing in on the origin of galactic cosmic rays using multimessenger information. Phys. Rep. 872 , 1 (2020).\n- 2. M Oka, et al., Electron power-law spectra in solar and space plasmas. Space Sci. Rev. 214 (2018).\n- 3. LB Wilson, et al., The need for accurate measurements of thermal velocity distribution functions in the solar wind. Front. Astron. Space Sci. 9 , 369 (2022).\n- 4. F Cruz, et al., Electron acceleration by wave turbulence in a magnetized plasma. Nat. Phys. 14 , 475 (2018).\n- 5. RM Magee, et al., Direct observation of ion acceleration from a beam-driven wave in a magnetic fusion experiment. Nat. Phys. 15 , 281 (2019).\n- 6. E Hartouni, et al., Evidence for suprathermal ion distribution in burning plasmas. Nat. Phys. 19 , 1 (2022).\n- 7. L Boltzmann, Vorlesugnen über Gastheorie . (Leipzig: J. A. Barth), (1896).\n- 8. LD Landau, Transport equation in the case of Coulomb interaction. Zh. Eksp. Teor. Fiz. 7 , 203 (1936).\n- 9. R Balescu, Irreversible processes in ionized gases. Phys. Fluids 3 , 52 (1960).\n- 10. A Lenard, On Bogoliubov's kinetic equation for a spatially homogeneous plasma. Ann. Phys. 10 , 390 (1960).\n- 11. NA Krall, AW Trivelpiece, Principles of Plasma Physics . (McGrawHill, New York), (1973).\n- 12. AFA Bott, SC Cowley, AA Schekochihin, Kinetic stability of Chapman-Enskog plasmas. J. Plasma Phys. 90 , 975900207 (2024).\n- 13. E Marsch, Kinetic physics of the solar corona and solar wind. Living Rev. Sol. Phys. 3 , 1 (2006).\n- 14. V Pierrard, M Lazar, Kappa distributions: theory and applications in space plasmas. Sol. Phys. 267 , 153 (2010).\n- 15. V Zhdankin, GR Werner, DA Uzdensky, MC Begelman, Kinetic turbulence in relativistic plasma: from thermal bath to nonthermal continuum. Phys. Rev. Lett. 118 , 055103 (2017).\n- 16. D Verscharen, KG Klein, BA Maruca, The multi-scale natrue of the solar wind. Living Rev. Sol. Phys. 16 , 5 (2019).\n- 17. D Lynden-Bell, Statistical mechanics of violent relaxation in stellar systems. Mon. Not. R. Astron. Soc. 136 , 101 (1967).\n- 18. H Ye, PJ Morrison, Action principles for the Vlasov equation. Phys. Fluids B 4 , 771 (1992).\n- 19. I Arad, PH Johansson, A numerical comparison of theories of violent relaxation. Mon. Not. Roy. Astron. Soc. 362 , 252 (2005).\n- 20. I Arad, D Lynden-Bell, Inconsistency in theories of violent relaxation. Mon. Not. Roy. Astron. Soc. 361 , 385 (2005).\n- 21. Y Levin, R Pakter, TN Teles, Collisionless relaxation in non-neutral plasmas. Phys. Rev. Lett. 100 , 040604 (2008).  \n- 22. M Assllani, D Fanelli, A Turchi, T Carletti, X Leoncini, Statistical theory of quasistationary states beyond the single water-bag case study. Phys. Rev. E 85 , 021148 (2012).\n- 23. Y Levin, R Pakter, FB Rizzato, TN Teles, FP Benetti, Nonequilibrium statistical mechanics of systems with long-range interactions. Phys. Rep. 535 , 1 (2014).\n- 24. AA Schekochihin, et al., Gyrokinetic turbulence: a nonlinear route to dissipation through phase space. Plasma Phys. Control. Fusion 50 , 124024 (2008).\n- 25. GL Eyink, Cascades and dissipative anomalies in nearly collisionless plasma turbulence. Phys. Rev. X 8 , 041020 (2018).\n- 26. ML Nastac, et al., Phase-space entropy cascade and irreversibility of stochastic heating in nearly collisionless plasma turbulence. Phys. Rev. E 109 , 065210 (2024).\n- 27. T Tatsuno, et al., Nonlinear phase mixing and phase-space cascade of entropy in gyrokinetic plasma turbulence. Phys. Rev. Lett. 103 , 015003 (2009).\n- 28. V Zhdankin, Generalized entropy production in collisionless plasma flows and turbulence. Phys. Rev. X 12 , 031011 (2022).\n- 29. RA Fonseca, et al., OSIRIS: A Three-Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based Accelerators. in Computational Science - ICCS 2002 . (Springer, Berlin-Heidelberg), p. 342 (2002).\n- 30. KV Roberts, HL Berk, Nonlinear evolution of a two-stream instability. Phys. Rev. Lett. 19 , 297 (1967).\n- 31. RL Morse, CW Nielson, One-, two-, and three-dimensional numerical simulation of two-beam plasmas. Phys. Rev. Lett. 23 , 1087 (1969).\n- 32. IB Bernstein, JM Greene, MD Kruskal, Exact nonlinear plasma oscillations. Phys. Rev. 108 , 546 (1957).\n- 33. IH Hutchinson, Electron holes in phase space: what they are and why they matter. Phys. Plasmas 24 , 055601 (2017).\n- 34. IH Hutchinson, Kinetic solitary electrostatic structures in collisionless plasma: phase-space holes. arXiv , 2407.08539 (2024).\n- 35. ML Nastac, RJ Ewart, J Juno, M Barnes, AA Schekochihin, Universal fluctuation spectrum in Vlasov-Poisson turbulence. in preparation (2024).\n- 36. RJ Ewart, ML Nastac, AA Schekochihin, Non-thermal particle acceleration and power-law tails via relaxation to universal Lynden-Bell equilibria. J. Plasma Phys. 89 , 905890516 (2023).\n- 37. PH Chavanis, J Sommeria, R Robert, Statistical mechanics of two-dimensional vortices and collisionless stellar systems. Astrophys. J. 471 , 385 (1996).\n- 38. RJ Ewart, A Brown, T Adkins, AA Schekochihin, Collisionless relaxation of a Lynden-Bell plasma. J. Plasma Phys. 88 , 925880501 (2022).\n- 39. PH Chavanis, Coarse-grained distributions and superstatistics. Phys. A 359 , 177 (2006).\n- 40. C Birdsall, A Langdon, Plasma Physics via Computer Simulation . (CRC Press, United Kingdom), (1985).\n- 41. M Touati, et al., Kinetic theory of particle-in-cell simulation plasma and the ensemble averaging technique. Plasma Phys. Control. Fusion 64 , 115014 (2022).\n- 42. U Banik, A Bhattacharjee, W Sengupta, Universal non-thermal power-law distribution functions from the self-consistent evolution of collisionless electrostatic plasmas. arXiv , 2408.07127 (2024).\n- 43. CS Gardner, Bound on the energy available from a plasma. Phys. Fluids 6 , 839 (1963).\n- 44. S Tremaine, M Henon, D Lynden-Bell, H-functions and mixing in violent relaxation. Mon. Not. Roy. Astron. Soc. 219 , 285 (1986).\n- 45. CD Arrowsmith, et al., Laboratory realization of relativistic pair-plasma beams. Nat. Commun. 15 , 5029 (2024).\n- 46. K Qu, NJ Fisch, Creating and detecting observable QED plasmas through beam-driven cascade. Phys. Plasmas 31 , 062102 (2024).\n- 47. PH Chavanis, Quasi-stationary states and incomplete violent relaxation in systems with long-range interactions. Phys. A 365 , 102 (2006).\n- 48. G Knorr, Time asymptotic statistics of the Vlasov equation. J. Plasma Phys. 17 , 533 (1977).\n- 49. PH Diamond, SI Itoh, K Itoh, Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas . (Cambridge University Press, Cambridge), (2010).\n- 50. S Servidio, et al., Magnetospheric multiscale observation of plasma velocity-space cascade: Hermite representation and theory. Phys. Rev. Lett. 119 , 205101 (2017).\n- 51. A Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds' numbers. Doklady Akad. Nauk SSSR 30 , 301 (1941).\n- 52. N Rostoker, Fluctuations of a plasma (I). Nucl. Fusion 1 , 101 (1961).\n- 53. GK Batchelor, Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity. J. Fluid Mech. 5 , 113 (1959).\n- 54. DN Hosking, D Wasserman, SC Cowley, Metastability of stratified magnetohydrostatic equilibria and their relaxation. arXiv , 2401.01336 (2024).\n- 55. G Gloeckler, LA Fisk, GM Mason, ME Hill, Formation of power law tail with spectral index -5 inside and beyond the heliosphere. AIP Conf. Proc. 1039 , 367 (2008).\n- 56. LA Fisk, G Gloeckler, The case for a common spectrum of particles accelerated in the heliosphere: Observations and theory. J. Geophys. Res. Space Phys. 119 , 8733 (2014).\n- 57. L Yang, et al., Quiet-time solar wind suprathermal electrons of different solar origins. Astrophys. J. Lett. 896 , L5 (2020).\n- 58. A Spitkovsky, Particle acceleration in relativistic collisionless shocks: Fermi process at last? Astrophys. J. 682 , L5 (2008).\n- 59. L Sironi, A Spitkovsky, Relativistic reconnection: an efficient source of non-thermal particles. Astrophys. J. 783 , L21 (2014).\n- 60. GR Werner, DA Uzdensky, Nonthermal particle acceleration in 3D relativistic magnetic reconnection in pair plasma. Astrophys. J. Lett. 843 , L27 (2017)."}