Out-of-equilibrium mean-field dynamics of a model for wave-particle interaction

The out-of-equilibrium mean-field dynamics of a model for wave-particle interaction is investigated. Such a model can be regarded as a general formulation for all those applications where the complex interplay between particles and fields is known to be central, e.g., electrostatic instabilities in plasma physics, particle acceleration and free-electron lasers (FELs). The latter case is here assumed as a paradigmatic example. A transition separating different macroscopic regimes is numerically identified and interpreted by making use of the so-called violent relaxation theory. In the context of free-electron lasers, such a theory is showed to be effective in predicting the saturated regime for energies below the transition. The transition is explained as a dynamical switch between two metastable regimes, and is related to the properties of a stationary point of an entropic functional.


I. INTRODUCTION
Mean-field models have been widely studied as paradigmatic representatives of the important class of systems subject to long-range coupling. In the simplest scenario, N particles are made to interact in one dimension, subject to a varying field which is self-consistently sensitive to the individual trajectories. A global network of connections is hence driving the dynamics of every constituting element, as it certainly happens for more realistic settings where e.g. gravity or unscreened Coulomb interactions are at play.
The dynamics of mean-field models displays intriguing features. Particles may be trapped in intermediate (out-of-equilibrium) states, whose duration diverges with the number of constitutive elements, and which substantially differ from the corresponding thermodynamic equilibrium configuration. These metastable states are often termed in the literature Quasi-Stationary States, hereafter QSS, and bear an extraordinary conceptual importance as they potentially corresponds to the solely experimentally accessible regimes, for a wide range of applications ranging from celestial mechanics to plasma physics.
Interestingly, the evolution of the QSS is intimately governed by the discreteness of the medium being investigated. More specifically, QSS are believed to correspond to stationary stable solution of a Vlasov model, invoked as the continuous analogue of the discrete N -particles dynamics.
Within this context, the Hamiltonian Mean Field model (HMF) [1] has often been referred to as the bench-mark model for elaborating onto the QSS emergence. This is a one-dimensional Hamiltonian describing the evolution of N rotors coupled via a mean-field cosinuslike potential. The QSS in the HMF setting have been explained by resorting to a maximum entropy principle, pioneered by Lynden-Bell in astrophysical context, and fully justified from first principles [2]. Here, the supposedly relevant Vlasov picture enters the description as a Fermionic contribution to an entropy functional. Besides, out-of-equilibrium phase transitions are predicted to occur, separating between distinct macroscopic QSSs.
The Lynden-Bell protocol, also termed violent relaxation theory, was also argued to apply to other mean-field models and, indeed, it showed effective in predicting the saturate intensity of a free-electron laser (FEL). FELs are lasing devices consisting of a relativistic beam of charged particles, interacting with a co-propagating electromagnetic wave. The interaction is assisted by the static and periodic magnetic field generated by an undulator. FELs admit a mean-field description in term of the so-called Colson-Bonifacio model [3,4], which captures the essence of the collective wave-particle dynamics. However, no detailed study has been carried out for the FEL case aiming at unravelling the possible existence of out-of-equilibrium transition of the type mentioned above. Are these transitions ubiquitous in mean field dynamics and, in this case, can we provide a consistent intrepretative framework for their emergence? This paper is dedicated to answering such questions, and makes reference to the specific FEL setting. We also stress that the Colson-Bonifacio model of FEL dynamics can be regarded as a general formula-tion for all those applications where the complex interplay between particles and fields is well known to be central, e.g. electrostatic instabilities in plasma physics [5]. As a closing remark, we notice that, on the practical implication side, by disposing of reliable predictive tools on the system evolution, one can aim at guiding the system towards different experimental regimes.

