Phase-space entropy cascade and irreversibility of stochastic heating in nearly collisionless plasma turbulence

We consider a nearly collisionless plasma consisting of a species of `test particles' in 1D-1V, stirred by an externally imposed stochastic electric field. The mean effect on the particle distribution function is stochastic heating. Accompanying this heating is the generation of fine-scale structure in the distribution function, which we characterize with the collisionless (Casimir) invariant $C_2 \propto \iint dx dv \, \langle f^2 \rangle$. We find that $C_2$ is transferred from large scales to small scales in both position and velocity space via a phase-space cascade enabled by both particle streaming and nonlinear interactions between particles and the stochastic electric field. We compute the steady-state fluxes and spectrum of $C_2$ in Fourier space, with $k$ and $s$ denoting spatial and velocity wavenumbers, respectively. Whereas even the linear phase mixing alone would lead to a constant flux of $C_2$ to high $s$ (towards the collisional dissipation range) at every $k$, the nonlinearity accelerates this cascade by intertwining velocity and position space so that the flux of $C_2$ is to both high $k$ and high $s$ simultaneously. Integrating over velocity (spatial) wavenumbers, the $k$-space ($s$-space) flux of $C_2$ is constant down to a dissipation length (velocity) scale that tends to zero as the collision frequency does, even though the rate of collisional dissipation remains finite. The resulting spectrum in the inertial range is a self-similar function in the $(k,s)$ plane, with power-law asymptotics at large $k$ and $s$. We argue that stochastic heating is made irreversible by this entropy cascade and that, while collisional dissipation accessed via phase mixing occurs only at small spatial scales rather than at every scale as it would in a linear system, the cascade makes phase mixing even more effective overall in the nonlinear regime than in the linear one.

We consider a nearly collisionless plasma consisting of a species of 'test particles' in one spatial and one velocity dimension, stirred by an externally imposed stochastic electric field-a kinetic analogue of the Kraichnan model of passive advection.The mean effect on the particle distribution function is turbulent diffusion in velocity space-known as stochastic heating.Accompanying this heating is the generation of fine-scale structure in the distribution function, which we characterize with the collisionless (Casimir) invariant C2 ∝ dxdv ⟨f 2 ⟩-a quantity that here plays the role of (negative) entropy of the distribution function.We find that C2 is transferred from large scales to small scales in both position and velocity space via a phase-space cascade enabled by both particle streaming and nonlinear interactions between particles and the stochastic electric field.We compute the steady-state fluxes and spectrum of C2 in Fourier space, with k and s denoting spatial and velocity wavenumbers, respectively.In our model, the nonlinearity in the evolution equation for the spectrum turns into a fractional Laplacian operator in k space, leading to anomalous diffusion.Whereas even the linear phase mixing alone would lead to a constant flux of C2 to high s (towards the collisional dissipation range) at every k, the nonlinearity accelerates this cascade by intertwining velocity and position space so that the flux of C2 is to both high k and high s simultaneously.Integrating over velocity (spatial) wavenumbers, the k-space (s-space) flux of C2 is constant down to a dissipation length (velocity) scale that tends to zero as the collision frequency does, even though the rate of collisional dissipation remains finite.The resulting spectrum in the inertial range is a self-similar function in the (k, s) plane, with power-law asymptotics at large k and s.Our model is fully analytically solvable, but the asymptotic scalings of the spectrum can also be found via a simple phenomenological theory whose key assumption is that the cascade is governed by a 'critical balance' in phase space between the linear and nonlinear time scales.We argue that stochastic heating is made irreversible by this entropy cascade and that, while collisional dissipation accessed via phase mixing occurs only at small spatial scales rather than at every scale as it would in a linear system, the cascade makes phase mixing even more effective overall in the nonlinear regime than in the linear one.

I. INTRODUCTION
Understanding the nature of turbulent cascades in nearly collisionless space and astrophysical plasmas is an outstanding problem [1][2][3][4][5][6][7] with a diverse range of applications, from solving the coronal-heating problem [7] to interpreting radiation emission from accretion disks around black holes [8][9][10].A major difficulty distinguishing such turbulence from its fluid counterparts lies in the fact that fluctuations evolve in the six-dimensional phase space of single-particle positions and velocities.Dissipation (in the sense of irreversible entropy production) and any transfer of energy between particles and fields is formally reversible.
It was realized by [22,23] that entropy production in the long-time limit of a nearly collisionless plasma must remain finite even as the collision frequency ν tends to zero.This idea crystallized in [2,[24][25][26], where the notion of entropy cascade in the context of gyrokinetics was introduced.Below the Larmor scale, in a phase-space inertial range between the injection range at large scales and the collisional dissipation range at small scales, the (negative) entropy of the perturbed distribution function cascades in both position and velocity space via a nonlinear perpendicular phase-mixing mechanism [27].Because this is a constant-flux cascade, the turbulent heating rate in a gyrokinetic plasma is finite and independent of ν even as ν → 0 + , analogous to so-called dissipative anomalies [28,29] in hydrodynamics, where viscous dissipation is finite in infinite-Reynolds-number turbulence.
Entropy cascade outside the gyrokinetic approximation is a frontier topic just beginning to be explored [30][31][32][33][34].In this paper, we study turbulent dissipation and entropy cascade via an analytically solvable model introduced by [35].We consider the '1D-1V' electrostatic, full-f Vlasov equation with a model diffusive collision operator for a test-particle species, where instead of the electric field being self-consistently determined by Poisson's equation, it is externally determined to be a stochastic Gaussian, white-noise source with a specified spatial correlation function.Physically, this model represents the evolution of a low-density minority species in a multi-component plasma whose dielectric response is dominated by the other, more abundant species.This model is the plasma-kinetic analogue of the Kraichnan [36] model of passive advection, where a scalar field [37], such as temperature or concentration of dye, is passively advected by an externally determined random flow, allowing for analytical calculations of the passive-scalar statistics.The Kraichnan model has been dubbed the 'Ising model' [38] of fluid turbulence because it is a solvable model that exhibits many properties also present in real systems, and so serves as a theoretical laboratory to study turbulence [39].It is in this spirit that we investigate our solvable model of kinetic plasma turbulence.
We decompose the distribution function into its mean and fluctuating parts, f = ⟨f ⟩+δf , where ⟨...⟩ denotes ensemble averaging over realizations of the random electric field.In our model, particles are stochastically accelerated by the electric field and undergo random walks in velocity space, resulting in bulk heating of ⟨f ⟩ [40].In the context of magnetized plasmas, this phenomenon is often referred to as stochastic heating [6,41,42].Accompanying this heating is the generation of velocity and spatial structure in the perturbed distribution function δf via linear phase mixing and nonlinear interactions between particles and the turbulent electric field [35].
We characterize this structure via the collisionless (Casimir) invariant C 2 = (1/L) dxdv ⟨f 2 ⟩/2, where L is the system size.C 2 has been considered before as a measure of phase-space structure and as a cascaded quantity in kinetic plasma turbulence [30,32,35,[43][44][45][46][47], and is closely related to the part of the traditional entropy, S = − dxdv f log f , associated with the perturbed distribution function and entering additively in the freeenergy invariant of δf gyrokinetics [2,22,24,48,49].The conservation of C 2 is broken by particle collisions, and, in particular, when collisions are modeled as a linear diffusion operator in velocity space, as they are in this paper, the collisional dissipation of C 2 is negativedefinite.Because the time irreversibility of our system can be tracked via non-conservation of C 2 , it can be used as a generalized (negative) entropy of the distribution function.
The diffusion of ⟨f ⟩ by the electric field has the side effect of injecting δC 2 = (1/L) dxdv ⟨δf 2 ⟩/2 fluctuations at large scales.These are then cascaded to small scales in both position and velocity space, where they are ultimately dissipated by collisions.This phase-space cascade of δC 2 is due to both linear phase mixing and nonlinear interactions between particles and the stochastic electric field.We analyze this cascade by computing the steadystate Fourier spectrum and fluxes of δC 2 in both position and velocity space, with dual variables (wavenumbers) k and s, respectively.
In the absence of nonlinearity, linear phase mixing advects the spectrum from low to high |s|, giving rise to an 'inertial range' in s where there is a constant flux of δC 2 from injection to dissipation scales, at every k [50,51].The resulting steady-state spectrum is flat in the inertial range, with an exponential cutoff at a ν-dependent collisional scale s ν .
Under the Kraichnan model, we find that the nonlinear term in the evolution equation for the Fourier spectrum becomes a fractional Laplacian operator [52][53][54] in k space, which leads to anomalous diffusion [55,56].Whereas fractional Laplacians (and fractional derivatives in general) are usually introduced in an ad-hoc manner to model systems with anomalous diffusion, both in plasma physics [57][58][59] and a wide variety of other contexts [60], here it emerges naturally as a result of our assumptions about the electric field.
The linear phase mixing and the turbulent fractional diffusion intertwine the position-and velocity-space cascades in such a way that the resulting spectrum is a self-similar function in the (k, s) plane, with power-law asymptotics at large k and s.Even though the Kraichnan model is fully analytically solvable, we can also recover these asymptotic scalings via a phenomenological theory whose key assumption is that the cascade is governed by a 'critical balance' in phase space between the linear and nonlinear time scales.
The δC 2 flux has components in both k and s directions.The flow of δC 2 occurs along outward unwinding spirals in (k, s) space.This circuitous route to dissipation scales is due to the nonlinearity generating modes that can linearly propagate from high to low |s|, called 'phase-unmixing' modes [61], a stochastic generalization of the textbook phenomenon of plasma echo [62,63].The net result after adding together contributions to the flux from the phase-mixing and phase-unmixing modes is that δC 2 is cascaded to both high s and high k simultaneously, and collisional dissipation only occurs at scales comparable to, or smaller than, the dissipation length and velocity scales, both of which tend to zero when the collision frequency ν does.Integrating over velocity (spatial) wavenumbers, the flux of δC 2 in k (s) space is constant down to these dissipation scales, beyond which perturbations are thermalized by collisions.The rate of collisional dissipation is finite and independent of the collision frequency as ν → 0 + -a clear analytical example of a 'dissipative anomaly' in a kinetic system.This turbulent dissipation ultimately mediates the irreversibility of the stochastic heating.
The rest of this paper is organized as follows.In Section II, we introduce the Vlasov-Kraichnan model.In Section III, we construct a phenomenological theory that captures the asymptotic scalings of the Fourier spectrum of δC 2 .Then, in Section IV, we directly calculate the 1-D phase-space fluxes of δC 2 in Fourier space in a statistical steady state.In Section V, we calculate the inertialrange Fourier spectrum and its corresponding 2-D fluxes in (k, s) space.Finally, in Section VI, we conclude our results and discuss their implications.Supplementary calculations and discussions are exiled to Appendices A-E.

II. KRAICHNAN MODEL FOR A 1D-1V
ELECTROSTATIC PLASMA We consider a test-particle species composed of particles with charge q and mass m in a 1D periodic box of length L, and subject to an external electric field E. At t = 0, we assume the particle distribution function f to have no spatial variation, but we keep its velocity dependence generic, only assuming that f is square-integrable and has finite kinetic energy.We denote the number density of the distribution function as n 0 and the initial thermal velocity as v th,0 = 2T 0 /m, where T 0 is the initial temperature of the particles.An example initial condition with these properties is a Maxwellian, The Vlasov equation for the particle distribution function is where C[f ] is the collision operator, and we have absorbed q and m into the definition of E, denoting qE/m → E.
We assume E to be a Gaussian white-noise field, with zero mean and correlation function where ⟨...⟩ denotes ensemble averaging over realizations of E and δ is a Dirac delta distribution.We assume E to be statistically homogeneous and isotropic in space, so D(x, x ′ ) = D(r), where r = |x − x ′ |.We choose where k ∈ (2π/L)Z and Here, D is a constant diffusion coefficient with dimensions (length 1−α ) × (time −3 ), 0 < α ≤ 2, and L E and η represent the integral length scale and dissipation scale, respectively, of the stochastic electric field.The electric field is chosen so that its correlation function (3) has a power-law spectrum ∝ |k| −(α+1) in the inertial range 1/L E ≪ k ≪ 1/η; note that ( 5) is defined in a similar way to the velocity field in the fluid passive scalar Kraichnan model [39,64].We identify two distinct regimes: α < 2 and α = 2.When α < 2, the field is multiscale, reminiscent of turbulent fields in fully developed turbulence.When α = 2, the spectrum is sufficiently steep that the field is effectively single-scale.This case is known as the Batchelor regime [65].While there are important differences between the two regimes [39], some of which we will discuss, many of the properties of the model considered in this paper will be qualitatively the same in both regimes.
It will be useful to decompose the distribution function into its mean and fluctuating parts: We make no assumption of δf being small compared to ⟨f ⟩.
The effect of collisions will be to wipe out fine-scale velocity-space structure in the distribution function.To model this in the simplest possible way, we ignore collisions between our test-particle species and the other species in the plasma, and represent collisions within the test species as a linear diffusion in velocity space acting only on δf , viz., where ν is the collision frequency (multiplied by v 2 th,0 ), which we consider to be vanishingly small, taking ν → 0 + .It is not a problem that (7) does not conserve energy or vanish on a Maxwellian because collisions will only matter for the parts of δf with sharp gradients in v [66].

A. Stochastic heating
We first work out the effect of the turbulent electric field on the mean distribution function.Ensemble averaging (2) over realizations of the stochastic electric field, we get The equation for δf is then To compute the ensemble average of the nonlinear term in (8), we apply the Furutsu-Novikov theorem [67,68] for splitting correlators.For a Gaussian field E that depends on variables q, this theorem states that where F [E] is a differentiable functional of E. Formally integrating (9) with respect to time, we have that Combining (10) and (11), we have Therefore, (8) becomes where D 0 = D(0) is the 'turbulent collisionality.'We have dropped the streaming term in (8) because our initial condition is spatially homogeneous, so ⟨f ⟩ at future times does not depend on x.The solution to (13) is The mean kinetic-energy density ⟨K⟩ of this distribution function is where K 0 is the initial kinetic-energy density.For the Maxwellian initial condition (1), ( 14) becomes simply with a growing thermal speed: As particles get stochastically accelerated by the electric field, they undergo Brownian random walks in velocity space, leading to bulk heating of the distribution function [40], with secularly growing kinetic energy (15).In magnetized plasmas, this phenomenon is usually referred to as stochastic heating [6,41,42].

B. C2 budget: injection and dissipation
Because the stochastic electric field continuously heats the distribution function, viz., (15), energy is not a conserved quantity in our system.However, what is conserved in the absence of collisions is the quadratic quantity where C 2 can only change via collisions.Using (2) and ( 7), we have Thus, when collisions are approximated as a linear diffusion in velocity space, they provide negative-definite dissipation of C 2 .Because −C 2 is conserved in the absence of collisions and positive-definitely increased by collisions, we can interpret it as a 'generalized entropy' of the distribution function.We will henceforth refer to −C 2 and entropy synonymously.Stochastic heating, which is a collisionless process, is accompanied by the decrease of C 2,0 .Indeed, using (13) and integrating by parts, we have generically, and for a Maxwellian in particular, where T = mv 2 th /2, with v th given by (17).To work out the δC 2 budget, we can combine (18) and (21), giving (24) Thus, without collisions, as C 2,0 decreases as a result of the stochastic heating of ⟨f ⟩, δC 2 increases to maintain entropy balance.Once δf has developed sharp enough gradients, collisions dissipate δf , increasing the total entropy.The irreversiblity of stochastic heating therefore hinges on the collisional dissipation of δf .
If the δC 2 fluctuations evolve faster than the mean C 2,0 and reach a quasi-steady state (as we shall argue that they do), then (24) becomes In the case where the initial condition is Maxwellian, combining ( 23) and ( 22) and substituting this expression into (25) yields a direct balance between the heating rate of ⟨f ⟩ and the collisional dissipation of δf .The velocity-space gradients of ⟨f ⟩ inject δC 2 at large scales, and collisions dissipate δC 2 at small scales.As in any prototypical turbulent system, the steady-state balance (25) between injection and dissipation at such disparate scales can hold if there is a constant-flux cascade bridging them.This is precisely what we will find in the following sections, viz., a phase-space cascade of δC 2 in both position and velocity space.This cascade will be our focus for the rest of this paper.
Before continuing, we note that in Appendix A, we discuss alternative formulations of the thermodynamics of our system in terms of other collisionless invariants, and why in this work we have chosen to study C 2 over those other invariants.

C. Evolution of Fourier spectrum
We analyze how δC 2 is partitioned between scales in phase space via its Fourier spectrum.We define and where k and s are dual variables to x and v, respectively.Fourier transforming (9), we get that δ f (k, s) satisfies We define the Fourier spectrum as which satisfies Parseval's theorem, has the budget equation (equivalent to (24) in spectral space) and satisfies the evolution equation Equation ( 33) is important, and its analysis will be the focus of the rest of this paper.There are several steps required to derive (32) and ( 33) from ( 29), primarily using (10) to calculate correlation functions involving the electric field.These steps are technical, so we detail them in Appendix B. The source term in (32) and ( 33) is which for a Maxwellian initial condition evaluates to where v th is given by (17).
In (33), (−∆ k ) α/2 is a fractional Laplacian [52][53][54] of order α/2 (in k space), and κ is a turbulent diffusion coefficient ∝ LD, for which the exact expression is given in Appendix B. To obtain this term, we have taken the limit η → 0 + in (5), which is analogous to taking the zeroviscosity limit in hydrodynamic turbulence, and kL E ≫ 1, which is a convenient limit to take to focus on how the distribution function is stirred in the 'inertial range' of the stochastic electric field.
The fractional Laplacian is a non-local integral operator, generalizing the Laplacian to non-integer order.Its Fourier transform satisfies where Note that we have converted sums over k into integrals (multiplied by the inverse step size L/2π).
Mathematically, the fractional Laplacian describes abstract 'particles' with coordinates (k, s) undergoing random jumps in k space, so-called Lévy flights, which leads to superdiffusion [55,56].In the limit α → 2 − , the fractional Laplacian becomes (minus) a regular Laplacian [52], corresponding to regular Brownian motion.This limit, the Batchelor regime, was studied before in [35], although we shall present some new conclusions about it below.
In the following sections, we will solve for steady-state solutions of (33).To this end, we further assume that the source is localized at small (k, s) and injects δC 2 at a constant rate, These assumptions may appear questionable in view of ( 24) and ( 22) (is the source really constant in time?) and of ( 34) and ( 5) (is it really local?).Before continuing, we address these concerns.
Regarding the constancy of the source, we observe that in (24), the injection rate of δC 2 via the stochastic heating of ⟨f ⟩, using (22) and dimensional analysis, scales as This scaling is particularly obvious, from (23), for the case of a Maxwellian initial condition.Comparing (39) and (38), we see that the injection rate is, in fact, timedependent, ∝ t −3/2 , even though in our calculations we would like to treat ε as a constant.However, if the cascade of δC 2 is set up quickly compared to the decay of the source (which, as will be discussed in Section VI B, it is on a ν-independent time scale when α < 2 and a time scale ∝ | log ν| when α = 2), it is reasonable to treat the source as approximately constant on the time scale of the δf evolution, and then solve for the steadystate spectrum.This becomes an ever-better approximation at long times, because the rate of change of ( 39) is a decreasing function of time.Also note that, as long as ε > 0, its actual value is not important; because the spectrum-evolution equation (B8) is linear, it has no special amplitude scale, and so the size of ε is just a global modifier of the spectrum's amplitude.
Regarding the locality of the source, it is easily satisfied in s space, as the velocity derivatives of ⟨f ⟩ become ever smaller over time as the distribution function stochastically heats.This is clearly seen in the case of a Maxwellian initial condition, where the source ( 35) is a Gaussian in s, with a width ∼ 1/v th .In k space, the source is not, in fact, truly localized, since D(k) is a power law, viz., (5), so fluctuations are injected at all scales where the electric field exists.However, this turns out not to be fatal to our formalism: we will neglect the source in our solving for the inertial-range spectrum in Sections III and V and argue a posteriori in Section V B that this choice was justified.

III. PHENOMENOLOGICAL THEORY OF PHASE-SPACE CASCADE
In Sections IV and V, we perform a detailed analysis of (33).However, the key results of those sections can, in fact, be intuited via a much simpler route, which we follow here first.
As discussed in Section II B, stochastic heating can only truly be made irreversible by particle collisions.For collisions to be relevant even in the limit ν → 0 + , we conjecture that δC 2 undergoes a cascade in both position and velocity space.Analogously to phenomenological theories of hydrodynamic turbulence [28,69], we assume that there exists an inertial range in phase space, whereby the flux ε is processed between scales, from the injection to dissipation range.Under this assumption, our task is then to ascertain the form of the spectrum F in (k, s) space in the inertial range.
Fluctuations of δf develop fine-scale structure in velocity space via linear phase mixing, which manifests in (33) as advection of F in s at the rate k.Dimensionally, the phase-mixing time is, therefore, Likewise, δf develops fine-scale structure in position space via nonlinear mixing by the stochastic electric field, which manifests in (33) as fractional diffusion of F in k with the diffusion coefficient κs 2 .Dimensionally, using (36), the turbulent-diffusion time is, therefore, The ratio of these time scales (taken to the power 1/(α + 1) for analytical convenience and in anticipation of the results of Section V) is We conjecture that the structure of F in (k, s) space is governed by the parameter ξ.In particular, we assume that the spectrum is a product of power laws in k and s, with different scaling exponents depending on whether ξ is small or large: When ξ ≪ 1, the turbulent-diffusion time is much shorter than the phase-mixing time, so we expect the inertial-range spectrum to satisfy, to lowest order in ξ, and therefore, to be independent of k.Consequently, a = 0 in (43).When ξ ≫ 1, the phase-mixing time is much shorter than the turbulent-diffusion time.Then, to lowest order in 1/ξ, we expect the spectrum to satisfy and, therefore, to be independent of s, the same as in the linear regime [50,51].Consequently, d = 0 in (43).
To find b and c, we invoke our initial assumption of a constant-flux cascade.As in any Kolmogorov-style theory, we need a prescription for the cascade time τ c , which is the typical time for δC 2 to be transferred across phasespace scales (ℓ, u), where ℓ ∼ 1/k and u ∼ 1/s.We conjecture that the cascade time is set by the phase-mixing and turbulent-diffusion time scales ( 40) and ( 41), and that the latter two must balance along the path of the cascade: In wavenumber space, this condition is satisfied when ξ ∼ 1, i.e., This is a kinetic, phase-space analogue of the criticalbalance conjecture in magnetohydrodynamic turbulence [70,71].Assuming a constant flux of δC 2 in position space, using the ℓ scaling in (46), and using (31) to relate fluctuations in real space and wavenumber space, we get where δf ℓ is the characteristic amplitude of δf at spatial scale ℓ.Let us compare this result to the 1-D k spectrum that follows from (43).We assume (and verify a posteriori ) that b > 1, so that the integral of F over s is dominated by the region ξ ≳ 1, i.e., This gives F (kL0, sv th ) vs. (kL0, sv th ), for α = 2, using the piecewise scalings (54).We use normalized units sv th and kL0, where L0 is defined in (91) (see the discussion at the beginning of Section V B).The black line is the critical-balance line (47), sv th ∼ kL0.In the τp ≪ τ d (ξ ≫ 1) region, the spectrum scales as k −2 , and in the τ d ≪ τp (ξ ≪ 1) region, the spectrum scales as s −2 .
We can now perform the same exercise in velocity space, viz., where δf u is the characteristic amplitude of δf at velocity scale u.On the other hand, using (43) and assuming that the integral over k of F is dominated by the region ξ ≲ 1, or we have Requiring consistency between ( 53) and ( 51) yields b = 6/(α + 1).
Assembling the scalings that we have surmised above, we find that the spectrum (43) is For α = 2, this was derived, by means of a formal solution, in [35], but the phenomenological argument and physical interpretation presented here are new.Note that the dimensional factors in (54) come from using Parseval's theorem (31) together with demanding that the integrals of the spectrum in ( 50) and ( 53) satisfy the constant-flux relations in position and velocity space, respectively, viz., (48) and (51).As an example, we plot the spectrum (54) for α = 2 in Fig. 1.
In summary, we have constructed a phenomenological theory according to which δC 2 undergoes a phase-space cascade in both position and velocity space.At large k (ξ ≫ 1), the spectrum is phase-mixing-dominated and has a power-law scaling in k.At large s (ξ ≪ 1), the spectrum is turbulent-diffusion-dominated and has a powerlaw scaling in s.These two regimes are separated in phase space by the critical-balance region ξ ∼ 1, where the linear and nonlinear time scales are comparable.The 1-D k and s spectra are dominated by contributions from this critical-balance region.
Since the collision operator is diffusive in velocity space and ν is assumed to be small, collisional dissipation must necessarily occur at fine velocity-space scales (large s).It is therefore unsurprising that a kinetic system with injection of a quadratic invariant exhibits a constant-flux cascade of that invariant in velocity space.However, there is no a priori scale in position space where the dissipation must happen, so it is nontrivial that there also exists a cascade in position space.This inertial range in position space emerges due to the nonlinear field-particle interactions in (33), which mix position and velocity space [35].
A reader interested in how the above results are obtained more rigorously should read Sections IV and V, where we solve (33) properly for the spectrum and fluxes of δC 2 .A reader interested only in the big picture can skip straight to Section VI.

IV. PHASE-SPACE FLUXES
We now verify the phenomenological theory presented in Section III by solving (33).It is useful to write this equation in flux-gradient form: where ∇ = (∂/∂k, ∂/∂s) is a gradient operator in the (k, s) space, and the flux Γ = Γk , Γs has components which is clear from (33), and The latter expression can be derived using (36), viz., by writing the Fourier transform of the nonlinear term as Inverse-Fourier transforming the term in brackets in ( 58) yields (57) [53].In steady state, (55) reads In Sections IV A and IV B, we compute Γk and Γs integrated over s and k, respectively.Then, in Section V, we compute the spectrum F (k, s) and analyze the full flux Γ in (k, s) space.This order of presentation may seem awkward, but is, in fact, necessary.This is because the integrated fluxes inform us of the nature of the solution of ( 59).This solution is in turn needed to compute the 2-D flux in (k, s) space.

A. Constant flux in k
Let us integrate (59) over all s.The s flux vanishes at s → ±∞, and we are left with an equation for ĝ where ).This integrated source Ŝ injects ĝ, which is diffused by the nonlinear term and dissipated by collisions.
The source ( 34) is peaked at low k, with characteristic width L −1 E .To analyze the behavior of ĝ(k) in the region kL E ≫ 1, we approximate Ŝ(k) ≈ (ε/L) δ(k).Fourier transforming (60) and solving for g(r), the Fourier transform of ĝ(k), yields where we have defined the 'Kolmogorov' (dissipation) scale as For α = 2, (62) simplifies to For α < 2, (62) does not have a simple closed-form expression, but it can be manipulated into an integral where the exponential in the integrand is decaying rather than oscillatory: This will be useful in what follows.Deriving ( 65) from (62) requires some work, which is done in Appendix C. The solution (62) implies that the rate of δC 2 injection by the source equals the rate of δC 2 dissipation by collisions.Indeed, using (32) and (62), and converting the sums over k into integrals, the collisional dissipation can be written in terms of ĝ: as expected.Since ( 62) is a Green's function solution to (60), it is straightforward to show that this balance also holds for arbitrary Ŝ(k).Importantly, (66) applies even in the limit ν → 0 + ; emergence of such finite collisional dissipation in the collisionless limit is known as a dissipative anomaly [29].This result, although perhaps obvious from the steady state of (32), demonstrates constructively that the steady-state solution to ( 55) is well defined in the collisionless limit.
Using the solution (61), we can now compute the 1-D k flux of δC 2 , viz., from (57), For α = 2, this flux reduces to where we used the solution (64).For α < 2, The derivation of (69) can be found in Appendix C. When kℓ ν ≪ 1, both ( 68) and ( 69) are constant and equal to sgn(k) ε/2 (in this limit, the integral over z in ( 69) is π/2, which is also shown in Appendix C).This result, one of the most important of this paper, means there is a constant-flux cascade in position space, viz., there exists an inertial range, 1/L E ≪ k ≪ 1/ℓ ν , unaffected directly by forcing or collisions, where δC 2 is transferred from larger to smaller scales.At kℓ ν ≳ 1, collisions become relevant.The dissipation rate at scales ≳ 1/k is Substituting (64) or ( 65) into (70) gives, for α = 2, or, for α < 2, (72) From ( 71) and ( 72), we see that the collisional dissipation is only order-unity when kℓ ν ≳ 1, i.e., below the Kolmogorov scale (63), which, therefore, deserves the name that we have given it.Past kℓ ν ≳ 1, the dissipation balances the δC 2 injected by the forcing ( D(k) ≃ ε when kℓ ν ≫ 1; derived in Appendix C).As discussed at the end of Section III, because the collision operator is diffusive in velocity space, there is no a priori scale in position space where the dissipation must happen.Rather, this dissipation range in position space forms because of the collisionless dynamics.Note that ℓ ν → 0 as ν → 0, so arbitrarily fine-scale structure in position space can be generated in the collisionless limit.Because of the constantflux cascade, all of the δC 2 injected at large scales reaches the dissipation range, no matter how small ℓ ν is.

B. Constant flux in s
Let us now integrate ( 59) over all k.The k flux vanishes at k → ±∞, resulting in the following equation for the s flux: Unfortunately, unlike (60), this equation is not closed and cannot be explicitly solved without knowing the spectrum itself.However, we can still learn a key lesson from it.In view of ( 34) and ( 35), the characteristic width over which dk Ŝ falls off in s is 1/v th .But for small ν, collisions can only be relevant when s is large.Balancing the phase-mixing term with the collision term in (59) tells us that collisions will start to matter when We know from Section IV A that collisional dissipation will start to occur around kℓ ν ∼ 1, so taking k ∼ 1/ℓ ν in (74) gives us the collisional velocity scale When su ν ≪ 1, we expect that we can drop the collision term in (73).Integrating the remaining terms in the equation from −s to s, where 1/v th ≪ s ≪ 1/u ν , yields, using (38), Because the distribution function is real, the Fourier spectrum satisfies Together with the definition (56) of Γs , this implies that dk Γs (−s) = − dk Γs (s), and so (76) becomes for 1/v th ≪ s ≪ 1/u ν .There is, therefore, also an inertial range in s where there is a constant flux of δC 2 from small to large |s|.Because the collision operator can only be relevant at large |s| when ν is small, the existence of a forced steady state does indeed require such a constant-flux inertial range in s.
We note that this result holds even in the absence of nonlinearity.Linear phase mixing can lead to the development of arbitrarily fine scales in velocity space, albeit at the price of a diverging δC 2 as ν → 0 + [51].However, nonlinearity is required for there to be a cascade in position space: note from (63) that ℓ ν → ∞ as κ → 0. Indeed, the solution of (60) for κ = 0 is simply ĝ = Ŝ(k)/2ν.The Vlasov equation is linear in this limit, so the spatial Fourier components of the distribution function are uncoupled from one another, and the k dependence of the Fourier spectrum is then simply set by the source.

V. PHASE-SPACE SPECTRUM
We now compute the spectrum in (59).From Sections IV A and IV B, we know that our model exhibits a phasespace entropy cascade in which the rate of δC 2 injection by the forcing at large scales is balanced by collisional dissipation at small scales, with the distribution function developing arbitrarily small scales in phase space in the limit ν → 0 + .This results in an inertial range in k (s) where the flux of δC 2 integrated over velocity (position) wavenumbers is constant.Therefore, in this inertial range, viz., for 1/L E ≪ k ≪ 1/ℓ ν and 1/v th ≪ s ≪ 1/u ν , we can meaningfully consider (59) in the absence of the forcing and collisional terms, viz., and seek a solution for the spectrum that supports constant fluxes in k and s.

A. Self-similar inertial-range solution
Because of the reality condition (77), we only need to solve (79) for F (k, s) on half of the (k, s) plane.We choose to solve for the spectrum in the upper half-plane, −∞ < k < ∞ and s ≥ 0 [72].
To deal with the fractional Laplacian, we Fourier transform this equation in k space.Using (36), we find that F (r, s) satisfies Because F (r, s) is the Fourier transform of F (k, s), which is purely real, F (r, s) must satisfy the reality conditions F (−r, s) = F * (r, s).We therefore only have to solve (80) for r > 0, for which |r| = r.
The inertial-range solution to (80) that does not depend on details of the forcing or dissipation ranges is a similarity solution [73]: (81) where β can be constrained by the fact that the flux L dk Γs must be constant in the inertial range, as per (78): The integration variable dual to y, has previously appeared in (42) as the ratio of the turbulent-diffusion and phase-mixing time scales.Substituting (81) into (80) gives us an ordinary differential equation for ϕ: To solve this equation, consider the transformation Then, (84) becomes a modified Bessel equation [74]: Therefore, the solution to (84) that vanishes at y → ∞ is where y > 0, K is the modified Bessel function of the second kind, and Λ is a constant.The reality condition ϕ(y) = ϕ * (−y) constrains us to pick the phase of Λ so that the real part of ϕ is even in y and the imaginary part is odd.This is accomplished by setting Λ = σ(e −iπ/4 ) 1/(α+1) , where σ > 0 is real and determined by the constraint (82).For generic α, we are unable to evaluate the integrals in (82) analytically, to compute σ, although it is straightforward to see that σ is finite.For α = 1, σ is easily computed: this calculation can be found in Appendix D.
Using the solution (87) for ϕ, we inverse-Fourier transform back to k space to obtain, where ξ is given by ( 83), and we have used the reality condition for ϕ.As we stated at the beginning of this calculation, this expression is valid only for s ≥ 0; the spectrum for negative s can be found by combining (77) and (89).For generic α, we are unable to simplify (89) further.For α = 1, (89) has a simple closed form, derived in Appendix D. For α = 2, the spectrum can be written in terms of incomplete Gamma functions: see [35].In our terms, (87) for α = 2 reduces to where Ai is the Airy function.This is clear by the fact that for α = 2, ( 84) is an Airy equation.

B. 2-D spectrum
To show what (89) looks like, we present contour plots of the normalized spectrum L 2 0 v 6/(α+1) th F (kL 0 , sv th ) for α = 1 (for which an explicit expression can be found in Appendix D) in Fig. 2. We use normalized units sv th and kL 0 , where This length scale is a natural choice, because the similarity variable (83) in these units is and, therefore, the spectrum, up to its amplitude, is independent of κ.
For ξ ∼ 1, the spectrum is a nontrivial function of k and s.For ξ small or large, it has asymptotics given by (54) (computed below), confirming the phenomenological theory presented in Section III.As discussed in Section III, the distinction between ξ small and large F (kL0, sv th ) vs. (kL0, sv th ), for α = 1.Note that the spectrum for s < 0 is given by the reality condition (77), F (k, s) = F (−k, −s).The black lines are the critical-balance lines (47), sv th = ±|kL0| 2/3 .In order to use logarithmic axes, we do not plot the spectrum along the k = 0 and s = 0 axes.We have also not plotted the spectrum near the (k, s) origin, as the similarity solution (89) diverges there.can be understood in terms of competition between the phase mixing and turbulent diffusion for control of the phase-space cascade.
To compute the asymptotics (54) from the full solution (89), the ξ ≪ 1 limit can be found by simply setting ξ = 0 in (89).For the ξ ≫ 1 limit, we use the following result from Fourier analysis.Suppose u(y) is a function that is smooth for y > 0 and has the scaling u(y) ∼ y λ−1 as y → 0 + , where λ > 0. Then as ξ → ∞ [75].Note that (87) has the series as y → 0 + , where A and B are real and positive constants.Combining ( 93) and ( 94), one finds that the contribution to the integral in (89) from the constant part of ( 94) is purely imaginary and so vanishes.The nextorder term from the linear part of (94) has a real component; the integral is therefore Note that, the 1-D k and s spectra, which can be found by integrating out the similarity variable ξ in (89), agree with (48) and (51).
Before continuing, it is instructive to assess the region of validity of (89) in the (k, s) plane.This spectrum is infrared-divergent and thus breaks down as (k, s) → (0, 0); information about the functional form of the source (34), which would regularize the spectrum at long wavelengths, is lost by construction in the similarity solution (however, (89) does contain information about the flux from the source via the constraint ( 82)).Of course, the inertial-range spectrum is no longer valid in the region where the source (34) is concentrated, viz., when sv th ≲ 1 and kL E ≲ 1.For these reasons, we have not plotted the spectrum near the origin in the (k, s) plane in Fig. 2, and likewise for the fluxes plotted in the following sections [76].
In addition to the solution (89) lacking an outer scale, the approximation that the nonlinear term is a fractional Laplacian to lowest order also breaks down when kL E ≲ 1.This is clear in the derivation of the nonlinear term in Appendix B 4, where the fractional Laplacian emerges as the lowest-order term when kL E ≫ 1.Yet, away from s = 0, ( 89) is, in fact, continuous across k = 0.By dropping the finite-kL E corrections in (33), we have shrunk the boundary layer kL E ≲ 1 to the point k = 0 [77].
We can also now address the concern of the lack of locality in k space of the source term (34) and whether dropping this term in the inertial range was justified, as discussed at the end of Section II C. Consider the two asymptotic regions, ξ ≪ 1 and ξ ≫ 1 (assuming also sv th ≫ 1, su ν ≪ 1, and kℓ ν ≪ 1).In the former region, the nonlinear term is dominant over the phase-mixing term, as per (44).Balancing the nonlinear term with (34) yields an inhomogeneous contribution to F ∝ ⟨ f ⟩ 2 .Since sv th ≫ 1, this term is strongly suppressed, e.g., exponentially so in the case of a Maxwellian initial condition [see (35)], so this term is negligible compared to the ξ ≪ 1 asymptotic in (54).In the the latter region (ξ ≫ 1), the phase-mixing term is dominant over the nonlinear term, as per (45).Balancing the phase-mixing term with (34) yields an inhomogeneous contribution to F ∝ D(k)/k, which is ∝ k −(2+α) when kL E ≫ 1.This contribution is subdominant to the homogeneous part of the spectrum, which scales like k −2 when ξ ≫ 1, viz., (54).We cannot explicitly estimate the inhomogeneous contribution to the spectrum from the source in the ξ ∼ 1 region, but the above analysis suggests it should be subdominant to (89).Therefore, even though the source is multiscale in k, there is still an inertial range in the position space.

C. 2-D flux: phase-space circulations
To gain insight about the pathways in phase space taken by F from injection to dissipation scales, it is informative to examine the vector flux Γ, which, in terms of the similarity solution (89), has the components (57), dy e −iξy y α−1 ϕ(y), (95) and (56),  To obtain (95), we used (81) and changed variables from r to y in the integral.Note that these expressions are valid only for s ≥ 0. For s < 0, combining ( 56), (57), and ( 77), we have that As was the case for the Fourier spectrum, the flux is a nontrivial function of k and s for ξ ∼ 1, but can be simplified when ξ is small or large.The asymptotics of the k and s components of the flux, which we derive below, are and respectively.Here, we have retained the signs of terms as well as dimensional factors in ( 98) and ( 99), but not order-unity constants.To evaluate the asymptotics of the fluxes for s < 0, these expressions must be combined with (97).We can derive these results in the same way as we did the asymptotics of the spectrum in Section V B. The asymptotics (99) for Γs come directly from combining ( 56) and (54).For Γk , the ξ ≪ 1 expansion in (98) comes from evaluating (95) at ξ = 0 (note that the integral is positive).For ξ ≫ 1, we can combine ( 93), (94), and (95).This gives, to lowest order, as ξ → ∞, The lowest-order k flux vanishes when α = 2, so (100) only gives the ξ ≫ 1 limit in (98) for α < 2. We need to go to next order for α = 2, which yields giving the ξ ≫ 1, α = 2 asymptotic in (98).We first discuss the Batchelor case (α = 2), where To Thus, when ξ ≪ 1, the flux is diffusion-dominated (dominated by its k component), and when ξ ≫ 1, the flux is phase-mixing-dominated (dominated by its s component).When ξ ∼ 1, the two fluxes are comparable; this can be seen in the plot of R in Fig. 3.The regions of the dominance of the two fluxes are, therefore, separated in phase space by the critical-balance line (47), the same as the phase-mixing and turbulent-diffusion time scales, viz., (42).We plot the (nondimensionalized) vector flux (102) in Fig. 4(a).Using the asymptotics (104) and noting the signs of the fluxes in ( 98) and ( 99), we see that in the diffusion-dominated region (ξ ≪ 1), Γk is negative when s > 0 and positive when s < 0, and in the phase-mixing-dominated region (ξ ≫ 1), Γs is positive when k > 0 and negative when k < 0. The flux therefore gives rise to counterclockwise circulation of δC 2 in (k, s) space.The sign changes in the components of the flux that enable this circulation occur at ξ = ξ 2 ≃ 0.747, where R has a zero (as can be seen in Fig. 3), and at ξ = 0, below (above) which R diverges positively (negatively).The first point corresponds to Γk changing sign in the top right and bottom left quadrants, while the second point corresponds to Γs changing sign from positive to negative between the top right and top left quadrants, as well as between the bottom left and bottom right quadrants.This latter effect occurs because, when sgn(ks) = −1, perturbations are phase unmixed rather than phase mixed, being advected to low |s| rather than high |s| [35,61].The phase-unmixing modes are a stochastic instantiation of the plasma-echo effect [62,63].
While the phase unmixing does undo the phase mixing, we will see in Section V E that the flux of the phasemixing modes outweighs that of the phase-unmixing ones, thus enabling the constant-flux cascade.
In Fig. 4(a), this manifests in the fact that the circulation swirls outward.Indeed, in the top right (bottom left) quadrant, below (above) the line ξ = ξ 2 where Γk = 0, there is a flux of δC 2 to both high |k| and high |s| simultaneously, toward the dissipation wavenumbers k ν ∼ 1/ℓ ν and s ν ∼ 1/u ν .

D. Non-local transport
We now examine the α < 2 cases.The important difference compared to the Batchelor regime is that the k flux is now non-local in k space.Note that Γk can be written as [53] [78] where c α is given by (B19).The derivation of this expression is given in Appendix B 5. The interpretation of ( 105) is that 'particles' (parcels of δC 2 ) cross the point (k, s) from points (k±p, s), with cumulative probability ∝ p −α .The particles undergo Lévy flights in k space, so the flux at a point (k, s) receives contributions not just from nearby particles taking small jumps but also from faraway ones taking large jumps.The net effect is an enhancement of Γk compared to Γs .For these cases, the ratio R of the nondimensionalized fluxes is Since Γk is positive when ξ ≫ 1 and negative when ξ ≪ 1, viz., (98), there is always an order-unity ξ α at which Γk = 0. Unlike the Batchelor case (103), the ratio of fluxes ( 106) is not a function solely of ξ.Therefore, (106) implies that along curves of constant ξ ∼ 1 ̸ = ξ α , the flux is always diffusion-dominated as sv th → ∞.Furthermore, in the asymptotic region ξ ≫ 1 (which for α = 2 was the phase-mixing-dominated region), we show in Appendix E that, as α gets smaller, the region where the flux is asymptotically phase-mixing-dominated shrinks.In fact, for α < 1/2, there is no asymptotic region at all where the flux is phase-mixing-dominated (except along the curve ξ = ξ α ).
As an example, we plot the (nondimensionalized) vector flux for α = 1 in Fig. 4(b) (explicit expressions for the components of the flux in this case can be found in Appendix D).In this case, apart from along the curve ξ = ξ 1 where Γk = 0, the flux is phase-mixing-dominated only in the region |sv th | ≲ 1 (irrespective of k), which in Fig. 4(b) we indicate with gray lines.
These results do not mean that there is no effective phase mixing for these cases.The ratio (42) of nonlinear and linear time scales gives the relative local transport of δC 2 in s versus k space.While the non-local transport of δC 2 in k space dominates over its flux in s, locally, the phase-mixing time (40) of F is still shorter than the diffusion time (41) when ξ ≫ 1.This dominant local phase mixing is what sets up the lowest-order, constant-flux-ins spectrum in the ξ ≫ 1 region, viz., (45) and (54).
The fluxes also obey critical balance.It is straightfor-ward to show, using ( 98) and ( 99) with calculations analogous to (50) and (53), that the 1-D fluxes ( 69) and ( 78) are dominated by contributions from the critical-balance region (47) (ξ ∼ 1).Even though the 2-D s flux is subdominant to the k flux, the fact that L dk Γs is constant in the s inertial range implies that phase mixing still provides an effective route to dissipation scales in velocity space.

E. Shell-averaged flux
To understand the net effect of having both phasemixing modes that propagate from low to high |s| and phase-unmixing modes that propagate from high to low |s|, it is useful to consider the flux shell-averaged in k, Note that in 1D, shell averaging amounts simply to adding together contributions from +k and −k.The flux ( 107) is defined so that the shell-averaged spectrum While the components of (107) depend on α, their (nondimensionalized) ratio does not (up to a prefactor).Note that, using (84) and integrating by parts, (95) can be rewritten as where the derivative of ϕ is taken at y → 0 + .Using ( 96) and ( 109), we get Γk (kL 0 , sv th ) − Γk (−kL 0 , sv th ) Γs (kL 0 , sv th ) + Γs (−kL 0 , sv th ) Remarkably, this expression is valid everywhere in the (k, s) plane, independent of ξ being small or large.The flux is radial at both large sv th and large kL 0 .We plot Γ for α = 1 in Fig. 5, which clearly exhibits this feature.
Note that the components of the shell-averaged flux are positive-definite.We are only able to show directly that this property holds for α = 1; this calculation can be found in Appendix D. An argument as to why the fluxes are positive-definite for the general case is as follows.Since (110) is positive, the components of (107) are either both positive or both negative for all (k, s) in the inertial range (no sign reversals are possible, otherwise the flux would not be divergence-free in the inertial range: see ( 59)).If the components were negative-definite, there would be a sink at the origin.However, there is a source at the origin, so the shell-averaged fluxes must therefore be positive-definite.This positive definiteness is important.The circulatory nature of the fluxes in Fig. 4 is due to phasemixing modes (with sgn(ks) = 1) and phase-unmixing modes (with sgn(ks) = −1) propagating in opposite directions in s.Since the s component of ( 107) is equal to k[ F (k, s) − F (−k, s)], the shell-averaged flux being positive-definite means that the spectral amplitudes of the phase-mixing modes are greater than those of the phase-unmixing modes.Therefore, by adding together the fluxes of the two modes, we are left with a net flux that points outward to both high k and high s toward the dissipation scales, in agreement with our analysis in Sections V C and V D. Γs (kL0, sv th ).For emphasis, we specify the magnitude of the flux by both the length and color of the vectors.The black curve is the critical-balance curve (47) for these parameters, sv th = (kL0/ξ1) 2/3 , where ξ1 = 1/ √ 3.

VI. CONCLUSION A. Summary
In this paper, we have presented a solvable model of kinetic plasma turbulence, in which the electric field is decoupled from the particle distribution function and taken to be an externally imposed Gaussian field, white-noise in time and power-law in k space.
The effect of this stochastic electric field on the mean distribution function is diffusion in velocity space, often referred to as stochastic heating [6,[40][41][42].The resulting energization of particles is a collisionless process.Indeed, the heating rate is set by the turbulent collisionality and is independent of ν [see (15)].However, the irreversibility of stochastic heating hinges on the presence of collisions.As ⟨f ⟩ heats, δf fluctuations are excited, transferring the minus 'entropy' C 2,0 = (1/L) dxdv ⟨f ⟩ 2 /2 into δC 2 = (1/L) dxdv ⟨δf 2 ⟩/2, which then cascades to small scales in both position and velocity space simultaneously.This cascade is then cut off by collisions at fine phase-space scales, thereby rendering the heating irreversible.The irreversibility of stochastic heating is therefore enabled by the entropy cascade.
We have analyzed this cascade in the Fouriertransformed (k, s) space, where an 'inertial range' forms to bridge the injection and dissipation scales of δC 2 .Integrated over k or s, the flux of δC 2 is constant in this inertial range.Importantly, there is no collisional dissipation at scales much larger than the dissipation scale (ℓ ν , u ν ) [see (63) and ( 75)], which tends to (0, 0) as the collision frequency does.
In the 2-D (k, s) space, the Fourier spectrum of δC 2 has a self-similar profile, with power-law asymptotics at high k and s, respectively.We find that these asymptotic scalings can be deduced by a phenomenological theory whose governing principle is that the cascade satisfies a critical balance in phase space between the time scales of linear phase mixing and turbulent diffusion.Because there is nothing in our phenomenological theory that is unique to 1D-1V, we also expect these ideas to apply in 2D-2V and 3D-3V.While the one-dimensional |k| and |s| spectra (48) and ( 51) should be the same, the two-dimensional |k|-|s| spectrum (integrated over angles) will not be the same as (54) because of the wavenumber Jacobian being different from unity in higher dimensions.

B. Fast cascade and the effectiveness of phase mixing
In a linear system, phase mixing acts as a route to collisional dissipation at every spatial scale [50,51], but in the model presented here, collisional dissipation is only non-negligible below the dissipation scale ℓ ν .Following the commonly held intuition that the effect of phase mixing (and Landau damping) on turbulent systems is that it steepens the spectrum (of, e.g., electromagnetic energy) at every scale by an amount set by the Landau-damping rate [27,[79][80][81][82][83][84], one may be led to conclude from our results that phase mixing is less effective in a nonlinear system than in a linear one.However, in fact, phase mixing is even more effective in a nonlinear system.This is because the presence of nonlinearity produces fine structure in position space and thus enhances the rate at which the distribution function develops fine phase-space gradients, reducing the time that it takes collisions to activate, τ ν .
To see this, note that in a steadily driven system, as considered in this paper, the total δC 2 in steady state divided by the injection rate ε gives a reasonable estimate for τ ν .For a linear system, restricting ourselves to a single k and noting that the spectrum is flat in s [85] up to a collisional cutoff s ν ∝ ν −1/3 given by the balance (74) [50,51], we find [86] For smaller and smaller ν, it takes longer and longer for collisions to dissipate the velocity-space structure of the distribution function, and in the collisionless limit, the amount of δC 2 stored in phase space diverges [51].This is consistent with a constant cascade time, set by the phase-mixing time, viz., τ c ∼ (k v th ) −1 (with k fixed).
In the presence of nonlinearity, when α = 2, the 1-D spectra (48) and ( 51) scale like k −1 and s −1 , respectively [35], the former scaling in agreement with the classical Batchelor [65] spectrum of a passive scalar advected by a single-scale flow.Integrating the 1-D spectrum up to the collisional cutoff (63) gives Although formally also divergent as ν → 0 + , ( 112) is asymptotically shorter than (111).In this case, the cascade time is also constant, viz., τ c ∼ κ −1/3 , as can be seen in ( 46), even though many k's are involved.When α < 2, the phase-space cascade is even more efficient.From (46), we have that τ c goes to zero as k and s go to infinity, so fluctuations 'turn over' faster to finer phase-space scales the deeper they are in the inertial range (compared to the constant cascade times in the linear and Batchelor regimes).As a result, the 1-D spectra (48) and (51) are steeper than k −1 and s −1 , so the total steady-state δC 2 and τ ν are independent of ν.
Thus, the presence of nonlinearity reduces the collision time from a (negative) fractional power of ν (111) in the linear regime to a time scale that is only logarithmic with ν (112) when α = 2, or one that is independent of ν when α < 2. To interpret this shortening of τ ν in the nonlinear regime, note that linear phase mixing still processes δC 2 from injection to dissipation scales in the inertial range, but so does, simultaneously, nonlinear mode coupling, in a critically-balanced fashion.The net result is fast dissipation, but only at wavenumbers k ≳ k ν ∼ 1/ℓ ν and s ≳ s ν ∼ 1/u ν , which, by construction, satisfy the critical-balance condition (47).The nonlinear cascade ensures that all the injected δC 2 flux at large phase-space scales is rapidly dissipated at small scales via collisions, no matter how small is ν.This reduction of τ ν is also an a posteriori justification of the assumption (25) of there being a separation of time scales between the time that it takes δC 2 to reach quasi-steady state and the time that it takes for the injection rate ε to decay due to the stochastic heating of ⟨f ⟩ [87].

C. Implications and outlook
As discussed in Section VI B, these results provide a conceptual understanding of the role of phase mixing in turbulent plasmas.Recent theoretical [35,61] and numerical [88,89] studies have suggested a statistical 'fluidization' of turbulent, collisionless plasmas by stochastic plasma echoes suppressing phase mixing in the (spatial) inertial range.This might have seemed to be at odds with Landau damping [90] clearly being identified in turbulent settings numerically [84,[91][92][93][94][95] in several works, as well in Magnetospheric Multiscale (MMS) mission observations of the turbulent plasma in the Earth's magnetosheath [96,97].While our model is quite simplified and does not include self-consistent electric fields and hence Landau damping, insofar as Landau damping and phase mixing are intimately related processes, our results indicate that these two seemingly contradictory sets of results can in fact be compatible.
This work also has implications for the relaxation of mean distribution functions in nearly collisionless plasmas.The existence of the entropy cascade implies that collisions will always dissipate fine-scale structure in the distribution function, even when ν is vanishingly small, but it does not necessarily imply that the rate by which the distribution function relaxes toward a Maxwellian is enhanced.This is clear in our model from the fact that, whereas δf develops sharp phase-space gradients, ⟨f ⟩ does not, so Coulomb collisions are never activated for it.Mean distributions in space plasmas can be highly non-Maxwellian, e.g., in the solar wind [1,6,98].Developing a theoretical formalism to predict the form of such non-equilibrium distribution functions is an outstanding problem.A direct consequence of entropy cascades is that theories of relaxation that assume phase-volume conservation [99][100][101][102][103] may not apply to nearly collisionless plasmas that are strongly turbulent.An alternative approach is to examine how the turbulent phase-space correlations of δf drive the evolution of the mean distribution function [43,104].
There is much opportunity to understand phase-space entropy cascades in nearly collisionless plasmas better, with theory, numerical simulations, laboratory experiments, and spacecraft data.With regards to the latter two, the works [30,[105][106][107] suggest that measuring entropy cascades in real plasmas is a realizable endeavor just beginning to be possible.M.L.N is grateful to I. Abel, T. Adkins, M. Allen, W. Clarke, S. Cowley, P. Dellar, J.-B.Fouvry, M. Kunz, R. Meyrand, F. Rincon, and J. Squire for helpful discussions related to this work.W.S would like to thank A. Bhattacharjee and H. Weitzner for useful discussions.We are also grateful to two anonymous referees, whose feedback improved the manuscript.M.L.N was supported by a Clarendon Scholarship.R.J.E was supported by a UK EPSRC studentship.W.S. was supported by a grant from the Simons Foundation/SFARI (560651, AB).The work of A.A.S. was supported in part by UK EPSRC (grant EP/R034737/1), STFC (grant ST/W000903/1), and the Simons Foundation via a Simons Investigator award.W.D.D. and M.L.N were supported by the US Department of Energy through grant DEFG0293ER54197 and Scientific Discovery Through Advanced Computing (SciDAC) grant UTA18000275.In this appendix, we discuss invariants of the Vlasov equation alternative to C 2 .It is straightforward to show that (2) conserves an infinite number of invariants, socalled Casimir invariants [21].Indeed, focusing on invariants defined in the averaged sense, for any smooth g(f ), the functional satisfies, using (7) for the collision operator, where the term with two derivatives on δf dominates when ν → 0 + .Note that while in this paper ⟨...⟩ denotes ensemble averaging with respect to the stochastic electric field, the arguments in this section hold for any averaging procedure for which one can decompose f = ⟨f ⟩ + δf , where ⟨δf ⟩ = 0.When ν = 0, every G[f ] is formally conserved.Furthermore, if g(f ) is convex, i.e., g ′′ (f ) ≥ 0 everywhere, then the corresponding Casimirs are negative-definitely dissipated by collisions.We label Casimirs with positivedefinite time evolution (which are minus the convex functionals of f ) as 'generalized entropies' (cf.entropy functions in hyperbolic partial differential equations [108]).This set includes −C 2 , as well as the traditional entropy S = − dxdv f log f .We define the relative entropy as which has the budget which is negative-definite if g(f ) is convex.Therefore, we can interpret (A4) analogously to (24) We have shown that C 2 is anomalously dissipated as ν → 0 + and is cascaded, i.e., exhibits an inertial range unaffected directly by collisions or forcing.In principle, invariants other than C 2 can also have these properties.A system mathematically similar to (2) in which this happens is the advection-diffusion equation for a scalar advected by a turbulent flow: [109] showed that the family of invariants dx θ 2n , where θ is a scalar field and n is a positive integer, satisfy constant-flux cascades; in contrast, for a passive scalar advected by a smooth, chaotic flow (Batchelor regime), only the quadratic invariant (n = 1) is cascaded.This is because for higherorder invariants, logarithmic correlations of the passive scalar give rise to injection of those invariants by the source at all scales, preventing the formation of an inertial range.
It is likely that a similar situation happens in the Vlasov-Kraichan model between the cases α < 2 and α = 2, but such a calculation is beyond the scope of this work.If it were true, then for α < 2, the cascade of C 2 does not necessarily hold deeper physical meaning than the cascade of any other convex functional of f .However, we still believe that C 2 is a particularly useful quantity in kinetic plasma turbulence.Because it is quadratic in f , C 2 is the only invariant (up to weight functions in the integrand of (A1)) that satisfies and so the relative entropy (A3) is a function solely of δf .This property is useful for conceptualizing the budget (A4) as a transfer of entropy between ⟨f ⟩ and δf , as any other Casimir invariant involving higher powers of f will necessarily involve cross terms containing both ⟨f ⟩ and δf .Furthermore, C 2 is the only invariant that lends itself to a simple Fourier analysis.For these reasons, in this paper, we have chosen to analyze phase-space turbulence using C 2 exclusively.
Appendix B: Derivation of the Fourier-spectrum equation In this appendix, we give a detailed derivation of (33).Multiplying (29) by δ f * (k, s), adding to the resulting equation its complex conjugate, and then ensemble averaging gives Note that the second term in the Fourier sum in (29) vanishes under multiplication by δ f * (k, s) and ensemble averaging.The source term Ŝ, defined in (34), comes from the second term on the right-hand side of (29).It can be found via application of the Furutsu-Novikov theorem (10) and using the fact that (3) implies where δ k,−k ′ is the Kronecker delta.
1. Derivation of (32) The δC 2 budget (32) in terms of F can be found by taking the time derivative of (31) and using (B1).Assuming the spectrum goes to zero at s → ±∞, the freestreaming term vanishes by integration over s.For the nonlinear term, note that the summand in Im k,p Ê(p)δ f (k − p, s)δ f * (k, s) , (B3) after taking p → −p and k → k − p, and applying reality conditions on the electric field, is equal to its own complex conjugate.Therefore, its imaginary part vanishes.What is left is injection by the source and dissipation by collisions, viz., (32).where (...) represents terms that will vanish after we take the functional derivative.Therefore, we get where H is the Heaviside step function, defined with the convention that H(0) = 1/2.Combining this expression with (B2) and (B4) yields Let us simplify the nonlinear term in (B8) further.We start with the case of α = 2.In the limit kL E ≫ 1, we suppose (and check a posteriori ) that ( 5) is sufficiently steep in wavenumbers that we can consider k ≫ p and Taylor-expand the summand of the wavenumber sum in (B8): This 'Batchelor approximation' was first used for the problem of passive-scalar mixing in fluids [65,110] and amounts to approximating the electric field as effectively single-scale.Substituting (B9) back into the sum in (B8), the first term vanishes because it is odd in p, and we are left with (B11) In the limit η → 0 + , (B11) is logarithmically divergent, being ∝ log (L E /η); without a small-scale cutoff, the approximation (B10) is invalid.This is because the k −3 spectrum in ( 5) is only marginally in the Batchelor regime [35,65].The Batchelor limit generically applies when the electric field is spatially smooth, corresponding to a rapidly decaying spectrum D(k).We choose the particular form of D(k) in (5) in order to match onto the Batchelor limit and fractional cases with one functional form of the correlation function, but we could just as well have picked a steeper (5) for α = 2 that would not have required a small-scale cutoff.Therefore, without loss of generality, we keep (B11) without modification.
4. The α < 2 cases: representation in terms of the fractional Laplacian For α < 2, it is convenient to manipulate the nonlinear term in (B8) in position space and then Fourier transform back to k space.We start by noting that (B12) We now take the limits η → 0 + and kL E ≫ 1.When η = 0, note that (4) is the kernel of the Bessel potential [64,111], and can therefore be written as where K is the modified Bessel function of the second kind, and Γ is the Gamma function.For α < 2 and r/L E ≪ 1, (B13) has a series expansion where We now Fourier transform back to k space.To lowest order in (r/L E ) 2−α ≪ 1 (equivalently, (kL E ) 2−α ≫ 1), the r α term is dominant over the r 2 term, which gives and p.v. means that the integral is defined in the principal-value sense.Thus, using (B17) and (B10), we have that (B8) becomes (33), where with D α given by (B16) and D 2 given by (B11).

Derivation of (105)
Here, we derive the form (105) of the non-local k flux in the α < 2 cases, as analyzed in Section V D.
We first rewrite (B18) as an integral over positive p: contour along the real line.However, this residue contribution vanishes when δ → 0, so whether the contour is deformed below or above the pole does not change the final answer.
Given the expression (69), we can take its asymptotics for kℓ ν ≪ 1, as analyzed in Section IV A. In this limit, the inertial-range flux is approximately equal to the rate of δC 2 injection.Changing variables to u = z α , we have

L 4 0 0
work out by what physical mechanism the phase-space cascade is enabled in different parts of the (k, s) plane, it is useful to compare the ratio of the k and s fluxes.Using the normalizations from Section V B, we have that the ratio of the dimensionless fluxes is R = Γk (kL 0 , sv th ) dy e −iξy yϕ(y) Re ∞ 0 dy e −iξy ϕ(y) , (103) which is a function solely of ξ.It has the asymptotics
. In the absence of collisions, as G[⟨f ⟩] decreases via stochastic heating, R[f ] increases to maintain the G[f ] balance.Once δf has developed sharp enough gradients, collisions dissipate the total G[f ].