The relaxation of a weakly collisional plasma which is of fundamental interest to laboratory astrophysical plasmas can be described by the selfconsistent BoltzmannPoisson equations with the LenardBernstein collision operator. We perform a perturbative linear and secondorder analysis of the BoltzmannPoisson equations and obtain exact analytic solutions which resolve some longstanding controversies regarding the impact of weak collisions on the continuous spectra the discrete Landau eigenmodes and the decay of plasma echoes. We retain both damping and diffusion terms in the collision operator throughout our treatment. We find that the linear response is a temporal convolution of two types of contribution a continuum that depends on the continuous velocities of particles crucial for the plasma echo and another consisting of discrete modes that are coherent modes of oscillation of the entire system. The discrete modes are exponentially damped over time due to collective effects or waveparticle interactions Landau damping as well as collisional dissipation. The continuum is also damped by collisions but somewhat differently than the discrete modes. Up to a collision time which is the inverse of the collision frequency inlineformulammlmathmmlmsubmmlmimmlmimmlmi mathvariantnormalcmmlmimmlmsubmmlmathinlineformula the continuum decay is driven by the diffusion of particle velocities and is cubic exponential occurring over a timescale inlineformulammlmathmmlmrowmmlmommlmommlmsubsupmmlmimmlmimmlmi mathvariantnormalcmmlmimmlmrowmmlmommlmommlmn1mmlmnmmlmommlmommlmn3mmlmnmmlmrowmmlmsubsupmmlmrowmmlmathinlineformula. After a collision time however the continuum decay is driven by the collisional damping of particle velocities and diffusion of their positions and occurs exponentially over a timescale inlineformulammlmathmmlmrowmmlmommlmommlmsubmmlmimmlmimmlmi mathvariantnormalcmmlmimmlmsubmmlmrowmmlmathinlineformula. This slow exponential decay causes perturbations to damp the most on scales comparable to the mean free path but very slowly on larger scales. This establishes the local thermal equilibrium which is the essence of the fluid limit. The longterm decay of the linear response is driven by the discrete modes on scales smaller than the mean free path but on larger scales is governed by a combination of the slowly decaying continuum and the least damped discrete mode. This slow exponential decay implies that the echo which results from the interference of the continuum response to two subsequent pulses is detectable even on scales comparable to the mean free path as long as the second pulse is introduced within a few phasemixing timescales after the first.
['Plasma Physics', 'Physics - Plasma Physics', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Solar and Stellar Astrophysics', 'High Energy Physics - Experiment']
Relaxation of weakly collisional plasma Continuous spectra discrete eigenmodes and the decay of echoes
{'Uddipan Banik ∗': 'Department of Astrophysical Sciences, Princeton University, 112 Nassau Street, Princeton, NJ 08544, USA School of Natural Sciences, Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA Perimeter Institute for Theoretical Physics, 31 Caroline Street N., Waterloo, Ontario, N2L 2Y5, Canada', '†': 'Department of Astrophysical Sciences, Princeton University,  \nAmitava Bhattacharjee 112 Nassau Street, Princeton, NJ 08540, USA', '(Dated: November 13, 2024)': 'The relaxation of a weakly collisional plasma, which is of fundamental interest to laboratory and astrophysical plasmas, can be described by the self-consistent Boltzmann-Poisson equations with the Lenard-Bernstein collision operator. We perform a perturbative (linear and second order) analysis of the Boltzmann-Poisson equations, and obtain exact analytic solutions which resolve some long-standing controversies regarding the impact of weak collisions on the continuous spectra, the discrete Landau eigenmodes, and the decay of plasma echoes. We retain both damping and di ff usion terms in the collision operator throughout our treatment. We find that the linear response is a temporal convolution of two types of contribution: a continuum that depends on the continuous velocities of particles (crucial for the plasma echo), and another, consisting of discrete modes that are coherent modes of oscillation of the entire system. The discrete modes are exponentially damped over time due to collective e ff ects or wave-particle interactions (Landau damping), as well as collisional dissipation. The continuum is also damped by collisions, but somewhat di ff erently than the discrete modes. Up to a collision time, which is the inverse of the collision frequency ν c, the continuum decay is driven by the di ff usion of particle velocities and is cubic exponential, occurring over a timescale ∼ ν -1 / 3 c . After a collision time, however, the continuum decay is driven by the collisional damping of particle velocities and di ff usion of their positions, and occurs exponentially over a timescale ∼ ν c. This slow exponential decay causes perturbations to damp the most on scales comparable to the mean free path, but very slowly on larger scales. This establishes the local thermal equilibrium (LTE), which is the essence of the fluid limit. The long-term decay of the linear response is driven by the discrete modes on scales smaller than the mean free path, but, on larger scales, is governed by a combination of the slowly decaying continuum and the least damped discrete mode. This slow exponential decay implies that the echo, which results from the interference of the continuum response to two subsequent pulses, is detectable even on scales comparable to the mean free path, as long as the second pulse is introduced within a few phase-mixing timescales after the first.', 'I. INTRODUCTION': 'The relaxation of plasma is a complicated many-body problem with rich phenomenology. Many astrophysical and space plasmas are nearly collisionless and are often described by the collisionless Boltzmann (also known as the Vlasov) equation, widely discussed in textbooks. However, no matter how infrequent the collisions, no plasma in nature or a terrestrial laboratory is strictly collisionless, so in practical terms one may also treat a collisionless plasma as the zero-collision limit of a weakly collisional plasma. Furthermore, if the collision operator is of the Fokker-Planck type, the collision frequency typically multiplies the highest velocity derivative in the kinetic equation. Taking the limit of zero collisions is then a subtle problem, requiring the use of singular perturbation theory. It should then not be surprising if the results obtained in the zero-collision limit of the weakly collisional problem di ff er qualitatively from those obtained by solving the collisionless problem. The physical e ff ect of weak collisions is to regulate  \nextreme filamentation in phase space, much in the same way as viscosity tends to regulate near-singular vortex structures in fluids governed by the Navier-Stokes equations.  \nThese considerations have deep implications for the thermodynamics of many-body systems. It is well known that the Vlasov equation does not possess an H-theorem, but the Boltzmann equation does. The Vlasov equation has a denumerably infinite number of Casimirs (invariants). How many of them (if any) survive in the weakly collisional limit? And how do they influence the turbulent relaxation of a weakly collisional plasma? The Vlasov equation, despite the occurrence of Landau damping, is time-reversible, as demonstrated by the phenomenon of plasma echoes. How do weak collisions a ff ect the time-evolution of echoes?  \nThe questions posed above do not just belong to the esoteric realm of mathematical theory but are germane to a variety of weakly collisional astrophysical and space plasma phenomena, which can be described using kinetic theory. For example, the problem of decay of plasma echoes in the presence of weak collisions, which is a fundamental problem in plasma physics, has important implications for the role echoes can play in competing with phase-mixing and Landau damping in plasma turbulence [1]. This, in turn, has significant implica-  \ntions for understanding the energetics of plasma turbulence and interpreting observations and numerical experiments pertaining to, for example, solar wind heating and acceleration [2-4], reconnection in a nearly collisionless magnetosphere [5], formation of coherent structures such as electron holes in phase space [6], or particle heating near a black hole [7]. Of course, most of these scenarios involve large perturbations and non-linear e ff ects, for which numerical simulations and theories of turbulent cascade [8] are the way forward. However, before addressing the fully nonlinear problem, we must have a holistic understanding of the perturbative relaxation of weakly collisional plasmas, especially the spatio-temporal evolution of the distribution function, coupled self-consistently with the Poisson equation (if the plasma dynamics is assumed to be electrostatic). This paper addresses some of the open issues on this topic.  \nUnderstanding the modal response of a weakly collisional plasma is critical to the development of a statistical description of the system (be it in equilibrium or non-equilibrium). Recently, there has been a resurgence of interest in the longtime behavior of weakly collisional kinetic plasmas. A rigorous linear theory of weakly collisional plasmas [9-11], which builds on the seminal work of [12], demonstrated that a plasma governed by the Lenard-Bernstein (LB) collision operator (a special case of the Fokker-Planck operator) oscillates in the singular limit of zero collisions as a complete set of discrete normal modes at long times. A striking finding of these studies was that the Case-van Kampen modes, which are known to be a complete set of normal modes for a purely collisionless plasma and constitute a continuous spectrum, are eliminated when weak collisions are introduced into the system. The new eigenmodes of a weakly collisional plasma constitute a discrete spectrum, and include the Landau solutions. As is well known, the Landau solutions are not eigenmodes of the collisionless Vlasov-Poisson system and are instead remnants of the initial perturbations in the time-asymptotic limit. However, in the presence of weak collisions, these become eigenmodes (referred to hereafter as Landau modes 1 ). For a stable plasma, these eigenmodes undergo a simple exponential decay with time. On the other hand, [13] showed in a classic paper that the linear response to perturbations decays cubic exponentially due to collisions (at least at short time). Although they do not make it explicit in their paper, their analysis only applies to the continuum part of the plasma response, as we clarify in section IV. Their findings imply that the plasma echo, which arises non-linearly from the continuum response, would rapidly decay due to electron-ion collisions. The [13] prediction for the temporal decay of the continuum response is, however, much more rapid than the exponential decay of the discrete eigenmodes, identified by [9, 10]. [14] revisited the [13] result from a boundary-layer analysis of the problem, rea ffi rming the results of [13].  \nIt thus appears that the relaxation of weakly collisional plasmas has produced somewhat conflicting results: the superexponential (cubic exponential) decay of the continuum response and the plasma echo according to [13] and the exponential decay of the discrete eigenmodes according to [912, 15]. This naturally raises the following questions. Which of these results is correct? And, if both are correct, are these results simply di ff erent limits of the same response? This paper addresses these very questions. We show that the exact solution to the problem of complete perturbative (linear and second order) relaxation of a weakly collisional plasma indeed yields the above results as di ff erent asymptotic limits in time. We show that including the damping (also known as dynamical friction / drag [16]) term in the collision operator modifies the [13] prediction of cubic exponential decay of the continuum response to an exponential decay at late times and large spatial scales. We find that the [13] prediction only holds for short times, i.e., t < 1 /ν c, and scales smaller than the mean free path, i.e., λ < λ c = σ/ν c or ν c < k σ , where ν c is the collision frequency, σ is the equilibrium velocity dispersion, λ c = σ/ν c is the mean free path, k is the wavenumber of the response and λ = 2 π/ k is the wavelength. At long times ( t > 1 /ν c) and on large scales ( λ > λ c or ν c > k σ ), the response damps exponentially and much slower, at a rate similar to that predicted by [9, 10]. The damping is most pronounced for ν c ∼ k σ or λ ∼ λ c. The slow exponential decay outside the mean free path implies that the distribution function (DF) approaches the equilibrium Maxwellian form rapidly only within the mean free path; this is consistent with the assumption of the establishment of a local thermodynamic equilibrium (LTE) in the fluid approximation (large ν c). Hence, we recover the fluid limit from a rigorous theory of weakly collisional plasmas.  \nWe remark that the presence of two timescales in a problem of kinetic relaxation is, by itself, not a novelty. For example, as pointed out to us by an anonymous referee, the problem has been extensively investigated in the context of inertial confinement fusion (ICF) experiments which require drastic revisions of classical transport theory, and also provide two temporal behaviors, one which is super-exponential and the other, exponential (see [17-22] and other references therein). However, the time-dependence as well as the parametric dependence of the kinetic (collisionless) super-exponential decay of the linear response on plasma parameters obtained in these works are quite di ff erent from our findings on the weakly collisional regime of the linear response and the non-linear plasma echoes. For example, the kinetic super-exponential decay of the density and temperature perturbations found by [20, 22] is Gaussian, whereas we find a more general form of the decay in the weakly collisional regime that can be Gaussian in one limit and cubic exponential (same as [13]) in another (see section VI for a detailed discussion on this). The slow exponential decay of the linear response and the nonlinear echoes in the fluid regime is, on the other hand, in good agreement with the slow decay of the linear response in this regime obtained by [17-22]. The classical collisional theory of plasma echoes [13] has su ff ered from the absence of this slow exponential decay, a defect we remedy in this  \npaper. Furthermore, a key aspect of our analysis that di ff ers from [17-22] is that we go beyond linear theory, and perform a second-order treatment of the kinetic equations to characterize the collisional damping of the plasma echo.  \nWe organize the paper as follows. In section II, we discuss the perturbative response theory of a weakly collisional plasma. In section III, we briefly sketch the solution of the linearized Boltzmann-Poisson equations assuming an LB collision operator. Section IV describes the detailed behaviour of the continuum response along with a complementary analysis using the Langevin formalism. In section V, we describe the discrete mode response, and discuss the damping rates and oscillation frequencies of the various modes. In section VI, we compute the complete linear response of a weakly collisional plasma, which is a temporal convolution of the continuum response and the discrete mode response. In Section VII, using second-order response theory, we discuss the impact of weak collisions on the observability of the plasma echo in experiments, which is shown to be a manifestation of the continuum response. We summarize our findings and its possible implications for a variety of applications in section VIII.', 'II. PERTURBATIVE RESPONSE THEORY FOR WEAKLY COLLISIONAL PLASMA': 'A plasma is characterized by the DF or phase space ( x , v ) density of charged particles, f ( x , v , t ). Typically, the electrons are much less massive and thus more mobile than the ions, and are the dominant drivers of plasma oscillations. Hence, we shall restrict ourselves to the study of the evolution of the electron DF, with the ions forming a stationary, homogeneous and cold background. The general equations governing the evolution of a (non-magnetized) electrostatic plasma are the Boltzmann-Poisson equations. The Boltzmann equation,  \n∂ f ∂ t + v · ∇ f + e m ∇ v f · ∇ ( Φ+Φ P) = C GLYPH<2> f GLYPH<3> , (1)  \nis a conservation equation for the DF. Here, e is the magnitude of electric charge and m is the electron mass. C [ f ] denotes the collision operator, Φ P is the electrostatic potential of an external perturbation, and Φ is the self-potential of the plasma, which is sourced by the electron DF via the Poisson equation:  \n∇ 2 Φ = e ϵ 0 Z d 3 v f -ne ! , (2)  \nwith ne the equilibrium electron (or ion) density.  \nIf the strength of the perturber potential, Φ P, is smaller than σ 2 , where σ is the velocity dispersion of the unperturbed equilibrium plasma (typically characterized by a Maxwellian DF), then the perturbation in f can be expanded as a power series in the perturbation parameter, ε ∼ | Φ P | /σ 2 , i.e., f = f 0 + ε f 1 + ε 2 f 2 + ... . The evolution equations for the linear order perturbation in the DF, f 1 (which we shall henceforth refer to as the linear response), and that in the plasma potential, Φ 1, are given by the following linearized form of the Boltzmann-Poisson equations:  \n∂ f 1 ∂ t + v · ∇ f 1 = -e m ∇ v f 0 · ∇ ( Φ P + Φ 1) + C GLYPH<2> f 1 GLYPH<3> , ∇ 2 Φ 1 = e ϵ 0 Z d 3 v f 1 . (3)  \nSimilarly, the evolution equations for the second order perturbations, f 2 and Φ 2, are given by  \n∂ f 2 ∂ t + v · ∇ f 2 = -e m ∇ v f 1 · ∇ ( Φ P + Φ 1) -e m ∇ v f 0 · ∇ Φ 2 + C GLYPH<2> f 2 GLYPH<3> , ∇ 2 Φ 2 = e ϵ 0 Z d 3 v f 2 . (4)  \nWe discuss the evolution of the linear response in sections III, IV, V and VI. We discuss the second order response in section VII in the context of the plasma echo phenomenon.  \nWeassume that the collisional scattering of electrons is predominantly driven by weak encounters with other electrons and ions, i.e., collisions only cause small deflections from the unperturbed collisionless orbits. Under this assumption, the collision operator, C [ f 1], can be approximated as the following Fokker-Planck operator:  \nC [ f 1] ≈ 1 2 ∂ ∂χ i Dij ∂ f 1 ∂χ j ! , (5)  \nwhere χ = ( x , v ) denotes the canonical phase space coordinates, and Dij denote the di ff usion tensor ( i , j = 1(1)3).', 'III. SOLVING THE LINEARIZED EQUATIONS': 'We solve the linearized problem (equations [3]) under a set of simplifying assumptions. First, we assume that the unperturbed system is homogeneous and characterized by an equilibrium DF, f 0( v ). Second, we restrict ourselves to 1D perturbations. Third, we assume that the majority of collisions are weak, i.e., only cause slight deflections of the electrons from their orbits. This entails that the collision operator can be simplified into the Fokker Planck form. One such operator that gives rise to a Maxwellian DF in equilibrium, given by  \nf 0( v ) = ne √ 2 πσ exp " -v 2 2 σ 2 # , (6)  \nis attributed to Lenard and Bernstein (LB) [12], and is given by  \nC GLYPH<2> f GLYPH<3> = ν c ∂ ∂ v v f + σ 2 ∂ f ∂ v ! , (7)  \nwhere ν c is the collision frequency and σ is the velocity dispersion of the system when it is in equilibrium. The Boltzmann equation with the LB operator is also known as the  \nmean-field Kramers equation [see 23, for a detailed discussion of the initial-value problem for this equation in the context of weakly collisional plasmas and self-gravitating systems]. The fundamental assumption behind this form of the collision operator is that the collisions are su ffi ciently short-range and random. The LB operator conserves the particle number but not the linear momentum and kinetic energy during collisions. This shortcoming is well-known in the plasma literature [9-12]. Despite this shortcoming, we expect the scaling behaviours of di ff erent aspects of plasma relaxation to be qualitatively robust (see section IV B for details).  \nUnder the above assumptions, the linearized Boltzmann equation given by the first of equations (3) becomes the Lenard-Bernstein (LB) equation in 1D, and can be written as  \n∂ f 1 ∂ t + v ∂ f 1 ∂ x = -e m ∂ f 0 ∂ v ∂ ∂ x ( Φ 1 + Φ P) + ν c ∂ ∂ v v f 1 + σ 2 ∂ f 1 ∂ v ! . (8)  \nTo solve this equation simultaneously with the Poisson equation, we take the Fourier transform of both sides in x and the Laplace transform in t . We perform a Laplace transform (and not a Fourier transform) in time because we are interested in the initial-value problem. We then obtain  \n( γ + ikv ) ˜ f 1 k = f 1 k (0) -ik GLYPH<16> ˜ Φ 1 k + ˜ Φ k GLYPH<17> e m ∂ f 0 ∂ v + ν c ∂ ∂ v v ˜ f 1 k + σ 2 ∂ ˜ f 1 k ∂ v ! . (9)  \nHere, ˜ f 1 k , ˜ Φ 1 k and ˜ Φ k are the Fourier-Laplace transforms of f 1, Φ 1 and Φ P respectively, i.e.,  \nQ ( x , t ) = 1 2 π i Z c + i ∞ c -i ∞ d γ Z ∞ -∞ d k ˜ Qk ( γ ) exp GLYPH<2> ikx + γ t GLYPH<3> , (10)  \nwhere Q refers to each of f 1, Φ 1 and Φ P, while ˜ Qk refers to each of ˜ f 1 k , ˜ Φ 1 k and ˜ Φ k . Here, c is a real number such that the γ integral is performed within its region of convergence. ˜ Φ 1 k is related to ˜ f 1 k through the Poisson equation,  \n˜ Φ 1 k = -e ϵ 0 k 2 Z d v ˜ f 1 k . (11)  \nLet us now non-dimensionalize the LB and Poisson equations (equations [9] and [11]) before proceeding further. Defining κ = k σ/ν c, η = γ/ν c, u = v /σ , ˜ g 1 κ = ˜ f 1 k ( σν c / ne ), g 0 = f 0( σ/ ne ), h κ = f 1 k (0)( σ/ ne ), A κ ( η ) = -( ν c e / m σ 2 ) ˜ Φ k ( η ), and the plasma frequency, ω P = p nee 2 / m ϵ 0, we have the following non-dimensional form of the linearized LB and Poisson equations:  \n( η + i κ u ) ˜ g 1 κ = h κ + i κ A κ ( η ) ∂ g 0 ∂ u + i κ ω P ν c ! 2 ∂ g 0 ∂ u Z d u ˜ g 1 κ + ∂ ∂ u u ˜ g 1 κ + ∂ ˜ g 1 κ ∂ u ! . (12)  \nThe above is an integro-di ff erential equation and thus further simplification is required before a solution is attempted. This crucial simplifying step involves taking the velocity Fourier transform of ˜ g 1 κ , h κ and g 0 (assuming that these tend to zero su ffi ciently fast as u tends to ±∞ ), i.e., expanding these as  \n˜ g 1 κ ( u ) = Z ∞ -∞ d w exp [ iwu ] ˜ G 1 κ ( w ) , h κ ( u ) = Z ∞ -∞ d w exp [ iwu ] H κ ( w ) , g 0( u ) = Z ∞ -∞ d w exp [ iwu ] G 0( w ) . (13)  \nThis implies that R d u ˜ g 1 κ ( u ) = 2 π ˜ G 1 κ (0). This reduces the integro-di ff erential equation (12) to the following first order ordinary di ff erential equation in w :  \n( w -κ ) ∂ ˜ G 1 κ ∂ w + ( w 2 + η ) ˜ G 1 κ = H κ ( w ) -       2 π κ ω P ν c ! 2 ˜ G 1 κ (0) + κ A κ ( η )        wG 0( w ) . (14)  \nThis can be integrated using the method of integrating factor [12] to yield the following solution for ˜ G 1 κ ( w ):  \n˜ G 1 κ ( w ) = exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) × Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 ×        H κ ( w \' ) -       2 π κ ω P ν c ! 2 ˜ G 1 κ (0) + κ A κ ( η )        w \' G 0( w \' )        . (15)  \nSubstituting w = 0 in the above, we get the following expression for ˜ G 1 κ (0):  \n˜ G 1 κ (0) = N ( κ, η ) D ( κ, η ) . (16)  \nHere,  \nN ( κ, η ) = ( -κ ) -( η + κ 2 ) Z 0 κ d w \' exp " w \' 2 2 + κ w \' # × GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 GLYPH<2> H κ ( w \' ) -κ A κ ( η ) w \' G 0( w \' ) GLYPH<3> ,  \nD ( κ, η ) = 1 + 2 π κ ω P ν c ! 2 ( -κ ) -( η + κ 2 ) × Z 0 κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' G 0( w \' ) = 1 + ω P κν c ! 2 GLYPH<2> 1 -( a -ζ ) ζ -a exp GLYPH<2> ζ GLYPH<3> γ ( a , ζ ) GLYPH<3> , (17)  \nwhere a = η + κ 2 , ζ = κ 2 , and γ ( a , ζ ) = R ζ 0 d z z a -1 exp [ -z ] is the lower incomplete gamma function. We have assumed that the unperturbed DF, g 0, is of Maxwellian form, i.e., g 0( u ) = exp h -u 2 / 2 i. √ 2 π , whose velocity Fourier transform is given by  \nG 0( w ) = 1 2 π Z ∞ -∞ d u exp [ -iwu ] g 0( u ) = 1 2 π exp h -w 2 / 2 i . (18)  \nNote that the time-dependent potential response of the system is nothing but Φ 1 k = -(2 π nee /ϵ 0)( G 1 κ (0) / k 2 ), where G 1 κ (0) is the inverse Laplace transform of ˜ G 1 κ (0). We shall compute this in section VI.', 'IV. LINEAR RESPONSE: CONTINUUM': 'Computing the total response, g 1 κ ( u , t ), requires us to take the inverse Fourier transform in u and the inverse Laplace transform in t of equation (15). As detailed in Appendix A, first we substitute ˜ G 1 κ (0) from equation (16) in the expression for ˜ G 1 κ ( w ) given in equation (15), and then write it as a sum of three terms, ˜ G (1) 1 κ ( w ), ˜ G (2) 1 κ ( w ) and ˜ G (3) 1 κ ( w ). Performing the inverse Laplace transform of ˜ G 1 κ ( w ) thus yields g 1 κ ( u , t ), which can be written as a sum of three terms that are the inverse Laplace transforms of ˜ G (1) 1 κ ( w ), ˜ G (2) 1 κ ( w ) and ˜ G (3) 1 κ ( w ) respectively: a continuum response, a discrete mode response (temporally convoluted with the continuum response), and a direct response to the perturber, also known as the \'wake\' (temporally convoluted with the discrete modes and the continuum). The continuum response depends on the continuous velocities of the particles, hence the name. The discrete mode response consists of coherent oscillations (modes) of the entire system. These oscillations occur at discrete frequencies, η n ( n = 0 , 1 , 2 , ... ), which satisfy the dispersion relation, D ( κ, η n ) = 0. The wake response has a temporal dependence that is similar to that of the external perturber.  \nLet us first compute the continuum part of the response, g cont 1 κ ( u , t ), which amounts to computing the inverse Laplace transform of ˜ G (1) 1 κ ( w ) (first of equations [A2]). Assuming that the unperturbed DF is Maxwellian, this yields (see Appendix A for details):  \ng cont 1 κ ( u , t ) = r 2 1 -e -2 ν c t exp " -u 2 2 GLYPH<0> 1 -e -2 ν c t GLYPH<1> # × exp GLYPH<20> -κ 2 GLYPH<18> ν c t -2 tanh ν c t 2 GLYPH<19>GLYPH<21> × exp GLYPH<20> -i κ u tanh ν c t 2 GLYPH<21> × F ( κ, ν c t ) , (19)  \nwhere  \nF ( κ, ν c t )  \n= Z ∞ -∞ d z exp         -       z +        r 2 tanh ν c t 2 κ -i u p 2 GLYPH<0> 1 -e -2 ν c t GLYPH<1>               2         × H κ GLYPH<16> κ + ze -ν c t GLYPH<17> . (20)  \nHere, H κ is the velocity Fourier transform of the initial DF perturbation, h κ .', 'A. Di ff erent models for h κ ( u )': 'Evidently, the temporal behaviour of the continuum response depends on the functional form of the initial perturbation, h κ ( u ). We compute the response for a few specific models of h κ ( u ).', '1. Gaussian pulse': 'Let us consider the case of an initial perturbation in the DF that is Gaussian or Maxwellian in velocities, but with a velocity dispersion di ff erent from the equilibrium velocity dispersion σ :  \nh κ ( u ) = h κ 0 √ 2 π s 2 exp " -( u -u 0) 2 2 s 2 # . (21)  \nHere, u 0 is the drift velocity of the initial pulse, and σ 0 is its velocity dispersion, with s = σ 0 /σ . The velocity Fourier transform of h κ ( u ) is given by  \nH κ ( w ) = h κ 0 2 π exp " -s 2 w 2 2 -iu 0 w # . (22)  \nIn this case, the continuum response in equation (19) becomes  \ng cont 1 κ ( u , t ) = h κ 0 p 2 πδ 2 u ( t ) exp       -GLYPH<0> u -u 0 e -ν c t GLYPH<1> 2 2 δ 2 u ( t )       × exp h -κ 2 δ 2 x ( t ) i × exp h -i κ GLYPH<16> ζ (1) x ( t ) u + ζ (2) x ( t ) u 0 GLYPH<17>i , (23)  \nwhere  \nζ (1) x ( t ) = GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<16> 1 -GLYPH<16> 1 -s 2 GLYPH<17> e -ν c t GLYPH<17> 1 -GLYPH<0> 1 -s 2 GLYPH<1> e -2 ν c t , ζ (2) x ( t ) = GLYPH<0> 1 -e -ν c t GLYPH<1> 2 1 -GLYPH<0> 1 -s 2 GLYPH<1> e -2 ν c t , δ 2 u ( t ) = 1 -GLYPH<16> 1 -s 2 GLYPH<17> e -2 ν c t , (24)  \nδ 2 x ( t ) = ν c t -1 -e -ν c t 2 3 -e -ν c t + GLYPH<16> 1 -s 2 GLYPH<17> GLYPH<0> 1 -3 e -ν c t GLYPH<1> 1 -GLYPH<0> 1 -s 2 GLYPH<1> e -2 ν c t .  \nFIG. 1. The continuum response of a weakly collisional plasma as a function of time for ν c /ω P = 0 . 01 , 0 . 1 and 1 indicated by the dark, medium dark and light blue lines, and three di ff erent values of k σ in the three panels. The solid lines indicate the response when the initial perturbation in the DF is of the equilibrium Maxwellian form ( s = σ 0 /σ = 1). The dashed lines denote the response for an initial perturbation that is Dirac delta in the velocities ( s = 0). The dotted lines indicate the continuum response predicted by [13]. Note that the s = 1 and s = 0 cases are very similar for ν c > k σ , but di ff erent, albeit with the same temporal scaling, otherwise. The continuum response decays as ∼ exp h -( k σ ) 2 ν c t 3 / 3 i for the s = 1 case and ∼ exp h -( k σ ) 2 ν c t 3 / 12 i for the s = 0 case at t < 1 /ν c, but switches over to a much slower ∼ exp h -GLYPH<16> ( k σ ) 2 /ν c GLYPH<17> t i decay at t > 1 /ν c, unlike [13] who only predict the former cubic exponential decay. Note that the collisional damping is most pronounced for ν c ∼ k σ , or in other words, when the perturbation scale, λ ∼ 1 / k , is comparable to the mean free path of collisions, λ c = σ/ν c.  \n<!-- image -->  \nLet us analyze separately the key features of the above response. The first factor in equation (23) denotes the damping and broadening of the initial perturbation in the velocity space due to the collisional di ff usion of particle velocities. The perturbation starts out being a Gaussian (in u ) with a velocity dispersion s (in units of σ ), then gradually damps away and spreads out to attain the equilibrium Maxwellian form with  \nvelocity dispersion σ .  \nThe second factor in equation (23) indicates the damping of the response (in κ ) due to collisions. The extent of damping depends on the wave-number k of the perturbation. The asymptotic temporal behaviour of this damping factor for the s , 0 case (we shall discuss the s = 0 case, which corresponds to an initial Dirac delta pulse, separately) is as follows:  \nexp h -κ 2 δ 2 x ( t ) i ≈        exp h -κ 2 ( ν c t ) 3 / 3 i = exp h -( k σ ) 2 ν c t 3 / 3 i , ν c t ≪ 1 , exp h -κ 2 ν c t i = exp h -GLYPH<16> ( k σ ) 2 /ν c GLYPH<17> t i , ν c t ≫ 1 . (25)  \nThe small time ( t ≪ 1 /ν c) behaviour of the above collisional damping factor is similar to the classic [13] result. This damping occurs cubic exponentially over a timescale, τ d , 1 = GLYPH<16> 3 /ν c( k σ ) 2 GLYPH<17> 1 / 3 . However, the long time ( t ≫ 1 /ν c) damping occurs exponentially, much slower than the [13] prediction. The timescale of this late time decay is τ d , 2 = ν c / ( k σ ) 2 . When ν c ∼ k σ , or in other words the scale of the perturbation, λ ∼ 1 / k , is comparable to the mean free path, λ c = σ/ν c, τ d , 1 is comparable to τ d , 2, and the corresponding response damps the most. This is evident in Fig. 1, which plots the amplitude of the continuum response as a function of time for three di ff erent values of ν c and three di ff erent values of k σ as indicated. Our result, computed using equation [23], is shown by the solid and dashed lines for the s = 1 and 0 cases respectively, and the corresponding [13] prediction is shown by dotted lines. The notable di ff erence between our prediction and that of [13] is a consequence of the fact that they neglected  \nthe damping / drag term in the LB operator, i.e., the first term on the RHS of equation (8). This approximation is justified only for t ≲ 1 /ν c, when collisions cause the particle velocities to di ff use / disperse rather than damp on an average. This is the di ff usive regime of collisions. However, collisions not only disperse the velocities of the particles but also damp out their individual mean velocities over time. This damping regime, which takes over at t ≳ 1 /ν c, is described by the first term on the RHS of equation (8), the very term that [13] neglect. By neglecting the damping term, [13] have limited themselves to the di ff usive regime of collisions, and therefore missed out on the crucial switchover of the decay to a slow exponential one in the long time damping regime. To obtain more physical insight into this, we shall discuss the collisional dynamics in the di ff usive and damping regimes in detail, in section IV B, using a stochastic analysis with the Langevin equation, which is complementary to the Fokker-Planck approach taken so far.  \nFIG. 2. The continuum response of a weakly collisional plasma as a function of x (in units of σ/ν c). The solid lines denote f 1( x , v , t ) (computed using equation [26] and normalized by its maximum value at t = 0) at seven di ff erent values of ν c t (ranging from 0, indicated by the blue / uppermost lines, to 3, indicated by the red / lowermost lines), while the dashed lines denote ρ 1( x , t ) (computed using equation [27] and normalized by its maximum value at t = 0). We adopt σ x , 0 = 0 . 5 σ/ν c, v = 0 . 5 σ , v 0 = σ and s = σ 0 /σ = 0 . 9, i.e., the initial DF perturbation is close to the equilibrium Maxwellian form. Both quantities have been normalized by their maximum value at t = 0. Note that both perturbations decay, broaden and slow down in x due to collisional di ff usion followed by damping, but the density perturbation decays faster due to the additional e ff ect of phase-mixing.  \n<!-- image -->  \nThe third factor in equation (23) is a phase factor that denotes the free streaming of the particles. It can be seen that  \nf cont 1 ( x , v , t ) = ne σ Z ∞ d k exp [ ikx ] g cont 1 κ ( u , t  \nat small time ( ν c t ≪ 1), this factor evolves as exp [ -i κ u ν c t ] = exp [ -ikvt ]. This behaviour is the same as that in the collisionless ( ν c → 0) case. At large time ( ν c t ≫ 1), the phase factor becomes exp [ -i κ ( u + u 0)], and therefore time-independent. This implies that the pulse comes to a halt after roughly the collision time, 1 /ν c, owing to the constant drainage of linear momentum from a particle via collisions. This is a generic feature of collisional damping in Fokker-Planck dynamics.  \nSo far we have seen that each k mode of the continuum response damps away and decelerates due to collisional diffusion and damping. The spatial response is an integrated response of all k modes. Since each k mode damps away at a rate that depends non-linearly on k , the spatial response, which is the sum of all k modes, will deviate from its initial functional form over time, i.e., the weakly collisional plasma acts as a dispersive medium. We shall demonstrate this by computing the spatial response for a given spatial distribution of the initial perturbation, i.e., a given h κ 0. Since the response for each k mode is Gaussian in k , as evident from equation (23) (for a Gaussian velocity distribution of the initial perturbation as we consider here), the integration over k is straightforward for a Gaussian (in k ) functional form of h κ 0, i.e., a Gaussian (in x ) functional form of the initial perturbation. The spatial response can be obtained by taking the inverse spatial Fourier transform of equation (23). For an initial Gaussian pulse with spatial dispersion σ x , 0, i.e., g 1( x , u , 0) = GLYPH<16> F / √ 2 π GLYPH<17> exp h -x 2 / 2 σ 2 x , 0 i exp h -( u -u 0) 2 / 2 s 2 i , for which h κ 0 = ( F / 2 π ) exp h -k 2 σ 2 x , 0 / 2 i exp h -( u -u 0) 2 / 2 s 2 i , we take the inverse Fourier transform of equation (23) and reestablish the dimensions to obtain the following form for the response:  \n-∞ ) = neF p 2 πσ 2 δ 2 u ( t ) exp       -GLYPH<0> v -v 0 e -ν c t GLYPH<1> 2 2 σ 2 δ 2 u ( t )       × 1 q 2 π GLYPH<16> σ 2 x , 0 + 2( σ/ν c) 2 δ 2 x ( t ) GLYPH<17> exp           -h x -1 ν c GLYPH<16> ζ (1) x ( t ) v + ζ (2) x ( t ) v 0 GLYPH<17>i 2 2 GLYPH<16> σ 2 x , 0 + 2 ( σ/ν c) 2 δ 2 x ( t ) GLYPH<17>           . (26)  \nConservation of particle number ensures that the pulse damps out as well as broadens due to collisions, as manifested by both x and v parts of the above response. The spatial broadening of the pulse occurs as ∼ GLYPH<16> 2 σ 2 / 3 ν 2 c GLYPH<17> ( ν c t ) 3 for t ≲ 1 /ν c, and as ∼ GLYPH<16> 2 σ 2 /ν 2 c GLYPH<17> ν c t for t ≳ 1 /ν c, while the spatial damping occurs as the inverse of this. If the width of the initial pulse is small, i.e., σ x , 0 ≪ σ/ν c, it damps and broadens in the former way for a long time, i.e., we are in the di ff usive regime. The typical damping timescale in this regime is ∼ GLYPH<16> 3 σ 2 x , 0 / 2 σ 2 ν c GLYPH<17> 1 / 3 . On the other hand, a well-extended initial pulse with σ x , 0 ≫ σ/ν c, falls in the damping regime, i.e., it quickly transitions to the latter phase of exponential damping  \nand broadening, which occurs over a relatively long timescale of ∼ GLYPH<16> σ 2 x , 0 / 2 σ 2 GLYPH<17> ν c. Since both timescales increase with σ x , 0 (for fixed ν c and σ ), we see that small-scale perturbations collisionally damp away faster than large-scale ones. Moreover, since the former timescale decreases with ν c (for fixed σ x , 0 and σ ) and the latter increases, perturbations damp slowly in both small and large ν c limits. The damping and broadening of the pulse due to collisions is visible in Fig. 2, which plots f 1 (solid lines) as a function of x at di ff erent points of time as indicated. Note that the pulse slows down over time due to collisional damping. The three key features of Fokker-Planck collisional dynamics are the damping, broadening and deceleration of a pulse, all three of which are manifest in Fig. 2.  \nLet us now study the evolution of the density perturbation, ρ cont 1 ( x , t ) = R d v f cont 1 ( x , v , t ). Upon integrating equation (26) over velocity, we obtain  \nρ cont 1 ( x , t ) = neF r 2 π GLYPH<20> σ 2 x , 0 + GLYPH<16> σ ν c GLYPH<17> 2 GLYPH<26> 2 δ 2 x ( t ) + ζ (1) x 2 ( t ) δ 2 u ( t ) GLYPH<27>GLYPH<21> × exp             -GLYPH<16> x -v 0 ν c GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<17> 2 2 GLYPH<20> σ 2 x , 0 + GLYPH<16> σ ν c GLYPH<17> 2 GLYPH<26> 2 δ 2 x ( t ) + ζ (1) x 2 ( t ) δ 2 u ( t ) GLYPH<27>GLYPH<21>             . (27)  \nThe density perturbation initially propagates at velocity v 0 and comes to a halt after t ∼ 1 /ν c. At the same time, it damps and spreads out. This is not only due to collisional damping but also phase-mixing, as long as the initial perturbation has some finite spread in velocities, or in other words particles are initially perturbed over a finite velocity range. The width of the Gaussian scales as ∼ q σ 2 x , 0 + σ 2 t 2 at small t , due to phase-mixing. At intermediate time, the width scales as ∼ q σ 2 x , 0 + σ 2 ν c t 3 / 3, which is an outcome of the collisional di ff usion of particle velocities. At late time, it scales as ∼ q σ 2 x , 0 + 2 GLYPH<0> σ 2 /ν c GLYPH<1> t , due to collisional damping. Due to  \nthe additional contribution of phase-mixing to the damping besides collisions, the density perturbation ρ cont 1 damps away faster than the fine-grained DF f cont 1 , as manifest in Fig. 2 that plots f cont 1 (solid lines) and ρ cont 1 (dashed lines) as functions of x .', '2. Dirac delta pulse': 'We now evaluate the response in the s → 0 limit of the initial Gaussian pulse. This corresponds to an initial pulse that is Dirac delta in the particle velocities, i.e., h κ ( u ) = h κ 0 δ ( u -u 0), whose Fourier transform is H κ ( w ) = h κ 0 exp [ -iu 0 w ] / (2 π ). In this case, the response is given by equation (23), with  \nζ (1) x ( t ) = ζ (2) x ( t ) = tanh ν c t 2 , δ 2 u ( t ) = 1 -e -2 ν c t , δ 2 x ( t ) = ν c t -2 tanh ν c t 2 . (28)  \nThe asymptotic temporal behaviour of the collisional damping factor is:  \nexp h -κ 2 δ 2 x ( t ) i ≈        exp h -κ 2 ( ν c t ) 3 / 12 i = exp h -( k σ ) 2 ν c t 3 / 12 i , ν c t ≪ 1 exp h -κ 2 ν c t i = exp h -GLYPH<16> ( k σ ) 2 /ν c GLYPH<17> t i , ν c t ≫ 1 . (29)  \nNote that the small time behaviour of the damping is the same as that in the general case given by equation (25), expect for the numerical factor of the exponent. The long time behaviour is exactly the same as equation (25).  \nThe phase factor, exp h -i κ ( ζ (1) ) x ( t ) u + ζ (2) ) x ( t ) u 0 i , scales as ∼ exp [ -i κ ( u + u 0) t / 2] = exp [ -ik ( v + v 0) t / 2] for ν c t → 0. In this limit, the velocity part of the response tends to δ ( v -v 0), and thus the phase-factor goes to exp [ -ikvt ]. This reproduces the form of the continuum response in the collisionless limit.', '3. Equilibrium Gaussian pulse': 'In the s → 1 limit, the initial perturbation, h κ ( u ), is of the equilibrium Maxwellian form, i.e., h κ ( u ) = h κ 0 exp h -( u -u 0) 2 / 2 i / √ 2 π , and its Fourier transform is equal to H κ ( w ) = h κ 0 exp h -w 2 / 2 -iu 0 w i / (2 π ). In this case, the response is given by equation (23), with  \nζ (1) x ( t ) = 1 -e -ν c t , ζ (2) x ( t ) = GLYPH<16> 1 -e -ν c t GLYPH<17> 2 , δ u ( t ) = 1 , δ 2 x ( t ) = ν c t -GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<0> 3 -e -ν c t GLYPH<1> 2 . (30)  \nThe phase factor is exp GLYPH<2> -ik GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<0> v + v 0 GLYPH<0> 1 -e -ν c t GLYPH<1>GLYPH<1>GLYPH<3> , which scales as exp [ -ikvt ] in the ν c t → 0 limit. This recovers the form of the continuum response in the limit of zero collisions. The perturbation is initially of the equilibrium Maxwellian form and therefore remains so throughout, since the Maxwellian DF is an equilibrium solution of the Boltzmann-Poisson equations with the LB collision operator. The amplitude of the response, however, damps away due to collisional di ff usion and damping. The damping factor has the same asymptotic behaviour as in the general case given by equation (25), i.e., the response damps away cubic exponentially till t ∼ 1 /ν c and transitions to a slower exponential decay afterwards.', 'B. Langevin approach: di ff usive and friction regimes': 'So far, we have seen that the collisional relaxation of a plasma occurs in two distinct regimes: the di ff usive and friction regimes. The nature of relaxation in these two regimes is very di ff erent. These distinct behaviours of plasma relaxation can be understood by taking a closer look at the stochastic dynamics of collisions. We assume that the collisions are typically short-range random encounters, and therefore mimic a stochastic process that randomly changes the phase space coordinates of the particles. If the collisions are random enough, they erase the dynamical memory (initial conditions) of the colliding particles over time, i.e., the collisional scatterings can be thought of as Markovian random walks of the particles in phase space. A key consequence of this random (or chaotic) nature of collisions is the phenomenon of molecular chaos , which refers to the fact that the momenta of two colliding particles, despite being correlated immediately after a collision, become uncorrelated soon enough due to their encounters with other particles 2 . This justifies the use of the Boltzmann equation as the governing kinetic equation. This also ensures that the collisions are approximately uncorrelated in time and can be modelled by a stochastic acceleration with a white noise power spectrum. Hence, in the same spirit as that of modeling Brownian motion [24, 25], we can model the collisional dynamics in a plasma using the Langevin equation, which is nothing but Newton\'s equation of motion with damping and stochastic forces:  \nd 2 x d t 2 + ν c d x d t = ξ ( t ) . (31)  \nHere, -ν cd x / d t is the collisional damping / deceleration, also known as friction, and ξ ( t ) is the stochastic acceleration due to collisions. We assume that ξ ( t ) has a zero ensemble average and a white noise or uncorrelated power spectrum, i.e.,  \n⟨ ξ ( t ) ⟩ = 0 , GLYPH<10> ξ ( t ) ξ ( t \' ) GLYPH<11> = D δ ( t -t \' ) , (32)  \nwhere D is a constant with the dimensions of velocity squared per unit time.  \nThe Langevin equation is a stochastic di ff erential equation which evolves the trajectory of each individual particle. The ensemble averages of the evolving phase space quantities and their dispersions, obtained by solving the Langevin equation, are the same as the first and second moments of these quantities, obtained with respect to the DF that evolves via the Fokker-Planck equation, as long as the initial DF perturbation is su ffi ciently close to the equilibrium Maxwellian form. In fact the Fokker-Planck collision operator corresponding to the  \nLangevin equation given in equation (31) is the LB operator. The correspondence between the Langevin and Fokker-Planck equations therefore justifies that the results obtained in this section are an alternative way to obtain the various scalings of the response that are computed using the Fokker-Planck (LB) equation in the previous sections.  \nBefore solving the Langevin equation in its full glory, let us first take a careful look at the applicability of this equation in modeling plasma collisions. The key assumption is that the collisions can be modelled as short-range random encounters, despite the fact that Coulomb interactions between charged particles are long-range. The impact parameter of effective collisions is typically the Debye length, making the collisions functionally short-range. Hence, as long as collisional relaxation is dominated by scales larger than the Debye length ( ω P ≳ k σ ), as is typically the case for a weakly collisional plasma, the collisions can be modelled as shortrange stochastic encounters, and the Langevin (equivalently, the Fokker-Planck) description in this paper should properly capture the essential physics of plasma collisions.  \nWith the aforementioned assumptions, let us now solve the Langevin equation to evaluate the mean trajectories of particles and the position and velocity dispersions around them. Integrating equation (31) once with respect to t , with the initial condition, v (0) = v 0, we obtain the following form for the velocity:  \nv ( t ) = d x d t = v 0 e -ν c t + Z t 0 d t \' e -ν c( t -t \' ) ξ ( t \' ) . (33)  \nSince the ensemble average of ξ ( t ) is zero, the ensemble averaged or mean velocity is equal to  \n⟨ v ( t ) ⟩ = v 0 e -ν c t . (34)  \nThe exponential decay of the mean velocity manifests damping or friction: the constant drainage of momentum from each particle due to collisions.  \nThe velocity dispersion is given by (see Appendix B for detailed derivation)  \nσ 2 ( t ) = GLYPH<28> GLYPH<16> v ( t ) -v 0 e -ν c t GLYPH<17> 2 GLYPH<29> = D ν c GLYPH<16> 1 -e -2 ν c t GLYPH<17> , (35)  \nwhere we have used the uncorrelated nature of the stochastic force, i.e., the second of equations (32). Note that the velocity dispersion initially increases linearly with time, i.e., σ 2 ( t ) ∼ 2 D t . This is the essence of velocity di ff usion or particle random walk in the velocity space, and is caused by the stochastic acceleration, ξ ( t ). The late-time dynamics, which sets in after t ∼ 1 /ν c, is governed by the frictional deceleration, -ν cd x / d t . This friction / damping drives σ 2 towards the equilibrium value, lim t →∞ σ 2 ( t ) = σ 2 . Thus, from the above equation, we infer that D = ν c σ 2 , which is the celebrated Einstein relation. D is the velocity di ff usion coe ffi cient, which  \nis related to the damping coe ffi cient or collision frequency, ν c, by the proportionality factor, σ 2 . Therefore, friction and di ff usion go hand-in-hand during collisions, and one cannot exist without the other. Consequently, neglecting the friction term in the collision operator [13, 14] provides an incomplete description of collisional dynamics, which may lead to erroneous conclusions about the long time decay of the response in weakly collisional plasmas.  \nLet us now evaluate the time dependence of the particle position. Integrating equation (33) with respect to t , with the initial condition that x (0) = 0, we obtain  \nx ( t ) = v 0 ν c GLYPH<16> 1 -e -ν c t GLYPH<17> + Z t 0 d t \' Z t \' 0 d t \'\' e -ν c( t \' -t \'\' ) ξ ( t \'\' ) . (36)  \nRecalling that ⟨ ξ ( t ) ⟩ = 0, we see that the mean position is equal to  \n⟨ x ( t ) ⟩ = v 0 ν c GLYPH<16> 1 -e -ν c t GLYPH<17> , (37)  \ni.e., the particle gradually slows down and traverses an overall distance of v 0 /ν c. The mean-squared position is given by (see Appendix B for detailed derivation)  \nσ 2 x ( t ) = *" x ( t ) -v 0 ν c GLYPH<16> 1 -e -ν c t GLYPH<17> # 2 + = σ ν c ! 2 " ν c t -GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<0> 3 -e -ν c t GLYPH<1> 2 # , (38)  \nwhere we have used the second of equations (32), i.e., ξ ( t ) is a white noise source. Initially, in the di ff usive regime, where the dynamics is predominantly driven by the stochastic acceleration ξ ( t ), the mean-squared position scales as σ 2 x ( t ) ∼ ν c σ 2 t 3 / 3. This is a consequence of the fact that the velocity dispersion, σ 2 ≈ d 2 σ 2 x / d t 2 , scales as ∼ 2 D t = 2 σ 2 ν c t in this regime. Finally, after t ∼ 1 /ν c, the mean-squared position scales as σ 2 x ( t ) ∼ ( σ 2 /ν c) t . This owes to the fact that the velocity dispersion, σ 2 ( t ) ≈ ν c d σ 2 x / d t , attains a nearly constant value, σ 2 , due to friction. In this late-time friction-dominated regime, the velocity evolution is governed by the frictional deceleration, but the evolution of the mean-squared position is still driven by the stochastic acceleration, ξ ( t ).  \nThe temporal dependence of the mean and mean squared velocity and position obtained from the above analysis is very similar to what shows up in our analysis of plasma relaxation using the LB operator in the previous sections. The exponentially decaying drift velocity obtained from the Langevin analysis (equation [34]), is the same as the drift velocity of the pulse obtained using the LB analysis (see equations [23] and 26). The slowing drift of the mean position is manifest in both the Langevin (equation [37]) and LB (equation [27]) analyses. The time-dependence of the velocity dispersion is also the same; compare σ 2 ( t ) in equation (35) with σ 2 δ 2 u ( t ) in equation (26), where δ 2 u ( t ) is given by the second of equations (28). And the mean squared position, computed using the Langevin  \nequation (equation [38]) is the same as the spatial width of the response computed using the LB analysis (note δ 2 x ( t ) in equation [30]). The prefactors of the above quantities slightly di ff er between the LB and Langevin analyses when the initial perturbation in the DF is very di ff erent from the equilibrium Maxwellian form. This is because the di ff usion coe ffi cients of the LB operator rely on the Einstein relation, and are therefore dictated by the equilibrium Maxwellian DF. The correspondence is exact only when the initial pulse is of the equilibrium Maxwellian form in velocities. Still, the asymptotic temporal behaviour of ⟨ x ⟩ , ⟨ v ⟩ , σ 2 x and σ 2 obtained using the Langevin analysis is the same as that obtained using the LB analysis in all cases. This proves that the LB operator, despite being only an approximate description of collisional dynamics, properly captures the essential physics of collisions, particularly in the small and large time limits. Put di ff erently, the di ff usive and friction regimes of collisional relaxation predicted by the LB operator in section IV A are not merely artefacts of the form of the collision operator chosen, but are qualitatively robust phenomena in a system governed by weak collisions.', 'V. LINEAR RESPONSE: DISCRETE MODES': 'So far, we have studied the first part of the plasma response, the continuum, which encapsulates the collisionally damped streaming of particles. The second part of the response on the other hand consists of discrete modes 3 that are collective coherent oscillations of the system. It can be computed by taking the inverse Laplace transform of ˜ G (2) 1 κ ( ω ) in equation (A2). The resulting response involves a convolution of the continuum and discrete mode responses. We defer the computation of the combined response to the next section, and only discuss the behaviour of the discrete modes in this section. The discrete modes oscillate at frequencies that follow the dispersion relation,  \nD ( κ, η ) = 1 + ω P κν c ! 2 GLYPH<2> 1 -( a -ζ ) ζ -a exp GLYPH<2> ζ GLYPH<3> γ ( a , ζ ) GLYPH<3> = 0 , (39)  \nwhere a = η + κ 2 , ζ = κ 2 , and γ ( a , ζ ) = R ζ 0 d z z a -1 exp [ -z ] is the lower incomplete gamma function. These frequencies are in general complex. The numerical computation of these frequencies is quite involved and the standard root-finding algorithms are found to be computationally expensive and imprecise. Therefore, we develop a novel root-finding algorithm to accomplish this. We rewrite the above dispersion relation as the following equations:  \nexp GLYPH<2> ζ GLYPH<3> Re GLYPH<2> ( a -ζ ) ζ -a γ ( a , ζ ) GLYPH<3> = 1 + κν c ω P ! 2 , exp GLYPH<2> ζ GLYPH<3> Im GLYPH<2> ( a -ζ ) ζ -a γ ( a , ζ ) GLYPH<3> = 0 , (40)  \nFIG. 3. The damping time, t d (in units of 1 / k σ ), of the discrete modes and the continuum as a function of ν c / k σ . The dark, medium dark and light blue dots indicate the discrete modes for ω P / k σ = 0 . 5 , 2 . 5 and 4 . 5 respectively. The white dashed lines denote the discrete modes of the Lenard-Bernstein (LB) branch. The rest belong to the Landau branch. The solid and dashed lines denote the damping time of the continuum response at t = 0 . 5 and 2 / k σ respectively. Note that at small time ( t ≲ 1 /ν c), the continuum response damps slower than the discrete modes. At later times, the discrete modes damp slower than the continuum for ν c < k σ . For ν c > k σ , the continuum response damps at the same rate as the zero-th order LB mode. The discrete modes of the Landau branch smoothly tend to the Landau modes in the ν c / k σ → 0 limit. At ν c ≈ 2 ω P, the Landau branch modes bifurcate and transition from being oscillatory (for ν c ≲ 2 ω P) to non-oscillatory (for ν c ≳ 2 ω P).  \n<!-- image -->  \nwhere Re and Im denote real and imaginary parts respectively. These equations are non-linear in Re( a ) and Im( a ), i.e., in Re( η ) and Im( η ). We simultaneously solve these non-linear equations in the following way. For each κ and ω P /ν c, we compute two sets of contours in the Re( η ) -Im( η ) plane. The first set of contours corresponds to those values of Re( η ) and Im( η ) that satisfy the first of the above equations. The second set corresponds to those values of Re( η ) and Im( η ) that satisfy the second equation. Now, we compute the intersection points of the first and second sets of contours, which are nothing but the simultaneous solutions of the above equations in the Re( η )-Im( η ) space. The η thus obtained, follow the dispersion relation given by equation (39), and are therefore the discrete mode frequencies.  \nThe discrete mode spectrum of a weakly collisional plasma has a rich structure. As shown by [12, 26, etc.], all modes in a plasma that is perturbed around an f 0 with ∂ f 0 /∂ v < 0 (e.g., the Maxwellian DF) are damped, i.e., have Re( γ ) < 0. The damping time is given by t d = 1 / | Re( γ ) | . Fig. 3 plots this damping time (in units of 1 / k σ ), as a function of ν c / k σ . The di ff erent colors denote di ff erent values of ω P / k σ as indicated. The white dashed lines trace out a particular family of  \nmodes called the Lenard-Bernstein or LB modes, which are non-oscillating (discussed in detail in section V B). The black solid and dashed lines indicate the damping time of the continuum response, defined as GLYPH<12> GLYPH<12> GLYPH<12> ∂ ln g cont 1 κ /∂ t GLYPH<12> GLYPH<12> GLYPH<12> -1 , at t = 0 . 5 and 2 times 1 / k σ , respectively. Since the continuum response scales as ∼ exp h -( k σ ) 2 ν c t 3 / 3 i at small time, and the modal response scales as ∼ exp GLYPH<2> -γ t GLYPH<3> throughout, the continuum response dominates over the discrete modes at short times. At long times, the discrete modes dominate on scales smaller than the mean free path ( k > ν c /σ ). However, on scales larger than the mean free path ( k < ν c /σ ), the dominant modes are the continuum as well as the least damped discrete mode, which is one of the non-oscillating LB modes (see section V B for details).  \nThe discrete modes in a weakly collisional plasma are diverse in nature and deserve a systematic categorization. We identify two distinct classes of discrete modes: (i) the Landau branch, which analytically continues to the Landau modes in the collisionless limit ( ν c ≪ k σ ) and are Landau damped, and (ii) the Lenard-Bernstein (LB) branch, which consists of an extra set of modes that appear in addition to the Landau-like modes for non-zero ν c. Wenow discuss these two families of discrete modes in detail.', 'A. Landau branch': 'The modes belonging to the Landau branch analytically continue to the Landau modes in the collisionless ( ν c / k σ → 0) limit. As evident from Fig. 3, the damping timescales of the Landau branch modes tend towards the Landau values in the ν c / k σ → 0 limit, unlike the LB branch modes, indicated by the white dashed lines, whose damping times tend to zero. The analytic continuation to the Landau values can also be seen by rewriting the dispersion relation given in equation (39) in the following form [12]:  \nD ( k , γ ) = 1 + ω 2 P ν c Z ∞ 0 d s GLYPH<0> 1 -exp [ -ν c s ] GLYPH<1> × exp " -γ + k 2 σ 2 ν c ! s + k 2 σ 2 ν 2 c GLYPH<0> 1 -e -ν c s GLYPH<1> # = 0 . (41)  \nIn the ν c / k σ → 0 limit, this reduces to the Landau dispersion relation [26],  \nD ( k , γ ) = 1 + ω 2 P Z ∞ 0 d s s exp " -γ s + k 2 σ 2 s 2 2 !# = 0 . (42)  \nThe solutions of this equation are the Landau modes, γ 0, which are complex. For each k , there are infinitely many discrete modes 4 . Of these, the least damped modes have the following asymptotic form [26]:  \nγ 0 = -2 k σ p ln ( k σ/ω P) ± i π 2 k σ √ ln ( k σ/ω P) , k σ ≫ ω P , = -r π 8 k σ GLYPH<18> ω P k σ GLYPH<19> 4 exp " -1 2 GLYPH<18> ω P k σ GLYPH<19> 2 # ± i ω P        1 + 3 2 k σ ω P ! 2        , k σ ≪ ω P . (43)  \nThe above form shows that, in the collisionless limit, the modes are oscillatory and very weakly damped on scales much larger than the Debye scale ( k ≪ k D = ω P /σ ). These are also known as Langmuir modes. On scales smaller than the Debye scale ( k ≫ k D), the modes are oscillatory but strongly (Landau) damped. Hence, oscillations in a collisionless plasma are damped by collective excitations only within the Debye scale.  \nThe Landau modes are modified in presence of collisions, i.e., for non-zero ν c. Let us compute the leading order correction to γ in presence of collisions, using linear perturbation theory in the neighborhood of the Landau modes, γ 0 (taking ν c ≪ min GLYPH<2> γ 0 , k σ GLYPH<3> ). We expand γ about γ 0 as γ = γ 0 + ν c δ 1. Then, we rewrite the dispersion relation given in equation (41) by only retaining terms that are zero-th and linear order in ν c, as follows:  \nD ( k , γ 0 + ν c δ 1) = 1 + ω 2 P Z ∞ 0 d s " s -ν c s 2 δ 1 + 1 2 -k 2 σ 2 s 2 6 !# × exp " -γ 0 s + k 2 σ 2 s 2 2 !# = 0 . (44)  \nUsing that γ 0 follow the Landau dispersion relation (equation [42]), and performing the s integrals, we obtain:  \nδ 1 = γ -γ 0 ν c = -1 2 + 1 6 P ( γ 0 / k σ ) Q ( γ 0 / k σ ) , (45)  \nwhere  \nP GLYPH<18> γ 0 k σ GLYPH<19> = " GLYPH<18> γ 0 k σ GLYPH<19> 3 + 4 γ 0 k σ # + " GLYPH<18> γ 0 k σ GLYPH<19> 4 + 6 GLYPH<18> γ 0 k σ GLYPH<19> 2 + 3 # r π 2 GLYPH<20> 1 + erf GLYPH<18> γ 0 k σ GLYPH<19>GLYPH<21> exp " 1 2 GLYPH<18> γ 0 k σ GLYPH<19> 2 # , Q GLYPH<18> γ 0 k σ GLYPH<19> = γ 0 k σ + " GLYPH<18> γ 0 k σ GLYPH<19> 2 + 1 # r π 2 GLYPH<20> 1 + erf GLYPH<18> γ 0 k σ GLYPH<19>GLYPH<21> exp " 1 2 GLYPH<18> γ 0 k σ GLYPH<19> 2 # . (46)  \nOn scales much larger than the Debye length ( k σ ≪ ω P), γ 0 / k σ → 0, and therefore, P / Q → 0. This implies that δ 1 → -1 2 , and therefore  \nγ ≈ γ 0 -ν c 2 (47)  \non large scales. In Fig. 3, this is evident from the linear decrease with ν c (at large ν c) of the damping time of the most weakly damped modes for the ω P / k σ = 4 . 5 and 2 . 5 cases, indicated by the light and medium dark blue dots respectively.  \nLet us now investigate the case of frequent collisions, i.e., when ν c ≳ k σ . As evident from Fig. 3, the modes remain oscillatory but become more strongly damped with increasing ν c, until a bifurcation occurs and they become non-oscillatory. At bifurcation, the complex frequencies of the modes become purely real and negative. The bifurcation occurs when ν c ∼ 2 ω P, i.e., when the mean free path, λ c = σ/ν c, is roughly half the Debye length, λ D = σ/ω P. Each bifurcation results in a less damped upper branch and a more damped lower branch. For the zeroth order modes (which analytically continue to the least damped Landau modes in the ν c / k σ → 0 limit), the damping timescale of the upper branch scales as t d ∼ ν c while that of the lower branch scales as t d ∼ 1 /ν c in the ν c ≫ 2 ω P limit. Hence, the zero-th order Landau branch mode is very weakly damped, albeit non-oscillatory for ν c ≫ 2 ω P. For the higher order modes, the damping timescales of both branches scale as t d ∼ 1 /ν c, i.e., these modes are strongly damped as well as non-oscillatory in the ν c ≫ 2 ω P limit.  \nThe above behaviour of the modes at large ν c / 2 ω P can be readily seen by analyzing the dispersion relation in the ν c / k σ → ∞ limit. We follow [12] to obtain the following form for the dispersion relation (for the zero-th order mode) in this limit:  \nD ( k , γ ) = 1 + ω P ν c ! 2 Z ∞ 0 d x x γ/ν c -1 (1 -x ) = 0 , (48)  \nwhich can be integrated to yield  \n1 + ω P ν c ! 2 " 1 γ/ν c -1 γ/ν c + 1 # = 0 . (49)  \nThis is a quadratic equation in γ/ν c and has the following solutions:  \nγ ν c = -1 2          1 ± s 1 -2 ω P ν c ! 2          . (50)  \nThe frequencies are complex conjugate, i.e., the zero-th order modes are damped but oscillating for ω P > ν c / 2. At ω P = ν c / 2, a bifurcation occurs and the frequencies become negative real, i.e., the modes become non-oscillatory and damped. For ν c ≫ 2 ω P, the upper branch mode has γ ≈ -ω 2 P /ν c and is weakly damped, while the lower branch one has γ ≈ -ν c and is strongly damped.  \nFrom Fig. 3 as well as the above analysis, we see that the discrete modes are more strongly damped than the Landau modes of the collisionless limit when ν c < 2 ω P. At ν c ≈ 2 ω P, each of these modes bifurcates into two non-oscillatory modes, one of which is less damped than the other. The less damped mode eventually becomes more weakly damped than the corresponding Landau mode at ν c ≫ 2 ω P. Therefore, collisions facilitate the collective, Landau-like damping of plasma perturbations for ν c < 2 ω P, i.e., when the mean free path, λ c = σ/ν c, exceeds half the Debye length, λ D = σ/ω P. On the other hand, for ν c ≫ 2 ω P, i.e., λ c ≪ λ D / 2, collisions oppose Landau damping, but make the modes non-oscillatory. Since the mean free path is smaller than the Debye length in this case, collisional scatterings destroy the coherence of the collective oscillations.', 'B. Lenard-Bernstein Branch': "The collision operator in the Boltzmann equation involves the highest velocity derivative and is therefore a singular perturbation to the collisionless Boltzmann equation. This observation foreshadows our finding that the discrete mode spectrum for the weakly collisional case is qualitatively di ff erent from the collisionless case. Therefore, in addition to the Landau branch modes discussed above, the LB operator introduces an additional set of modes, which we term LenardBernstein or LB modes. We indicate these using white dashed lines in Fig. 3. Unlike the Landau branch modes, the LB modes do not analytically continue to the Landau modes in the collisionless limit. Moreover, the LB modes are damped and non-oscillatory. As also found by [12], the frequencies of these modes are given by  \nγ ν c = -k σ ν c ! 2 -n , (51)  \nwhere n = 0 , 1 , 2 , ... . These modes are independent of the plasma frequency and therefore completely decoupled from the mean electrostatic field. These solely depend on ν c, and are thus pure manifestations of short range collisional scatterings in the plasma. The n , 0 modes are heavily damped, especially in large and small ν c / k σ limits, since γ ≈ -n ν c for ν c ≫ k σ and γ ≈ -( k σ ) 2 /ν c for ν c ≪ k σ . On the other hand, the n = 0 mode is heavily damped for small ν c / k σ but very weakly damped in the large ν c / k σ limit, since γ = -( k σ ) 2 /ν c throughout.  \nThe LB modes play a crucial role in plasma relaxation in the large ν c / k σ limit. The damping timescale of the LB modes  \nscales with ν c in the same way as that of the Landau branch modes in this limit. Also, in this limit, the damping time of the continuum response is the same as that of the n = 0 LB mode. Hence, it is the LB modes that dictate the relaxation of plasmas in the limit of frequent collisions, or equivalently, on scales larger than the mean free path.  \nVlasov simulations have demonstrated some of the phenomenology of wave damping described in this paper. [27] find in their simulation the survival of a non-oscillating and slowly decaying 'entropy mode' (also discussed by [21]) in the long-time limit. This bears qualitative resemblance to the LB modes [12] discussed in this paper. The long term, slow decay of the response in their simulation is very similar to that of the slowly decaying continuum and the n = 0 LB mode that we predict (although non-linear Landau damping e ff ects become important after approximately the bounce time or libration time for the wave-particle resonances). The pronounced damping of the response for ν c ∼ k σ , as predicted by our analysis, is seen by [28] in their Vlasov simulations, who note that the collision frequency suitable for preventing the 'numerical recurrence of initial conditions' (the echo; see section VII for a detailed discussion of how echoes damp) is a function of the wavenumber of the perturbation. We show that this wavenumber is nothing but ∼ ν c /σ , the inverse of the mean free path. Perturbations tend to damp away cubic exponentially within a mean free path, and undergo slow, exponential damping beyond.", 'VI. COMBINED LINEAR RESPONSE': 'In the above sections, we separately computed the continuum and discrete mode parts of the linear response. The combined linear response of a plasma is a temporal convolution of the two. This can be readily seen by studying the evolution of the potential perturbation, Φ 1 k ( t ). Note that Φ 1 k = -(2 π nee /ϵ 0)( G 1 κ (0) / k 2 ), where G 1 κ (0) is the inverse Laplace transform of ˜ G 1 κ (0), given in equation (16). To compute the inverse Laplace transform, we make use of the identity,  \nL -1 GLYPH<16> ˜ f ( γ )˜ g ( γ ) GLYPH<17> = 1 2 π i Z c + i ∞ c -i ∞ d γ ˜ f ( γ )˜ g ( γ ) exp GLYPH<2> γ t GLYPH<3> = Z t 0 d τ f ( τ ) g ( t -τ ) , (52)  \nwhere f and g are the inverse Laplace transforms of ˜ f and ˜ g respectively. For our problem, ˜ f and ˜ g correspond to N and 1 / D respectively, with N and D given by equations (17).  \nFor an initial Gaussian pulse, h κ ( u ), given by equation (21), we find that, for t ≳ 1 / k σ , the potential response is given by 5  \nFIG. 4. The combined temporal response of a weakly collisional plasma, including both the continuum and the discrete modes, in terms of | Φ 1 k | vs t (in units of 1 / k σ ). The initial perturbation in the DF is taken to be a Dirac delta function in v ( s = σ 0 /σ = 0). Di ff erent colors indicate the response for di ff erent values of ν c / k σ , and di ff erent panels correspond to di ff erent values of ω P / k σ as indicated. The initial phase of the decay, up to t ∼ 1 /ν c, is driven by the continuum. For ν c < k σ , the latter phase of the decay is driven by the discrete modes, while for ν c > k σ , the long term decay is governed by the continuum and the zero-th order LB mode, which is a non-oscillatory mode.  \n<!-- image -->  \nΦ 1 k ( t ) = -nee ϵ 0 h κ 0 k 2 X n exp GLYPH<2> γ nt GLYPH<3> D \' ( κ, γ n ) Z k σ t 0 d z exp " -GLYPH<18> γ n k σ + κ GLYPH<19> z + κ 2 2 GLYPH<16> 1 -e -z /κ GLYPH<17> nGLYPH<16> 1 -s 2 GLYPH<17> GLYPH<16> 1 -e -z /κ GLYPH<17> + 2 o -i κ u 0 GLYPH<16> 1 -e -z /κ GLYPH<17> # , (53)  \nwhere κ = k σ/ν c, s = σ 0 /σ , and u 0 = v 0 /σ with v 0 the drift velocity and σ 0 the velocity dispersion of the initial Gaussian (in v ) pulse. γ n are the frequencies of the discrete modes that follow the dispersion relation, D ( κ, γ n ) = 0, given in equation (39). D \' refers to the derivative of D with respect to γ n / k σ and is given by  \nD \' ( κ, γ n ) = -ω 2 P k σν c e ζ GLYPH<2> ζ -an γ ( an , ζ ) -Γ 2 ( an )( an -ζ ) 2 F 2( an , an ; an + 1 , an + 1; -ζ ) i , (54)  \nwhere an = γ n /ν c + κ 2 , ζ = κ 2 , and 2 F 2 refers to the confluent hypergeometric function of the second kind.  \nLet us investigate the temporal behaviour of Φ 1 k in some detail. It decays super-exponentially at small time and transitions to an exponential decay at large time. The nature of the initial super-exponential decay, however, depends on s = σ 0 /σ , i.e., on how cold the initial perturbation in the DF is relative to the equilibrium DF. Taylor expanding the argument of the exponential, we see that, at t ≲ 1 / k σ , Φ 1 k ( t ) decays as  \nΦ 1 k ( t ) ∼ exp " -k 2 σ 2 t 2 2 " s 2 + 2 3 -s 2 ! ν c t # # . (55)  \nFor a su ffi ciently cold perturbation with s ≲ √ 2 / 3, the response decays as ∼ exp h -GLYPH<16> 1 / 3 -s 2 / 2 GLYPH<17> k 2 σ 2 ν c t 3 i at small time. Therefore, a perturbation that is Dirac delta in v ( s = 0) initially decays as a cubic exponential, i.e., ∼ exp h -k 2 σ 2 ν c t 3 / 3 i . If, on the other hand, the initial perturbation is close to the equilibrium Maxwellian form ( s ≈ 1), the response decays as a Gaussian, i.e., ∼ exp h -k 2 σ 2 t 2 / 2 i . This key di ff erence in the initial decay owes to the fact that the potential (also density or temperature) response to a su ffi ciently cold perturbation decays at small time only due to collisional di ff usion, whereas that to a perturbation with a large spread in velocities decays due to phase-mixing as well. The Gaussian decay due to phase-mixing was noted by [20], [22] and many others, while the cubic exponential decay was noted by [13]. We show that these di ff erent decay scalings can be obtained from a general form (equations [53] and [55]). Over long time the response to an initial perturbation of any form decays exponentially at the rate of min h Re( γ 0) , ( k σ ) 2 /ν c i , under the combined action of collective e ff ects (discrete modes) and collisional damping, where γ 0 is the frequency of the least damped discrete mode.  \nFig. 4 plots Φ 1 k in units of neeh κ 0 /ϵ 0 k 2 as a function of time (in units of 1 / k σ ) for an initial Dirac delta pulse ( s = 0), for eight di ff erent values of ν c / k σ as indicated and for three di ff erent values of ω P / k σ in three di ff erent panels. The initial decay is cubic exponential, which is more pronounced for ν c < k σ . The long time decay is exponential, which manifests more prominently for ν c > k σ . The rate of this long time  \nexponential decay is given by min h Re( γ 0) , ( k σ ) 2 /ν c i . For ν c ≳ 2 ω P, i.e., when the mean free path is somewhat smaller than the Debye length, this decay rate is ∼ ( k σ ) 2 /ν c. The transition from the initial cubic exponential decay to the late time exponential decay is more pronounced on scales comparable to the mean free path ( ν c ∼ k σ ). On scales much smaller than the mean free path ( ν c ≪ k σ ), the small time cubic exponential decay is governed by the continuum while the long time exponential decay is governed by the least damped discrete mode (as evident from the temporal oscillations of the long time response). On scales larger than the mean free path ( ν c ≫ k σ ), the response undergoes a slow exponential decay [see also 17-22], guided by the weakly damped continuum and the least damped discrete mode (the zero-th order non-oscillatory LB mode), over a timescale ∼ ν c / ( k σ ) 2 . For ω P ≤ ν c / 2, the response becomes non-oscillatory, as evident from the absence of oscillations in the response for larger values of ν c. All in all, Fig. 4 illustrates the smooth changeover from the super-exponential decay of the response on scales smaller than the mean free path ( ν c < k σ ) and / or in the diffusive regime ( t < 1 /ν c), to the exponential decay on scales larger than the mean free path ( ν c > k σ ) and / or in the friction regime ( t > 1 /ν c).', 'VII. NON-LINEAR RESPONSE: PLASMA ECHO': 'So far, we have discussed how the super-exponential (cubic exponential) and exponential decays with two di ff erent timescales show up in the linear response of a weakly collisional plasmas as two asymptotic limits. Now we shall discuss the implications of these two timescales on the non-linear (second order) response.  \nThe key features of the continuum response in a weakly collisional plasma are: (i) it decays cubic exponentially with time for λ < λ c but exponentially for λ > λ c, and gets damped  \nthe most for λ ∼ λ c, and (ii) an externally imposed pulse slows down over time due to the collisional damping of particle velocities over a collision time, 1 /ν c. These behaviors not only impact the linear regime of plasma relaxation, but also have implications for the survival and persistence of coherent structures and non-linear phenomena. One such non-linear phenomenon is the plasma echo [13, 29-31], which attests to the time-reversibility of the Vlasov equation and arises from the continuum response. The fine-grained continuum response to a perturbation gives rise to a macroscopic disturbance, when a second perturbation is introduced into the system sometime after the first perturbation and interferes with the response to it.  \nLet us analyze the echo problem in some detail, and understand the impact of collisionality on this phenomenon. This requires us to evaluate the second order response of the plasma to an external perturbation. Since we are interested in the echo phenomenon, which is caused by the continuum response, we ignore the self-potentials Φ 1 and Φ 2 in this section 6 . We assume that the plasma is excited at t = 0 by an impulsive kick of sinusoidal electric field with wavenumber k 1. The plasma is impulsively kicked for a second time at t = T by another sinusoidal electric field with wavenumber k 2. The corresponding k Fourier mode of the potential is given by (see Appendix C for details)  \n˜ Φ k ( γ ) = t P [ A δ ( k -k 1) δ ( t ) + B ( k + k 2) δ ( t -T )] , (56)  \nwhere A and B are constants, and t P is the temporal pulsewidth that goes to zero, keeping At P and Bt P finite. The resulting second order response, g 2 κ = ( σ/ ne ) f 2 k is derived in detail in Appendix C, and is given by equations (C13) and (C12). The third term in equation (C13) is the relevant term for the plasma echo, and can be expressed as  \ng 2 κ ( u , t ) = -ab ( k 1 σ t P) ( k 2 σ t P) √ 2 π exp " -u 2 2 # δ ( k -( k 1 -k 2)) Θ ( t -T ) J κ, -κ 2 ( t , T ) , (57)  \nwhere a and b are the dimensionless strengths of the first and  \nsecond potential pulses respectively (see Appendix C for details), and J κ, -κ 2 is given by  \nFIG. 5. Plasma echo: Second order density perturbation of the k = ( k 1 -k 2) mode as a function of ∆ t / T , with ∆ t = t -T the time elapsed since the second pulse is introduced at T (see Appendix C and section VII for details). The first and second pulses have wavenumbers k 1 and k 2 respectively. Here we have adopted k 1 σ/ν c = 3 and k 2 σ/ν c = 2. Di ff erent lines indicate di ff erent values of ν c / k σ . The top and bottom panels correspond to k σ T = 1 and 5 respectively. The left panel indicates the density perturbation, r 2 κ , normalized by its maximum absolute value, and the right panel indicates r 2 κ/ ( ab ( k 1 σ t P) ( k 2 σ t P)), where a and b are the dimensionless strengths of the first and second potential pulses respectively, and t P is the temporal pulsewidth that goes to zero, keeping at P and bt P finite. The spikes in the density perturbation indicate the echos. These occur at ∆ t / T ≈ k 2 / ( k 1 -k 2) = 2 in the collisionless limit. The echo time approaches T with increasing ν c. Note that the echoes are detectable even for ν c / k σ ∼ 1, as long as k σ T is small, i.e., the initial response has not su ffi ciently decayed via phase-mixing.  \n<!-- image -->  \nσ  \n(b)  \nk  \nT  \n=  \n5  \nJ κ, -κ 2 ( t , T ) = " e -ν c(2 t -T ) + GLYPH<20> iue -ν c( t -T ) + κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 + κ 2 e -2 ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<21> GLYPH<20> iue -ν c t + κ e -ν c T GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 -κ 2 GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 1 -e -ν c(2 t -T ) GLYPH<17> GLYPH<21> # × exp         -κ 2         ν c ( t -T ) -GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> GLYPH<16> 3 -e -ν c( t -T ) GLYPH<17> 2                 exp         -κ 2 2         ν c T -GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 2 + e -2 ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17>GLYPH<17> 2                 × exp GLYPH<20> κκ 2 GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 GLYPH<21> × exp h -iu GLYPH<16> κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> -κ 2 e -ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17>GLYPH<17>i . (58)  \nEquations (57) and (58) above denote the second order response at the wavenumber, k = k 1 -k 2, which manifests the echo born out of the interference between the first and second pulses. Information about the initial growth of the echo is carried by the first, second and third factors of equation (58). The second and third factors encode information about the damping of the envelope of the echo that occurs over the two timescales discussed in section IV, τ d , 1 ∼ ( k σ ) -2 / 3 ν -1 / 3 c and τ d , 2 ∼ ν c / ( k σ ) 2 . The fourth factor is a phase-factor that scales as ∼ exp h -i v ν c GLYPH<16> ( k 1 -k 2) GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> -k 2 e -ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17>GLYPH<17>i (note that k = k 1 -k 2) and becomes independent of v when the argument is zero. This corresponds to a bunching of particles and an enhancement in the macroscopic density perturbation: the echo. It occurs at  \nt = t echo ≈ T + 1 ν c ln GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 1 + k 2 k 1 -k 2 GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> . (59)  \nThe approximation improves in accuracy for large ( k 1 -k 2) σ T . In the collisionless ( ν c T → 0) limit, the phase-factor scales as ∼ exp [ -iv (( k 1 -k 2) t -k 1 T )]. Accordingly, the echo occurs at t echo ≈ k 1 T / ( k 1 -k 2), which is also the ν c T → 0 limit of equation (59). This is the classic result for echo formation in a collisionless plasma [29]. In the strong collision ( ν c T → ∞ ) limit, t echo ≈ T + 1 ν c ln | k 1 / ( k 1 -k 2) | . Hence, the echo forms closer to T , the time of the second pulse (since the continuum response slows down due to collisional damping), and becomes increasingly di ffi cult to detect with increasing ν c.  \nThe echo involves the continuum response becoming independent of v , and therefore corresponds to a macroscopic disturbance in the plasma density. Keeping this in mind, the ( k 1 -k 2) mode of the second order density perturbation, r 2 κ = ρ 2 k /ρ 0, is computed by integrating equation (57) over u , and is given by equation (C15). The left panels of Fig. 5 plot r 2 κ , divided by its maximum absolute value, as a function of ∆ t / T = ( t -T ) / T . The right panels plot the absolute value of r 2 κ/ ( ab ( k 1 σ t P) ( k 2 σ t P)) as a function of ∆ t / T . The top and bottom panels correspond to ( k 1 -k 2) σ T = 1 and 5 respectively. We adopt k 1 σ/ν c = 3 and k 2 σ/ν c = 2. The spikes in the  \ndensity perturbation are the echoes. In the collisionless limit, the echoes occur at times close to ∆ t / T ≈ k 2 / ( k 1 -k 2) = 2. The time of the echo approaches T with increasing ν c, which is also evident from equation (59), and becomes indistinguishable from T in the ν c T →∞ limit. The right panels of Fig. 5 show that the response gets weaker with increasing ν c, decaying as ∼ k 1 k 2 σ 2 /ν 2 c in the ν c T → ∞ limit. The echoes become more pulse-like or localized in time for larger T . They decay super-exponentially for ν c < ( k 1 -k 2) σ but undergo a slow exponential decay for ν c > ( k 1 -k 2) σ . For large ( k 1 -k 2) σ T , the linear response prior to t = T , and consequently the echo, are heavily suppressed on scales comparable to the mean free path, i.e., for ν c ∼ ( k 1 -k 2) σ . This is because the damping timescale is shorter than the time of occurrence of the echo in this case. However, as long as ( k 1 -k 2) σ T is small, echoes are detectable in a weakly collisional plasma even when ν c / ( k 1 -k 2) σ ∼ 1, i.e., on scales comparable to the mean free path (see the top panels of Fig. 5). This is because, in this case, the initial response has not su ffi ciently phase-mixed and damped away due to collisions before the second pulse is introduced.', 'VIII. DISCUSSION AND SUMMARY': 'We have provided a general formalism for the relaxation of a weakly collisional plasma in the perturbative (linear and second order) regime. Following the seminal work of [12], we compute the linear response of a weakly collisional plasma using the perturbed Boltzmann-Poisson equations and demonstrate that the [13] continuum response and the discrete modes [9-12] are nothing but di ff erent limits of a general form of the response. We show that at late times and on scales larger than the Debye length (fluid regime), the continuum response undergoes a much slower decay than what [13] predicts. Unlike the [13] cubic exponential decay, this decay is a slow exponential just like the discrete modes, and is similar to the rates predicted by [9-11, 17-22]. This finding provides a rigorous justification for the establishment of LTE in the fluid limit. Finally, we compute the second order response of a weakly  \ncollisional plasma and demonstrate how the cubic exponential and exponential decays show up in the evolution of the plasma echo. Our findings about the collisional damping of the plasma echo have important implications for dissipation in turbulent plasmas [1].  \nWe describe a weakly collisional plasma using a perturbative analysis of the Boltzmann-Poisson equations. We model collisions using a collision operator of the Fokker-Planck type, called the Lenard-Bernstein (LB) operator [12], given by equation (7). The collision frequency, ν c, parametrizes the rate of collisions. The collision operator consists of two terms: a di ff usion term proportional to σ 2 ∂ 2 f /∂ v 2 and a damping / friction term proportional to ∂ ( v f ) /∂ v . Although one considers both terms to obtain the discrete mode spectrum in the standard theory [12], the analysis of the continuum spectrum and the echo has so far only been performed using the di ff usion term [13, 14]. This approximation simplifies the calculations but is inconsistent with the fluctuation-dissipation theorem. It is only valid in a limited spatio-temporal domain: at times smaller than the collision time and on scales smaller than the mean free path, i.e., only within the boundary layer. We include both the di ff usion and friction terms in the collision operator, and obtain the full response, which is a combination of the continuum response and the discrete modes. In fact, we find that the total response is a temporal convolution of the two. The continuum dominates the initial part of the response, while the long time response has contribution from both the discrete modes and the continuum.  \nCollisions a ff ect the discrete modes and the continuum in di ff erent ways. The discrete mode response exponentially damps away due to a combination of collective Landau damping and collisional e ff ects; the frequencies and damping rates of the discrete modes get shifted from the Landau values in presence of collisions. When ν c < 2 ω P, collisions enhance the damping rate of the discrete modes above the Landau damping rate of collisionless plasma. For ν c ≳ 2 ω P, collisions make the modes non-oscillatory and purely damped. For ν c ≫ 2 ω P, the damping rates of the modes become weaker than the corresponding Landau damping rates (this slow decay in the frequent collision / fluid regime has also been noted by [17-22]). Therefore, collisions facilitate the collective damping of the response when they are rare, but oppose damping when they are very frequent.  \nThe continuum part of the response evolves di ff erently with time than the discrete modes. At times smaller than the collision time, 1 /ν c, the continuum decays as ∼ exp h -GLYPH<0> t / t d , 1 GLYPH<1> 3 i over a damping timescale of t d , 1 ∼ ν -1 / 3 c ( k σ ) -2 / 3 , similar to the prediction by [13]. At t > 1 /ν c, however, it goes over to an exponential decay over a timescale of t d , 2 ∼ ν c / ( k σ ) 2 . The response damps away the most on scales comparable to the mean free path, i.e., for ν c ∼ k σ . On smaller scales the decay is super-exponential, as predicted by [13], but on scales larger than the mean free path, the decay is exponential and gets slower with increasing ν c. We supplement the FokkerPlanck analysis of the continuum response with a complementary analysis of stochastic collisional scattering using the Langevin equation, in the same spirit as the study of Brownian motion [24, 25]. We model the e ff ect of collisional di ff u-  \ntic acceleration term corresponding to white noise, and that of collisional damping by a friction term. The resulting velocity dispersion increases linearly with time and then saturates to the equilibrium value, σ 2 , after t ∼ 1 /ν c. Meanwhile, the dispersion in position increases as ∼ t 3 for t ≲ 1 /ν c and as ∼ t for t ≳ 1 /ν c. This is a classic feature of random walk or di ff usion in the velocity space, followed by damping or friction in the velocity space and di ff usion in the position space. The cubic exponential decay of the continuum response for t < 1 /ν c, followed by an exponential decay for t > 1 /ν c, is a consequence of collisional di ff usion followed by damping. We emphasize that the temporal nature of continuum decay obtained in this paper is not particularly sensitive to the Lenard-Bernstein form of the collision operator, but is a robust phenomenon (at least qualitatively) in any system governed by short range random / chaotic collisions.  \nThe general treatment of plasma relaxation discussed in this paper shows that the [13] scaling of continuum decay is only valid for small time and on scales smaller than the mean free path. This is expected since the [13] calculation neglects the damping term in the LB operator, and entirely relies on a boundary-layer analysis, as also demonstrated by [14]. Therefore, one should be cautious about the limited range of validity of the [13] result when addressing phenomena pertaining to collisional di ff usion in plasmas. One such phenomenon is the non-linear e ff ect of plasma echo, a macroscopic disturbance that occurs due to the interference of the continuum responses to two subsequent pulses of electric field with wavenumber k 1 and k 2 and has a wavenumber k = k 1 -k 2. On scales smaller than the Debye length, i.e., for small ν c ( ν c < k σ ), the echo damps cubic exponentially as [13] predicts, but on larger scales and / or large ν c ( ν c > k σ ), it damps far slower. However, with increasing ν c, the echo time approaches the time of initiation of the second pulse (since the continuum response slows down due to collisional damping), and becomes harder, albeit possible, to detect. Contrary to what [13] predict, we infer that the echo is detectable even for ν c as large as ∼ ( k 1 -k 2) σ , as long as the second pulse is introduced within a few phase-mixing timescales after the first.  \nThe essential features of plasma relaxation described in this paper can be understood by comparing the plasma frequency, ω P, the collision frequency, ν c, and the phase-mixing frequency, k σ , or alternatively by comparing three di ff erent length scales: the Debye length, λ D = σ/ω P, the mean free path, λ c = σ/ν c, and the perturbation scale, λ ∼ 1 / k . The selflimiting nature of collisions is manifest in the stagnated damping of both the discrete modes and the continuum beyond certain scales. The damping of the discrete modes occurs e ffi -ciently for ν c ≲ 2 ω P, but gets stagnated for ν c ≳ 2 ω P. In other words, collisions facilitate the collective Landau-like damping of the discrete modes when the mean free path, λ c = σ/ν c, is larger than half the Debye length, λ D = σ/ω P, and oppose collective damping when λ c ≲ λ D / 2. The discrete modes, however, become non-oscillating for λ c ≲ λ D / 2, since multiple collisional scatterings within the Debye scale destroy the coherence of the collective plasma oscillations. The collisional damping of the continuum response occurs rapidly for ν c ∼ k σ or λ ∼ λ c, but is arrested for ν c > k σ or λ > λ c.  \nThis transition of the damping rate at λ ∼ λ c owes to the fact that the collisional di ff usion of particle velocities is restricted within the mean free path. This means that an LTE is established only within the mean free path since the DF approaches the Maxwellian form locally (within λ c around a point), after a collision time, 1 /ν c. As a result, a local equation of state can be adopted and a fluid formalism implemented to study the evolution of plasma perturbations on scales above the mean free path. An ideal fluid description is, however, only valid as long as ν c ≪ 2 ω P, or the mean free path, λ c, is substantially larger than λ D / 2. When ν c ≳ 2 ω P, or λ c ≲ λ D / 2, the discrete modes are rendered non-oscillatory by collisions. Therefore, an ideal fluid description becomes questionable in this regime. To properly study the long time relaxation of plasmas in this regime, one should either use a collisional kinetic theory formalism as described in this paper [see also 17-22], or incorporate viscous dissipation in the corresponding fluid formalism.  \nWehave developed a general formalism for the perturbative (linear and second order) relaxation of weakly collisional plasmas in this paper. It would be interesting to investigate the e ff ect of weak collisions also on the non-linear relaxation of plasmas. The crucial role of collisional damping or friction in breaking the symmetry of wave-particle resonances (shifting and splitting of resonance lines) during the quasilinear relaxation of plasmas has been noted by [32]. Their analysis is, however, restricted to the domain where the collision frequency significantly exceeds the libration / bounce frequency of charged particles in the trapped region of the phase space. In future work, we intend to investigate quasilinear relaxation and non-linear relaxation / particle trapping when these two frequencies are comparable. In particular, we would like to address the following questions. How does weak collisionality a ff ect the trapping of ions and electrons in the non-linear Bernstein-Green-Kruskal or BGK [33] modes? Is the BGK mode stable in presence of collisions? Do collisions untrap the otherwise trapped particles, and if so, how fast? How do collisions impact the so-called sideband instability [34] of BGK modes? We envision developing a theory for phase space cascade and turbulence in a weakly collisional plasma, similar to that advocated by [8], but relaxing some of their  \nassumptions. We plan to extend the formalism presented in this paper to the study of the impact of weak collisions on self-gravitating systems, since such systems are governed by the Boltzmann-Poisson equations just like plasmas (with the crucial di ff erence that the right-hand-side of the Poisson equation carries the opposite sign relative to the plasma case). This formalism can be applied to the evolution of various astrophysical systems such as globular clusters and self-interacting dark matter halos, where collisions play an important role.', 'ACKNOWLEDGMENTS': "We wish to acknowledge Scott Tremaine, Chris Hamilton, Vin'ıcius Duarte, Martin Weinberg and Wrick Sengupta for insightful discussions and valuable suggestions. This research is supported by the National Science Foundation Award 2209471 and Princeton University.", 'Appendix A: Calculation of the linear response': 'The linear response of a weakly collisional plasma is computed by simultaneously solving the linearized Boltzmann and Poisson equations, which amounts to solving equation (12) for ˜ g 1 κ ( u ). Upon taking the Fourier transform in u , this becomes equation (14), which needs to be solved for ˜ G 1 κ ( w ), the Fourier transform of ˜ g 1 κ ( u ). The solution for ˜ G 1 κ ( w ) is given by equation (15), which involves ˜ G 1 κ (0). Therefore, substituting w = 0 in equation (15), we obtain the expression for ˜ G 1 κ (0) given in equations (16) and (17). Re-substituting the expression for ˜ G 1 κ (0) in equation (15) now yields the following expression for ˜ G 1 κ ( w ):  \n˜ G 1 κ ( w ) = ˜ G (1) 1 κ ( w ) + ˜ G (2) 1 κ ( w ) + ˜ G (3) 1 κ ( w ) , (A1)  \nwhere  \n˜ G (1) 1 κ ( w ) = exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 H κ ( w \' ) , ˜ G (2) 1 κ ( w ) = -2 π κ ω P ν c ! 2 N (2) ( κ, η ) D ( k , η ) × exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' G 0( w \' ) , ˜ G (3) 1 κ ( w ) = -κ A κ ( η ) D ( k , η ) × exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' G 0( w \' ) . (A2)  \nN (2) ( κ, η ) = ( -k ) -( η + κ 2 ) × Z 0 κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 H κ ( w \' ) , (A3)  \nHere,  \nand  \nD ( κ, η ) = 1 + 2 π κ ω P ν c ! 2 ( -κ ) -( η + κ 2 ) × Z 0 κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' G 0( w \' ) = 1 + ω P κν c ! 2 GLYPH<2> 1 -( a -ζ ) ζ -a exp GLYPH<2> ζ GLYPH<3> γ ( a , ζ ) GLYPH<3> , (A4)  \nwhere a = η + κ 2 , ζ = κ 2 , and γ ( a , ζ ) = R ζ 0 d z z a -1 exp [ -z ] is the lower incomplete gamma function. We have assumed that the unperturbed DF, g 0, is of Maxwellian form, i.e., g 0( u ) = exp h -u 2 / 2 i. √ 2 π , whose velocity Fourier transform is given by  \nG 0( w ) = 1 2 π Z ∞ -∞ d u exp [ -iwu ] g 0( u ) = 1 2 π exp h -w 2 / 2 i . (A5)  \nThe temporal response can be computed by taking the inverse Laplace transform and inverse Fourier transform of  \n˜ G 1 κ ( w ). The three di ff erent terms in the expression for ˜ G 1 κ ( w ) given in equation (A1) then respectively correspond to three di ff erent parts of the response: ˜ G (1) 1 κ corresponds to the continuum response, ˜ G (2) 1 κ corresponds to the discrete modes temporally convoluted with the continuum, and ˜ G (3) 1 κ corresponds to the direct response to the perturber or \'wake\', temporally convoluted with the continuum and the discrete modes.  \nComputation of the continuum part of the response amounts to taking the inverse Laplace transform and inverse Fourier transform (in u ) of ˜ G (1) 1 κ in equation (A2). The inverse Laplace transform of ( w -k ) -( η + κ 2 ) ( w \' -k ) ( η + κ 2 -1 ) (recall that γ = ν c η ) yields  \nL -1 GLYPH<20> ( w -κ ) -( η + κ 2 ) GLYPH<0> w \' -κ GLYPH<1> ( η + κ 2 -1 ) GLYPH<21> = 1 2 π i Z c + i ∞ c -i ∞ d η ( w -κ ) -( η + κ 2 ) GLYPH<0> w \' -κ GLYPH<1> ( η + κ 2 -1 ) exp GLYPH<2> ην c t GLYPH<3> = exp h -κ 2 ν c t i δ h w \' -GLYPH<16> we -ν c t + κ GLYPH<16> 1 -e -ν c t GLYPH<17>GLYPH<17>i . (A6)  \nNext, we integrate the inverse Laplace transform of ˜ G (1) 1 κ over w \' and take the inverse Fourier transform in u of the result to obtain the following expression for the continuum response:  \ng cont 1 κ ( u , t ) = Z ∞ -∞ d w exp [ iwu ] 1 2 π i Z c + i ∞ c -i ∞ d η exp GLYPH<2> ην c t GLYPH<3> ˜ G (1) 1 κ ( w , η ) (A7) = exp h -κ 2 ν c t i Z ∞ -∞ d w exp [ iwu ] exp " -w 2 2 -κ w # Z w κ d w \' exp " w \' 2 2 + κ w \' # δ h w \' -GLYPH<16> we -ν c t + κ GLYPH<16> 1 -e -ν c t GLYPH<17>GLYPH<17>i H κ ( w \' ) = r 2 1 -e -2 ν c t exp " -u 2 2 GLYPH<0> 1 -e -2 ν c t GLYPH<1> # exp GLYPH<20> -κ 2 GLYPH<18> ν c t -2 tanh ν c t 2 GLYPH<19>GLYPH<21> exp GLYPH<20> -i κ u tanh ν c t 2 GLYPH<21> F ( κ, ν c t ) , (A8)  \nwhere  \nF ( w , κ, ν c t ) = Z ∞ -∞ d z exp         -       z +        r 2 tanh ν c t 2 κ -i u p 2 GLYPH<0> 1 -e -2 ν c t GLYPH<1>               2         × H κ GLYPH<16> κ + ze -ν c t GLYPH<17> . (A9)  \nIn going from the second to third line of equation (A8), we have made the substitution, z = w -κ . The continuum response is evaluated in section IV A for di ff erent models of hk ( u ), whose velocity Fourier transform is Hk ( w ).', 'Appendix B: Calculation of velocity and position dispersions in the Langevin formalism': 'In section IV B, we gain physical insight into the workings of collisional damping and di ff usion by studying the evolution of the first and second moments of the continuum DF response, rather than that of the continuum response itself, which we study in section IV A. To achieve this, we solve the stochastic Langevin equation, given in equation (31), assuming a collisional friction proportional to -v and a white noise power spectrum for the stochastic acceleration due to collisions. Here we solve this equation to derive the velocity dispersion, σ 2 ( t ), and the position dispersion, σ 2 x ( t ).  \nThe Langevin equation (31) can be integrated once with respect to time to yield the following expression for the velocity of a particle:  \nv ( t ) = d x d t = v 0 e -ν c t + Z t 0 d t \' e -ν c( t -t \' ) ξ ( t \' ) . (B1)  \nThe velocity dispersion, σ ( t ), can now be computed as follows:  \nσ 2 ( t ) = GLYPH<28> GLYPH<16> v ( t ) -v 0 e -ν c t GLYPH<17> 2 GLYPH<29> = Z t 0 d t \' e -ν c( t -t \' ) Z t 0 d t \'\' e -ν c( t -t \'\' ) GLYPH<10> ξ ( t \' ) ξ ( t \'\' ) GLYPH<11> = D Z t 0 d t \' e -ν c( t -t \' ) Z t 0 d t \'\' e -ν c( t -t \'\' ) δ GLYPH<0> t \' -t \'\' GLYPH<1> = D Z t 0 d t \' e -2 ν c t \' = σ 2 GLYPH<16> 1 -e -2 ν c t GLYPH<17> . (B2)  \nIn going from the second to the third line, we have used the assumption of a white noise power spectrum for ξ ( t ), i.e., ⟨ ξ ( t ) ξ ( t \' ) ⟩ = D δ ( t -t \' ), and in going from the fourth to the fifth line, we have used the Einstein relation that D = ν c σ 2 .  \nEquation (B1) can be integrated one more time with respect to time to yield the position of a particle:  \nx ( t ) = v 0 ν c GLYPH<16> 1 -e -ν c t GLYPH<17> + Z t 0 d t \' Z t \' 0 d t \'\' e -ν c( t \' -t \'\' ) ξ ( t \'\' ) . (B3)  \nThe dispersion in the position, σ x ( t ), is computed as follows:  \nσ 2 x ( t ) = *" x ( t ) -v 0 ν c GLYPH<16> 1 -e -ν c t GLYPH<17> # 2 + = Z t 0 d t 1 Z t 1 0 d t 2 e -ν c( t 1 -t 2) Z t 0 d t 3 Z t 3 0 d t 4 e -ν c( t 3 -t 4) ⟨ ξ ( t 2) ξ ( t 4) ⟩ = D Z t 0 d t 1 Z t 1 0 d t 2 e -ν c( t 1 -t 2) Z t 0 d t 3 Z t 3 0 d t 4 e -ν c( t 3 -t 4) δ ( t 4 -t 2) , (B4)  \nwhere, in going from the second to the third line, we have used the assumption that ⟨ ξ ( t ) ξ ( t \' ) ⟩ = D δ ( t -t \' ). Noting that the integral over t 4 is non-zero only when t 2 ≤ t 3, and using the Einstein relation, D = ν c σ 2 , we perform the integral over t 4 to obtain  \nσ 2 x ( t ) = ν c σ 2 Z t 0 d t 1 e -ν c t 1 Z t 1 0 d t 2 e 2 ν c t 2 Z t t 2 d t 3 e -ν c t 3 . (B5)  \nPerforming the integral over t 3 yields  \nσ 2 x ( t ) = σ 2 Z t 0 d t 1 e -ν c t 1 Z t 1 0 d t 2 GLYPH<16> e ν c t 2 -e ν c(2 t 2 -t ) GLYPH<17> , (B6)  \nand performing the t 2 integral yields  \nσ 2 x ( t ) = σ 2 ν c Z t 0 d t 1 " 1 -3 2 e -ν c t 1 + 1 2 e -ν c( t + t 1) # . (B7)  \nFinally, performing the above t 1 integral yields the following expression for σ 2 x ( t ):  \nσ 2 x ( t ) = σ ν c ! 2 " ν c t -GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<0> 3 -e -ν c t GLYPH<1> 2 # . (B8)', 'Appendix C: Calculation of the second order response and the plasma echo': 'The linear response of a plasma consists of a free streaming continuum and damped discrete modes. When marginalized over the velocities, the continuum response phase-mixes away. The non-linear response is, however, fundamentally di ff erent in that the continuum response manifests a bunching of particles in the phase space at specific times, leading to macroscopic disturbances in the density. This phenomenon is known as the plasma echo, and can be readily seen at second order in perturbation.  \nThe evolution of the second order response, f 2, is governed by the following equation:  \n∂ f 2 ∂ t + v ∂ f 2 ∂ x = -e m ∂ f 1 ∂ v ∂ Φ P ∂ x + ν c ∂ ∂ v v f 2 + σ 2 ∂ f 2 ∂ v ! . (C1)  \nSince the echo phenomenon is caused by the continuum response, we ignore the self-potential Φ 1 of the plasma that leads to discrete modes. Non-dimensionalizing the above equation, then taking the Fourier transform in x and u = v /σ and the Laplace transform in t as detailed in section III, we ob-  \nin the following di ff erential equation for ˜ G 2 k ( w ), the FourierLaplace mode of f 2:  \n( w -κ ) ∂ ˜ G 2 κ ∂ w + ( w 2 + η ) ˜ G 2 κ ( w , η ) = -w Z d κ \' GLYPH<0> κ -κ \' GLYPH<1> Z d η \' A κ -κ \' ( η -η \' ) ˜ G 1 κ \' ( w , η \' ) . (C2)  \nHere, we have assumed that f 2( t = 0) = 0. The above di ff erential equation can be integrated using the method of integrating factor to yield the following solution for ˜ G 2 κ :  \n˜ G 2 κ ( w , η ) = -exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' × Z d κ \' GLYPH<0> κ -κ \' GLYPH<1> Z d η \' A κ -κ \' ( η -η \' ) ˜ G 1 κ \' ( w \' , η \' ) . (C3)  \nThe solution for ˜ G 1 κ is given by  \n˜ G 1 κ ( w , η ) = -κ A κ ( η ) exp " -w 2 2 -κ w # ( w -κ ) -( η + κ 2 ) Z w κ d w \' exp " w \' 2 2 + κ w \' # GLYPH<0> w \' -κ GLYPH<1> η + κ 2 -1 w \' G 0( w \' ) , (C4)  \nwhere we have assumed that H κ ( w ) = 0, i.e., f 1( t = 0) = 0.  \nFor simplicity, we assume hereon that the perturbing potential is impulsive in time. Let the perturbing potential be a combination of two impulses, one at t = 0 and the other at t = T . This implies that  \nΦ k ( t ) = t P h N (1) k δ ( t ) + N (2) k δ ( t -T ) i . (C5)  \nHere, N (1) k and N (2) k are the spatial Fourier transforms of the first and second pulses respectively, and t P is the temporal width of each pulse, which tends towards zero, keeping N (1) k t P and N (2) k t P finite. The Laplace transform of this, for t > T , is given by  \n˜ Φ k ( γ ) = t P h N (1) k + N (2) k e -γ T Θ ( t -T ) i . (C6)  \nIn non-dimensional form, we have  \nA κ ( η ) = ν c e m σ 2 ˜ Φ k ( γ ) = ν c t P h A κ + B κ e -ην c T Θ ( t -T ) i , (C7)  \nwhere A κ = ( e / m σ 2 ) N (1) k and B κ = ( e / m σ 2 ) N (2) k .  \nLet us first evaluate the linear response, G 1 κ ( w , t ), which is an essential ingredient of the second order response. Using the result for inverse Laplace transform from equation (A6), we get that  \nG 1 κ ( w , t ) = -κν c t P exp " -w 2 2 # × [ A κ L κ ( w , t ) + B κ L κ ( w , t -T ) Θ ( t -T )] , (C8)  \nwhere  \nL κ ( w , t ) = GLYPH<16> we -ν c t + κ GLYPH<16> 1 -e -ν c t GLYPH<17>GLYPH<17> × exp h -κ ( w -κ ) GLYPH<16> 1 -e -ν c t GLYPH<17> -κ 2 ν c t i . (C9)  \nNow, let us evaluate the second order response. We inverse Laplace transform equation (C3) to obtain  \nG 2 κ ( w , t ) = -ν c t P exp " -w 2 2 -κ w # Z w κ d w \' exp " w \' 2 2 + κ w \' # w \' Z d κ \' GLYPH<0> κ -κ \' GLYPH<1> × h A κ -κ \' exp h -κ 2 ν c t i δ h w \' -GLYPH<16> we -ν c t + κ GLYPH<16> 1 -e -ν c t GLYPH<17>GLYPH<17>i G 1 κ \' ( w \' , 0) + Θ ( t -T ) B κ -κ \' exp h -κ 2 ν c ( t -T ) i δ h w \' -GLYPH<16> we -ν c( t -T ) + κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17>GLYPH<17>i G 1 κ \' ( w \' , T ) i . (C10)  \nPerforming the w \' integral and the velocity Fourier transform  \nof the above, we obtain the second order response as follows:  \ng 2 κ ( u , t ) = 1 √ 2 π exp " -u 2 2 # ( ν c t P) 2 Z d κ \' κ \' GLYPH<0> κ -κ \' GLYPH<1> " A κ -κ \' A κ \' I κ ( t ) + Θ ( t -T ) B κ -κ \' GLYPH<16> B κ \' I κ ( t -T ) + A κ \' J κ,κ \' ( t , T ) GLYPH<17> # , (C11)  \nwhere  \nI κ ( t ) = " e -2 ν c t + GLYPH<18> iue -ν c t + κ GLYPH<16> 1 -e -ν c t GLYPH<17> 2 GLYPH<19> 2 # exp " -κ 2 ν c t -GLYPH<0> 1 -e -ν c t GLYPH<1> GLYPH<0> 3 -e -ν c t GLYPH<1> 2 !# exp h -i κ u GLYPH<16> 1 -e -ν c t GLYPH<17>i ,  \nJ κ,κ \' ( t , T ) = h e -ν c(2 t -T ) + GLYPH<20> iue -ν c( t -T ) + κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 -κ \' e -2 ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<21> GLYPH<20> iue -ν c t + κ e -ν c T GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 + κ \' GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 1 -e -ν c(2 t -T ) GLYPH<17> GLYPH<21>GLYPH<21> × exp         -κ 2         ν c ( t -T ) -GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> GLYPH<16> 3 -e -ν c( t -T ) GLYPH<17> 2                 exp         -κ \' 2         ν c T -GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 2 + e -2 ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17>GLYPH<17> 2                 × exp GLYPH<20> -κκ \' GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> 2 GLYPH<21> × exp h -iu GLYPH<16> κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> + κ \' e -ν c( t -T ) GLYPH<16> 1 -e -ν c T GLYPH<17>GLYPH<17>i . (C12)  \nNote that I κ ( t ) = J κ,κ \' = 0( t , T = 0).  \nThe above response is a convolution over k , and can be simplified if we adopt a perturber potential that is sinusoidal in x . Assuming that Φ P ∼ exp [ ik 1 x ] at t = 0 and Φ P ∼ exp [ -ik 2 x ] at time T , we have that A κ = a δ ( k -k 1) and B κ = b δ ( k + k 2). Substituting these in equation (C11) and performing the integration over κ \' yields the following form for the second order response:  \ng 2 κ ( u , t ) = 1 √ 2 π exp " -u 2 2 # × " a 2 ( k 1 σ t P) 2 δ ( k -2 k 1) I κ ( t ) + b 2 ( k 2 σ t P) 2 δ ( k + 2 k 2) Θ ( t -T ) I κ ( t -T ) -ab ( k 1 σ t P) ( k 2 σ t P) δ ( k -( k 1 -k 2)) Θ ( t -T ) J κ, -κ 2 ( t , T ) # . (C13)  \nThe first and second terms denote parts of the response with wavenumbers 2 k 1 and -2 k 2 respectively. The third term is a cross-term between the k 1 and k 2 modes and manifests an echo of wavenumber k 1 -k 2 that occurs slightly after t = T . The  \nphase factor in this cross-term becomes independent of v , or in other words, phase space bunching occurs at  \nt = t echo ≈ T + 1 ν c ln GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 1 + k 2 k 1 -k 2 GLYPH<16> 1 -e -ν c T GLYPH<17> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> . (C14)  \nIn the collisionless ( κ = k σ/ν c → ∞ , κ 2 = k 2 σ/ν c → ∞ ) limit, the phase factor scales as ∼ exp [ -iv (( k 1 -k 2) t -k 1 T )].  \nAccordingly, the ν c T → 0 limit of equation C14 implies t echo ≈ [1 + k 2 / ( k 1 -k 2)] T = k 1 T / ( k 1 -k 2). In the strong collision or ν c T → ∞ limit, t echo -T ≈ 1 ν c ln | k 1 / ( k 1 -k 2) | . Hence, it becomes increasingly di ffi cult to observe the echo with increasing ν c.  \nThe echo is a manifestation of the bunching of particles in the phase space. This implies that the macroscopic density perturbation of the ( k 1 -k 2) mode (obtained by integrating the cross term in g 2 κ from equation [C13] over u ), given by  \nr 2 κ = ρ 2 k ρ 0 = -ab ( k 1 σ t P) ( k 2 σ t P) δ ( k -( k 1 -k 2)) Θ ( t -T ) × κ GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> h κ e -ν c T GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> -κ 2 GLYPH<16> 1 -e -ν c T GLYPH<17>i × exp h -κ 2 h ν c ( t -T ) -GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17>i -κ 2 2 h ν c T -GLYPH<16> 1 -e -ν c T GLYPH<17>i + κκ 2 GLYPH<16> 1 -e -ν c( t -T ) GLYPH<17> GLYPH<16> 1 -e -ν c T GLYPH<17>i (C15)  \npeaks around t = t echo (see Fig. 5).  \n- [1] A. A. Schekochihin, J. T. Parker, E. G. Highcock, P. J. Dellar, W. Dorland, and G. W. Hammett, Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence, Journal of Plasma Physics 82 , 905820212 (2016).\n- [2] R. Bruno and V. Carbone, The Solar Wind as a Turbulence Laboratory, Living Reviews in Solar Physics 2 , 4 (2005).\n- [3] E. Marsch, Kinetic Physics of the Solar Corona and Solar Wind, Living Reviews in Solar Physics 3 , 1 (2006).\n- [4] G. G. Howes, A prospectus on kinetic heliophysics, Physics of Plasmas 24 , 055907 (2017), arXiv:1705.07840 [physics.plasmph].\n- [5] J. L. Burch, T. E. Moore, R. B. Torbert, and B. L. Giles, Magnetospheric Multiscale Overview and Science Objectives, Space Science Reviews 199 , 5 (2016).\n- [6] I. H. Hutchinson, Electron holes in phase space: What they are and why they matter, Physics of Plasmas 24 , 055601 (2017).\n- [7] A. Chael, M. Rowan, R. Narayan, M. Johnson, and L. Sironi, The role of electron heating physics in images and variability of the Galactic Centre black hole Sagittarius A*, MNRAS 478 , 5209 (2018), arXiv:1804.06416 [astro-ph.HE].\n- [8] M. L. Nastac, R. J. Ewart, W. Sengupta, A. A. Schekochihin, M. Barnes, and W. D. Dorland, Phase-space entropy cascade and irreversibility of stochastic heating in nearly collisionless plasma turbulence, arXiv e-prints , arXiv:2310.18211 (2023), arXiv:2310.18211 [physics.plasm-ph].\n- [9] C. S. Ng, A. Bhattacharjee, and F. Ski ff , Kinetic Eigenmodes and Discrete Spectrum of Langmuir Oscillations in a Weakly Collisional Plasma, in APS Division of Plasma Physics Meeting Abstracts , APS Meeting Abstracts, Vol. 41 (1999) p. CP1.36.\n- [10] C. S. Ng, A. Bhattacharjee, and F. Ski ff , Complete Spectrum of Kinetic Eigenmodes for Plasma Oscillations in a Weakly Collisional Plasma, Phys. Rev. Lett. 92 , 065002 (2004).\n- [11] C. S. Ng and A. Bhattacharjee, Landau Modes are Eigenmodes of Stellar Systems in the Limit of Zero Collisions, Astrophys. J. 923 , 271 (2021), arXiv:2109.07806 [astro-ph.GA].\n- [12] A. Lenard and I. B. Bernstein, Plasma Oscillations with Di ff u-  \nsion in Velocity Space, Physical Review 112 , 1456 (1958).  \n- [13] C. H. Su and C. Oberman, Collisional Damping of a Plasma Echo, Phys. Rev. Lett. 20 , 427 (1968).\n- [14] R. W. Short and A. Simon, Damping of perturbations in weakly collisional plasmas, Physics of Plasmas 9 , 3245 (2002).\n- [15] C. Black, K. Germaschewski, A. Bhattacharjee, and C. S. Ng, Discrete kinetic eigenmode spectra of electron plasma oscillations in weakly collisional plasma: A numerical study, Physics of Plasmas 20 , 012125 (2013), https: // pubs.aip.org / aip / pop / articlepdf / doi / 10.1063 / 1.4789882 / 13427281 / 012125 1 online.pdf.\n- [16] S. Chandrasekhar, Dynamical Friction. I. General Considerations: the Coe ffi cient of Dynamical Friction., Astrophys. J. 97 , 255 (1943).\n- [17] A. R. Bell, Electron energy transport in ion waves and its relevance to laser-produced plasmas, The Physics of Fluids 26 , 279 (1983), https: // pubs.aip.org / aip / pfl / articlepdf / 26 / 1 / 279 / 12262372 / 279 1 online.pdf.\n- [18] E. M. Epperlein, R. W. Short, and A. Simon, Damping of ionacoustic waves in the presence of electron-ion collisions, Phys. Rev. Lett. 69 , 1765 (1992).\n- [19] V. Y. Bychenkov, J. Myatt, W. Rozmus, and V. T. Tikhonchuk, Kinetic theory of ion acoustic waves in a plasma with collisional electrons, Phys. Rev. E 52 , 6759 (1995).\n- [20] A. V. Brantov, V. Y. Bychenkov, and W. Rozmus, Relaxation of a thermal perturbation in a collisional plasma, Plasma Physics Reports 32 , 337 (2006).\n- [21] A. V. Brantov, V. Y. Bychenkov, and W. Rozmus, Electrostatic response of a two-component plasma with coulomb collisions, Phys. Rev. Lett. 108 , 205001 (2012).\n- [22] W. Rozmus, A. Brantov, C. Fortmann-Grote, V. Y. Bychenkov, and S. Glenzer, Electrostatic fluctuations in collisional plasmas, Phys. Rev. E 96 , 043207 (2017).\n- [23] P.-H. Chavanis, Initial value problem for the linearized mean field Kramers equation with long-range interactions, European Physical Journal Plus 128 , 106 (2013), arXiv:1306.1225 [cond-  \nmat.stat-mech].  \n- [24] A. Einstein, Uber die von der molekularkinetischen Theorie der Warme geforderte Bewegung von in ruhenden Flussigkeiten suspendierten Teilchen, Annalen der Physik 322 , 549 (1905).\n- [25] D. S. Lemons and A. Gythiel, Paul Langevin\'s 1908 paper \'On the Theory of Brownian Motion\' [\'Sur la th\'eorie du mouvement brownien,\' C. R. Acad. Sci. (Paris) 146, 530-533 (1908)], American Journal of Physics 65 , 1079 (1997).\n- [26] L. D. Landau, On the vibrations of the electronic plasma, J. Phys. (USSR) 10 , 25 (1946).\n- [27] J. W. Banks, S. Brunner, R. L. Berger, and T. M. Tran, Vlasov simulations of electron-ion collision e ff ects on damping of electron plasma waves, Physics of Plasmas 23 , 032108 (2016), https: // pubs.aip.org / aip / pop / articlepdf / doi / 10.1063 / 1.4943194 / 14894612 / 032108 1 online.pdf.\n- [28] O. Pezzi, E. Camporeale, and F. Valentini, Collisional effects on the numerical recurrence in Vlasov-Poisson simula-  \n- tions, Physics of Plasmas 23 , 022103 (2016), arXiv:1601.05240 [physics.plasm-ph].\n- [29] R. W. Gould, T. M. O\'Neil, and J. H. Malmberg, Plasma Wave Echo, Phys. Rev. Lett. 19 , 219 (1967).\n- [30] V. N. Pavlenko, Echo phenomena in plasmas, Soviet Physics Uspekhi 26 , 931 (1983).\n- [31] L. K. Spentzouris, J.-F. Ostiguy, and P. L. Colestock, Direct measurement of di ff usion rates in high energy synchrotrons using longitudinal beam echoes, Phys. Rev. Lett. 76 , 620 (1996).\n- [32] V. N. Duarte, J. B. Lestz, N. N. Gorelenkov, and R. B. White, Shifting and splitting of resonance lines due to dynamical friction in plasmas, Phys. Rev. Lett. 130 , 105101 (2023).\n- [33] I. B. Bernstein, J. M. Greene, and M. D. Kruskal, Exact Nonlinear Plasma Oscillations, Physical Review 108 , 546 (1957).\n- [34] C. B. Wharton, J. H. Malmberg, and T. M. O\'Neil, Nonlinear E ff ects of Large-Amplitude Plasma Waves, Physics of Fluids 11 , 1761 (1968).'}