II. THE FEL MODEL
The Colson-Bonifacio model for the FEL dynamics describes the coupled evolution of the electrons with a copropagating wave. The equations read: where θ j stands for the particle phase with respect to that of the optical wave, p j being its conjugate normalized momentum. The complex quantity A = A x + iA y represents the transverse field and N the number of electrons composing the electron bunch [27].
Herez labels the longitudinal position along the undulator, and it effectively plays the role of time. The intensity of the laser field is I = A 2 x + A 2 y . As it can be seen from the last of Eqs. (1), the bunching term, b = 1 N j e −iθj , is the source of wave amplification. The bunching quantifies the degree of localization of the electrons in the generalized space of their associated phases. The above discrete system of equations admits a Hamiltonian formulation to which we shall make reference as to the N-body model. In the N → ∞ limit, system (1) converges to the following Vlasov-wave set of equations [6]: Eqs.
(2) can be simulated numerically, thus allowing us to monitor the evolution of the phase space distribution function f (θ, p) along thez axis. In our implementation, we adopt the semi-Lagrangian method [7], associated with a cubic spline interpolation [8]. Results of the numerical integration are also checked versus N -body simulations and shown to return a perfect matching on relatively short time scale, for large enough values of N . On longer times, finite-N corrections do matter. The discrete system is in turn sensitive to intrinsic granularity effects, stemming from the intimate finiteness of the simulated medium, and progressively migrate from the Vlasov state towards the deputed equilibrium configuration. When increasing its size, the system spends progressively more time in the Vlasov-like, outof-equilibrium regime. Formally, in the N → ∞ limit, it never reaches equilibrium, being permanently trapped in the QSS.
As previously anticipated, our study is hence ultimately concerned with the emergence of QSSs, in a context where particles and waves evolve self-consistently. We shall be in particular interested in elucidating the occurrence of out-of-equilibrium phase transitions via dedicated numerical simulations, and substantiate our claims analytically. In doing so, we will virtually extend the conclusion of [9] to a broad spectrum of potentially relevant applications, beyond the specific case under inspection. Among others, it is again worth mentioning plasma physics: A formulation equivalent to model (2) is in fact often invoked, when studying the collective effects of beam-plasma dynamics [5].

III. ON THE INITIAL CONDITIONS AND THEIR SUBSEQUENT DYNAMICAL EVOLUTION
Let us turn to discussing our results, as obtained via numerical integration of (2). In order to make contact with the investigations reported in [9], we shall employ in the following a two-dimensional water-bag initial condition in phase space, which can be seen as a rough approximation of a smooth Gaussian profile. A (rectangular) water-bag is formally parametrized by two quantities, namely the semi-width of the spanned interval in phase, ∆θ, and its homologue in the momentum direction, ∆p. The corresponding expression for f can be cast in the form (see also The initial conditions can be also characterized by defining where b 0 is the initial bunching, and ǫ the initial average kinetic energy per particle. Notice that we thus access all possible values of the bunching b 0 ∈ [0, 1] by properly tuning ∆θ, and all positive energies ǫ by varying ∆p. We here limit our discussion to the case of vanishing initial optical field, I 0 ≃ 0, the relevant parameters' space being therefore solely bound to the plan (b 0 , ǫ). The initial condition here selected is hence solely controlled by these two parameters. In other words, we are specializing on a given bidimensional subset, and deliberately ignore the third, in principle available, direction of the reference parameter space. Quantifying the role of such an additional degree of freedom ultimately amounts to investigate the so-called seeded configuration [12,13] and will be the subject of future work. Let us start by discussing the simplest scenario, where the initial beam of particles is uniformly distributed over [−π; π]. From a physical point of view, this amounts to specialize to the case of Self-Amplified Spontaneous Emission (SASE, [14,15]), where b 0 = 0 and no seed is applied externally. Such a choice was also considered by Barré et al. [6] and Curbis et al. [10], where the dependence of the system evolution on the energy was numerically monitored, within the N -body discrete viewpoint. Interestingly, b 0 = 0 is a stationary solution of the Vlasov system: a local perturbative calculation can hence be straightforwardly implemented so to investigate its inherent stability [11]; the calculations are detailed in appendix A. For ǫ < 0.315, an instability occurs: both the wave intensity and the bunching factor rapidly grow, before relaxing towards an oscillating plateau. The average value of I reached in the oscillating regime is called the saturated intensityĪ. This behaviour is displayed in Fig. 1, where the simulations with ǫ > 0.315 do not show an amplification of I. This is a well-known property, indeed correctly reproduced by our numerical simulations, and which first signals the existence of phase transitions, of the type depicted in [9]. To further corroborate our guess on the b 0 = 0 behaviour, we turn to measuring the saturated intensityĪ as function of the energy ǫ, wherē I stands for the mean of I after saturation. It is here computed during four oscillations of I.
As shown in Fig. 2,Ī rapidly shrinks, when increasing the energy ǫ, until a critical value is reached where a sudden transition toĪ ≃ 0 is observed, bearing the characteristic of a first order phase transition. This is a further point of contact with the analysis carried on in [9] for the HMF toy model.
Motivated by these findings, and to push the analogy with the HMF setting, we consider bunched initial distributions. From a physical point of view, this choice is relevant to the case of FELs working in the so-called harmonic generation regime [12,13]. Particles' positions are here initially assigned so to uniformly span a limited portion of the allowed support, symmetric with respect to the origin, controlling the associated bunching via Eq. (4). The inhomogeneous (b 0 = 0) distribution in phase space is by nature non-stationary. Vlasov dynamics can however smooth it to a homogeneous, (b = 0, I = 0) possibly non water-bag, distribution or evolve to a bunched situation. The saturated mean-field average intensityĪ vs. the energy parameter ǫ is represented in Fig. 2, showing the newly collected data for different values of b 0 > 0 to the reference profile relative to b 0 = 0. In all cases the intensity is shown to decrease, as the energy increases. Importantly, for small values of b 0 , an abrupt transition is observed, which can be naively interpreted as of the firstorder type. For larger values of b 0 , the observed transition becomes smoother, such as for a second-order one. A substantially identical scenario holds for the bunching, which evolve towards an asymptotic plateaub, also sensitive to the ǫ and b 0 parameters, see Fig. 3. This scenario points towards a unifying picture on the emergence of out-of-equilibrium phase transitions within the considered class of mean-field Hamiltonian model. As previously anticipated, the Lynden-Bell theory of violent relaxation was successfully applied to the HMF problem, allowing one to gain a comprehensive understanding on the out-of-equilibrium phase transition issue, including a rather accurate characterization of the associated transition order. In the following section we set down to apply the Lynden-Bell argument to the present case, benchmarking the theory to numerical experiments. Before ending this section, we briefly discuss the phase-space structures resulting from the Vlasov-based simulations. Two phase-space portraits are enclosed in Fig. 4, and refer to different values of the energy ǫ, respectively below (upper panel) and above (lower panel) the critical transition energy relative to the selected (fixed) b 0 amount. When the system evolves towards a state atĪ = 0, then f (θ, p) shows a large resonance. At variance, in the opposite regime, a holeresonance dipole structure is observed, see also [16]. This observation seems to suggest that the out-of-equilibrium phase transition materializes via a bifurcation of invariant structures, an observation that was recently made for the HMF model [17]. Results of further investigations on this specific topic will be presented in a separate contribution [18].

IV. ON THE VIOLENT RELAXATION THEORY
In his work on self-gravitating systems, Lynden-Bell suggested [2] that the collisionless dynamics governed by the Vlasov equation tends to maximize a Fermionic entropy. The latter is obtained from the classical definition, where the counting of the microscopic states, compatible with a given macroscopic configuration, results from a combinatorial calculation, and is sensitive to the underlying Vlasov dynamics. The method was successfully employed in the study of the HMF model [9,19,20] and also applied to predict the quasi-stationary amplitude of the FEL wave [6,10]. In these works, however, the analysis just focused on the unstable regime (Ī = 0): no attempt was in fact made to reconcile it, with the high energy homogeneous state, via the phenomenon of out-of-equilibrium phase transitions. More recently, Yamaguchi [21] used the Lynden-Bell approach to predict the core of the gravitational sheet model, demonstrating its adequacy within the field for which it was originally conceived. In the following we shall review the main steps of the derivation of the violent relaxation theory, applied to the FEL setting. Starting from a water-bag, the entropy to be maximized [2] can be cast in the form [22] s(f ) = − dp dθ f (3) andf is the coarse grained distribution function. Following Barré et al. [6,23], maximizing the functional (5), results in the following set of equations where ζ = exp (−2Aβ sin θ), F 0 (y) = Here β and x are (rescaled) Lagrange multipliers and ultimately stem from the conservation of mass, momentum and energy. A, β and x are calculated by solving Eqs. (6) numerically via a Newton-Raphson method. The resulting (real) value of A is expected to return an estimate of the laser intensity at (Vlasov) saturation,Ī = A 2 , while f (θ, p) is: f (θ, p) = f 0 1 1 + x e β(p 2 /2+2A sin θ+A 2 p+A 4 /2) . (7) For A = 0 (namelyĪ = 0) the optimization problem (6) reduces to: An x value exists which solves the above equation for any choice of ∆θ. The homogeneous state is a stationary solution of the Lynden-Bell entropy, and so a potentially attractive state of the Vlasov dynamics. Additional inhomogeneous solutions (A = 0, or, equivalently,Ī = 0) might however emerge from investigating the full system (6). The homogeneous and inhomogeneous solutions will be referred to as to LB0 and LBA, respectively (see Fig.  5). The forthcoming discussion will focus on how to discriminate between the two, and eventually predict the asymptotic fate of the system.

V. INTERPRETING THE OUT-OF-EQUILIBRIUM TRANSITION AS A DYNAMICAL SWITCH BETWEEN LBA AND LB0
Let us focus first on the LBA solution. In Fig. 6, we report the value ofĪ as it follows from Eqs. (6), for different choices of the energy and initial bunching. The predicted intensity is shown to decrease, when the energy gets larger, but no transition is observed, in contradiction with the results of our numerical simulations. As previously stressed, the homogeneous LB0 state is also solution of the optimization problem (6) and could in principle prevail over the former. To shed light on this issue, we calculated the entropy values S A and S 0 , associated to LBA and LB0, respectively. Results of the computations are shown in Fig. 7, where the dependence on the energy ǫ is monitored for various choices of b 0 . Surprisingly, and at odd with what happens for the HMF model [9], S A is always larger than S 0 . The two curves do not cross each other and the LBA configuration is entropically favoured. Let us note that for higher values of ǫ,Ī decreases and the two solutions get close to each other, also from the point of view of the entropy, without crossing. The observed transition could possibly stem from a purely dynamical mechanism, and this would justify the discrepancy between the simulation output and the statistical prediction. More specifically, we here argue that, depending on the selected initial conditions, the system explores a local basin of attraction and struggles to find its way to the deputed, global maximum of the entropy. To clarify this point, we focus on a single numerical simulation, assuming the system to be initialized in a LB0 state. The dynamics can progressively take the system towards the LBA configuration, respecting the maximization of the Lynden-Bell entropy. The opposite is not possible, and a simulation started in the LBA state will certainly not evolve to the LB0. However, dynamical effects might be also at play and interfere with ideal situation here schematized, by virtually blocking the system in the neighborhood of an initially assigned LB0 conformation. Is this the correct scenario? And how can one explain the observed transition when starting from waterbag initial conditions? These issues are addressed in the following, where the stability of LB0 and LBA is investigated via direct Vlasov simulations.
In Fig. 8 the dynamical evolution of the intensity I is depicted, for three different classes of initial conditions, relative to the same choice of ǫ and b 0 . The LBA is indeed stable, no deviation from the initial configuration being observed as an effect of the Vlasov dynamics. Conversely, the LB0 condition proves unstable, and the intensity converges towards an oscillating plateau. Interestingly, the LB0 and water-bag (WB) evolutions are qualitatively similar, and, moreover, display the same average asymptotic value for the intensity I. Even more important, the asymptotic value corresponds to the LBA (maximum en-  To further elucidate the analogies between LB0 and WB, and also clarify the stability issue, we performed a detailed campaign of simulations, aimed at generalizing the results of Fig. 8. Results of the investigations are reported in Figs. 9 and 10, whereĪ is represented versus ǫ, for two selected b 0 amounts. The LBA initial condition is always stable under Vlasov dynamics. It is in fact a global maximum of the Lynden-Bell entropy, which, in this respect, proves adequate to describe the system at hand. The observed LB0 evolution is by far more complex. For large values of the energy, LB0 is stable. The stability is eventually lost when reducing the energy parameter: A transition materializes and the LB0 evolves towards the LBA state. In the vicinity of the transition, the LB0 initial condition approaches an asymptotic configuration, which slightly differs from the LBA one, and possibly results from a balance between opposing dynamical strengths. The WB evolution mimics that of LB0, the two curves returning a pretty close correspondence. In practice, during a short transient, an initially bunched WB expands (almost) ballistically, the particles being essentially transported by their own initial velocities, so visiting the whole interval [−π; π]. The obtained distribution can be approximated by a homogeneous LB0 state (the field has not yet developed, its intensity being effectively negligible), which in turn explains the observed correspondence. We shall however emphasize that this mechanism applies to relatively small b 0 amounts. This fact is testified in Fig. 10: The agreement between LB0 and WB evolution is shown to worsen, when compared to that of Fig. 9. This is understood as follows: Starting from a high degree of bunching, the induced field opposes the natural ballistic contribution, by further enhancing the tendency to form a coherent clump of particles.
In summary, our calculation returns two stationary points of the Lynden-Bell entropy. The first, which we termed LBA, corresponds to a inhomogeneous (laser on) configuration and it is a global maximum of the entropy. The second, labelled LB0, is homogeneous (laser off). For sufficiently large energies, the system can be locally trapped in the vicinity of the LB0. This happens if the system is initiated close enough to a LB0 state, as e.g. in the case of WB with moderate b 0 values. For smaller ǫ, the LB0 loses stability and the system departs towards the entropically favoured LBA state.
Having detected no additional stationary points of the fermionic entropy, other than LB0 and LBA, we interpret LB0 and LBA as a saddle point and a global maximum, respectively. While LB0 could also be in principle a global minimum, this hypothesis is invalidated by the fact that there exists (at least) a function f , compatible with the system dynamics, which yields a lower entropy value. Consider in fact : This is not a stationary point of Lynden-Bell's entropy but represents one of the admissible equilibria of the dynamical system (2), in its coarse grained perspective as implied by the theory. The Lynden-Bell entropy associated to Eq. (9) can be straightforwardly computed via Eq. (5) and it is found to be always smaller than S 0 , the value associated to the LB0 configuration, for all b 0 and ǫ. Based on the above, and as previously anticipated, we can exclude the possibility for LB0 to be a global minimum. Following our deductive reasoning, we hence suggest that LB0 is instead a saddle-point and the observed transition is consequently interpreted to stem from a local modification of LB0 stability properties or morphological characteristics (e.g. width/flatness of the stability basin). A detailed analytical characterization of LB0 stability is at present missing and could eventually help clarifying the underlying scenario.

VI. CONCLUSIONS
Using the case of a FEL as a paradigmatic example, in this paper we described the out-of-equilibrium dynamics of a mean-field model for wave-particle interaction. Our numerical investigation moves from the Vlasov version of the model, which rigorously applies to the continuous limit and is believed to constitute the correct interpretative framework to elaborate on the QSS peculiarities. Working within this context, and assuming a specific class of initial conditions, we identified a switch between different macroscopic regimes. Such a transition is ultimately controlled by the nominal energy value and by the initial particles bunching, in qualitative agreement with what was previously observed for the HMF model.
The Lynden-Bell violent relaxation theory is developed with reference to the FEL setting to quantitatively substantiate our findings. We numerically characterized the stability of the stationary points of the Fermionic entropy functional. A sudden change in their characteristic is related to the occurrence of the observed transition, supporting the adequacy of the violent-relaxation theory to describe the states reached by the dynamical system.
As a final comment, and beside stressing the unifying picture that is here brought forward, we emphasize that the transition here predicted can be in principle observed in real devices. We regard this as a rather important series of experiments which could eventually result in a direct proof on the existence of QSS for wave-particle systems.
where the symbol * refers to the complex conjugate and in general ω ∈ C. Making use of the above ansatz in the linearized system of equations returns the following consistency equation: often referred to as to the dispersion relation. To determine whether a given distribution f 0 (p) is stable or unstable, one can solve the above dispersion relation and estimate the sign of the imaginary part of ω. Depending on the sign the field grows exponentially (instability) or oscillates indefinitely (stability). If the selected initial condition is parametrized via an adjustable parameter, one can then calculate the corresponding theshold value which discriminates between stable and unstable regimes. This analytical procedure can be persecuted in simple cases, as the one addressed in this paper (the waterbag).
For more complicated situations one can resort to the celebrated Nyquist method, first introduced in plasma physics [25] (see also [26]). For the case at hand, f 0 (p) takes the form: where Θ stands for the Heaviside function. Inserting (A9) into (A8), carrying out the integral explicitly and looking for the value of ∆p which sets the transition between complex and real ω, leads to ∆p ≃ 1.37 or equivalently ǫ = (∆p) 2 /6 ≃ 0.315.