Quantum and classical phase-space dynamics of a free-electron laser

C. Moritz Carmesin ,1,2 Peter Kling ,3,1 Enno Giese ,1 Roland Sauerbrey,2 and Wolfgang P. Schleich 1,3,4 1Institut für Quantenphysik and Center for Integrated Quantum Science and Technology (IQST), Universität Ulm, Albert-Einstein-Allee 11, D-89081, Germany 2Helmholtz-Zentrum Dresden-Rossendorf e.V. Bautzner Landstraße 400, D-01328 Dresden, Germany 3Institute of Quantum Technologies, German Aerospace Center (DLR), Söflinger Straße 100, D-89077 Ulm, Germany 4Hagler Institute for Advanced Study at Texas A&M University, Texas A&M AgriLife Research, Institute for Quantum Science and Engineering (IQSE) and Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843, USA

In a quantum mechanical description of the free-electron laser (FEL), the electrons jump on discrete momentum ladders, while they follow continuous trajectories according to the classical description. In order to observe the transition from quantum to classical dynamics, it is not sufficient that many momentum levels are involved. Only if additionally the initial momentum spread of the electron beam is larger than the quantum mechanical recoil, caused by the emission and absorption of photons, the quantum dynamics in phase space resembles the classical one. Beyond these criteria, quantum signatures of averaged quantities like the FEL gain might be washed out.

I. INTRODUCTION
Usually, an FEL is considered as a device that can be fully described within classical physics. However, there exists a "quantum regime" [1][2][3] where quantum mechanics is indeed mandatory for an accurate description of the FEL dynamics.
In this article, we analyze the transition from quantum to classical in a low-gain FEL by contrasting the dynamics of an electron in phase space with the corresponding classical description. We find that the occurrence of quantum effects depends on the quantum mechanical recoil, caused by the absorption and emission of photons: A small recoil energy, compared to the coupling to the fields, and a small recoil momentum, compared to the initial momentum spread, are necessary to observe a classical evolution of the Wigner function. Furthermore, we study quantum corrections to the FEL gain.

A. Historical overview
In his ground-breaking article [4] in 1971, John Madey had already formulated a quantum theory for the FEL even before the first classical theories emerged [5,6]. His approach relies on a perturbative solution for the electron wave function and was lateron refined for example in Refs. [7][8][9][10].
On first sight, this model perfectly describes the transition between the quantum and the classical regime of the FEL: in the former one the resonances for photon emission and absorption are well separated; by taking the limith → 0, however, these resonances overlap and the difference of photon emission and absorption turns into a derivative. Since in this case all terms with Planck's constanth drop out, the expression for the FEL gain is purely classical [11].
Nevertheless, it was soon realized [12] that the correct description of this transition is more subtle. An electron in the FEL emits many photons during the interaction. The firstorder perturbation theory in Madey's work, however, includes only single-photon processes. Although Madey derived the correct result for the FEL gain, we strictly speaking cannot employ his method.
This puzzling result has led to a variety of different approaches towards a quantum theory for the FEL [10,[12][13][14][15][16][17][18]. For example, in Ref. [14], it was argued that the higher-order contributions due to multiphoton transitions cancel similar to the elementary model of a classical and fixed electron current that is coupled to a quantized radiation field [19].
In Ref. [10], however, the problem was considered from a more practical point of view. Although the perturbative expansion of the quantum state does not converge, the corresponding expansion for the observable of interest nevertheless may converge. Hence, in such a situation, the results from standard perturbation theory can be used to calculate the corresponding expectation value, regardless of the underlying physical mechanism.

B. Wigner function
A new facet to these topics was added by formulating a quantum theory of the FEL in terms of the Wigner function [20][21][22]. Indeed, the Wigner function is a perfect choice to study the transition from classical to quantum: On one hand it contains all information of the quantum state [23]. One the other hand, this description of quantum mechanics [24][25][26][27][28] is as close as one can get to classical phase space [28].
The dynamics of an electron in the FEL can be interpreted as one-dimensional motion of a particle with mass m in a periodic potential [6]. For convenience we choose the dimensionless representations θ ≡ 2kz + const. and ℘ ≡ p/ √ U 0 m for the position z along the wiggler axis and its conjugate momentum p, respectively, see Appendix A. Here, k denotes the wave number of the laser field as well as of the wiggler field in the Bambini-Renieri frame [29], while U 0 is the height of the periodic potential and includes the amplitudes of both fields.
The quantum state of a single electron in Wigner representation W = W(θ,℘; τ ) evolves according to the quantum Liouville equation [25,27] from Eq. (B1), where τ ≡ 2kt √ U 0 /m denotes a dimensionless version of the time t. The left-hand side of this partial differential equation describes the free time evolution while corresponds to the periodic potential. For a derivation see Appendices A and B.
Here we have introduced the parameter α ≡ U 0 /(2hω r ) as the ratio of potential height and recoil energyhω r . The recoil frequency ω r ≡ (2hk) 2 /(2mh) is associated with the recoil 2hk the electron experiences when scattered from a laser and a wiggler photon. While the recoil itself is the origin of gain in the FEL [14], its discrete nature is responsible for the emergence of quantum effects [2]. Finally, the dynamics of the normalized amplitude ε = ε(τ ) of the laser field follows from Maxwell's equations resulting in a semiclassical model for the FEL dynamics.

C. Quantum versus classical
In Ref. [2], we have identified α as the quantum parameter that governs the transition from the classical limit α 1 to the quantum regime α 1 of the FEL. For large values of α, the quantum Liouville equation at first sight reduces [20] to the Boltzmann equation [5] ∂ ∂τ for a classical phase-space distribution function f cl = f cl (θ,℘; τ ). While the free part is the same as in the quantum model, the potential term L (1) cl contains only one derivative with respect to ℘ instead of an infinite sum. Thus we apparently obtain L (1) → L (1) cl in the limit α → ∞. However, the situation is more involved: only increasing α does not necessarily lead to the classical regime. In Eq. (2), we cannot simply truncate the series if we do not take the magnitude of the derivatives of W into account [28,30].
Let the Wigner function be characterized by the width ℘ ≡ p/ √ U 0 m in momentum space. Due to the derivatives, powers of p will appear in the denominator when L (1) acts on W. We thus infer the scaling for the mth term of the series after recalling the definition of α. Indeed, if p is of the order of the recoilhk (that is ℘ ∼ α −1/2 ), all higher contributions of the series are of the same order as the "classical" contribution W/ ℘ (for m = 0) and we must not neglect them. Hence, to obtain the classical limit it is not sufficient that the quantum mechanical recoil is small compared to the height of the potential. In addition, it has to be small compared to the momentum spread of the electron beam. For a pure state, for example a Gaussian wave packet [31,32], decreasing the momentum uncertainty is inevitably connected to an increase in position uncertainty. In a true classical limit, however, momentum and position uncertainties can be chosen independently of each other [28]. This feature only emerges if the quantum state of an electron is described by a statistical mixture rather than by a pure state. In other words, the classical uncertainties for momentum and position are larger than the intrinsic quantum ones from the uncertainty principle.
It is convenient to assume that the initial Wigner function for the electron beam is given by the product of an arbitrary momentum distribution ρ = ρ(℘) and a uniform distribution in θ direction. This choice in position space is in accordance with classical FEL theory [5,33,34] since we cannot control the exact positions of each one of the electrons that are distributed over several wiggler wavelengths [33,35]. In contrast to a pure state, we can choose the width of ρ without affecting the distribution for θ [28]. However, the evolution of a "classical state" does not ensure that we are in the classical regime for all times [28,36]. Local modulations of the Wigner function on the scale of hk may emerge so that quantum effects reappear. We therefore expect that ultimately other effects enforce classicality, namely sources of decoherence [37], for example induced by space charge or spontaneous emission [38]. Alternatively, one can argue that a "classical observer" is unable to resolve these fine quantum signatures and is therefore limited to an averaged measurement result [39].
We emphasize that the use of a statistical mixture does not imply a many-particle theory. In the low-gain regime [33], the motions of all electrons decouple from each other and every electron interacts separately with the fields. Each electron in a bunch is a copy of itself, initially distributed according to W(θ,℘; 0), and we interpret the N single-particle interactions as N repetitions of the same experiment. At the end we calculate the observable quantities, like the FEL gain, by averaging over all possible outcomes weighted by the timeevolved Wigner function W(θ,℘; τ ).
For example, we derive in Appendix A the semiclassical equation of motion for the dimensionless electric field ε with the constant χ that includes the initial electron density n e and the initial strengths of wiggler and laser fields. The dynamics of ε follows from Maxwell's equations with the electron current being determined by the Wigner function. From Eq. (5), it is evident why it is difficult to observe quantum effects in the FEL radiation: even if the Wigner function shows distinct quantum signatures, these features might be washed out when we average over position and momentum.

D. Outline
The remainder of this article is structured in the following way. In Sec. II, we compare the time evolution in phase space obtained from a quantum and a classical theory for the same initial state. We first employ a perturbative expansion valid for the small signal-limit before we resort to numerics in order to treat also longer interaction times. Next, we consider in Sec. III the gain of the FEL. Again, we cover the small-signal limit as well as the evolution for longer times. Finally, we discuss our results in Sec. IV.
For completeness we derive in Appendix A the semiclassical model of the FEL used throughout this article. Further, we have moved the perturbative calculation for the Wigner function to Appendix B. In Appendix C, we introduce energy eigenstates in terms of Mathieu functions and employ them to derive a formal expression for the Wigner function used for the numerical computations.

II. EVOLUTION OF WIGNER FUNCTION
The quantum Liouville equation constitutes a partial differential equation to which we cannot find an exact, analytic solution in general since it contains infinitely many derivatives with respect to ℘. However, before we turn to a numerical approach, we asymptotically solve Eq. (1) for small times. This short-time limit corresponds to the small-signal regime of an FEL.

A. Small-signal limit
In Appendix B, we perform the asymptotic expansion of W and solve Eq. (1) order by order. This perturbative approach resembles the procedure in Ref. [5] for the Boltzmann equation and thus differs from the usual perturbative solution of the Schrödinger equation [4,[7][8][9], which is restricted to single-photon processes. While the former procedure is allowed as long as τ 1, the latter one is only valid for √ ατ 1 and thus breaks down quickly in the classical regime due to α 1. We choose the initial state from Eq. (4) and assume that the initial momentum distribution ρ is Gaussian with mean valuē ℘ and standard deviation ℘. In Appendix B, we derive the expression . In all cases, we have chosen the values ε = 1, τ = 0.01,℘ = π , and ℘ = 0.1 for the normalized field amplitude, the dimensionless time, the mean initial momentum, and the initial momentum spread, respectively. Although α = 16 1 indicates that we are close to the classical regime the corresponding curve for W (1) significantly differs from the classical result. Only after increasing α to α = 400, we observe an agreement between the quantum and the classical theory. As apparent from Eq. (6) the important parameter that governs this transition is not α but rather the ratiohk/ p of recoilhk and momentum spread p. Indeed, for α = 16 this parameter is withhk/ p = 1.25 of the order of unity and we cannot neglect quantum corrections while α = 400 leads to the decreased valuehk/ p = 0. 25 1.
for the first-order contribution to the unperturbed Wigner function, where f (1) cl from Eq. (B7) denotes the corresponding solution for the classical Boltzmann equation (3).
The quantum corrections Q = Q(ξ ) are given by the series of Hermite polynomials which depend only on the relative momentum ξ ≡ (℘ −℘)/( √ 2 ℘), but not on position or time. The analysis of Eq. (6) reveals the importance of the ratio of recoil and momentum spread: The terms of the series scale with powers ofhk/ p and only if this parameter is small, each term of the series decreases and the Wigner function W (1) approaches the classical distribution function f (1) cl .
In Fig. 1, we plot W (1) as a function of the relative momentum for different values of the quantum parameter but for a fixed momentum spread and compare it to the classical result f (1) cl . Although α = 16 indicates that we are close to the classical regime, W (1) does not agree with the classical curve. Following the discussion above, the important parameter is not α but ratherhk/ p = 1/ 2 √ α ℘ which is in this case of the order of unity, that ishk/ p = 1.25. Indeed, increasing the quantum parameter to α = 400 leads tohk/ p = 0.25 and thus to negligible quantum corrections in Fig. 1.
In conclusion, the phrase "small recoil" in the present context means: small compared to the initial momentum spread. Only in this limit we can neglect the quantum corrections to the Wigner function.

B. Longer times
The small-signal regime allows for a perturbative treatment, which is not possible for longer times. Hence, we have to determine the time evolution of the Wigner function numerically. In contrast, the solution of the classical equation of motion can be given in a closed form [41]. In Fig. 2, we compare our results for the Wigner function to the classical distribution function for different initial momentum widths ℘ and for different values of the quantum parameter α.
The relevant-classical and quantum mechanical-phasespace dynamics takes place inside the classical separatrix, which separates open and closed trajectories, because here the largest changes of momentum appear. The classical evolution of a phase-space distribution is a rotation inside the separatrix. Since we consider an anharmonic oscillator, the angular velocity depends on position and momentum and decreases going from the center to the separatrix. As a consequence, the phase-space distribution is stretched during the evolution, since the inner parts move faster than the outer ones, see the right column of Fig. 2.
The time evolution of a Wigner function that is initially uniform in space can be expressed as for details see Appendix C 2. This expression can be understood in the following way: The momentum of the electrons changes by integer multiples s of the discrete recoil 2hk. The counter-intuitive occurrence of half-integer multiples s/2 of 2hk in Eq. (7) comes from the interference between two adjacent momentum levels [42]. Note that the momentum changes of s2hk correspond to the emission and absorption of |s| photons each. Therefore multiple scattering events are included in our low-gain model. We interpret the prefactors w s as scattering amplitudes for the shifted parts of the Wigner function. They depend not only on time and position but also on the momentum itself. As a consequence of this momentum dependence, only a fraction of the distribution is selected to participate in the interaction. This effect is also known as "velocity selectivity" in atomic Bragg diffraction [43].
In the quantum regime [2], that is, α 1 and p/2hk 1, the summation in Eq. (7) breaks down to a few terms, as only few momentum levels are involved. In this extreme limit, single-photon processes dominate the dynamics. Since the initial momentum spread is sufficiently small, the shifted distributions do not overlap and hence appear as discrete lines in phase space.
In Fig. 2(b), where the time-evolved Wigner function for α = 1/3 and ℘ = 0.1 is shown, we indeed observe only three momentum levels, that is the initial momentum p = 0 and the levels p = 2hk and p = −2hk, as well as two interference patterns between these levels. From the marginal distribution in Fig. 2(a), we recognize that these intermediate momenta do not contribute to the momentum distribution. As expected from the valuehk/ p = 8.66, there is no agreement with the classical distribution in Fig. 2

(d).
With increasing α, more and more momentum levels intersect with the separatrix and therefore contribute to the sum in Eq. (7), that is, multiphoton processes occur. The phase-space structure starts to resemble the classical dynamics but with an additional fine structure due to the interference between many levels, see Fig. 2(c). Here we havehk/ p = 1.58, which is still outside the range for a classical evolution. From the momentum distribution in Fig. 2(a), we recognize several single peaks corresponding to the distinct momentum levels, while the classical distribution is spread out rather homogeneously over the area enclosed by the separatrix.
We infer from these plots that a large α by itself is not sufficient to obtain a classical time evolution, which underlines the results from the preceding section. A broad initial momentum distribution is also not enough, as we can see in Fig. 2(f) (hk/ p = 0.87) and (m) (hk/ p = 0.43) for α = 1/3. Here the discreteness of the few momentum levels is washed out due to the overlap of the shifted contributions [2]. The shape of each distribution is however quite different from the classical counterpart in Figs. 2(h) and 2(l), respectively, since only few momentum levels are involved.
The combination of a broad initial distribution and the inclusion of many momentum levels, that is a large ℘ and a large α, eventually leads to a Wigner function which seems to resemble the classical distribution function. By comparing Fig. 2(l) with Fig. 2(k), where the valuehk/ p = 0.08 matches the condition for a classical evolution, we can hardly recognize any difference between these two distribution functions. The momentum distributions agree even better, see Fig. 2 At first sight, Fig. 2 completely confirms our expectations: For a small value ofhk/ p the Wigner function resembles its classical counterpart. However, if we go beyond this purely visual comparison, we indeed observe significant differences. For this purpose, we consider the numerical deviation W(θ,℘, τ ) − f cl (θ,℘, τ ) of both distributions. From Fig. 3, we notice that the differences increase in the course of time, even if the conditions for a classical evolution are satisfied, that is α = 10 and ℘ = 2 leading tohk/ p = 0.08.
For short times, the deviations have a low amplitude, see the color bar of Fig. 3(a). The structure in momentum direction can be identified as the Hermite polynomials of the approximated analytical solution from Eq. (6). For longer times, the anharmonic dynamics leading to dispersion of the distribution in proximity to the separatrix becomes more prominent. At the points in phase space where the distribution is narrow, the deviation increases by orders of magnitude, see In all cases, we have chosen the time τ = π , which classically corresponds to half of a rotation near to the center of the potential. Each row corresponds to a different initial momentum width, that is ℘ = 0.1, 1.0, and 2.0. The marginals on the left show the initial momentum distribution (dotted) and the evolved ones for α = 1/3 and 10 as well as the corresponding classical result (solid lines). The second and third column show the evolved Wigner functions for α = 1/3 and α = 10, the rightmost column displays the classical phase-space distribution. In the case of small α and ℘ the discrete structure in momentum is visible in (b), where we havehk/ p = 8.66. With an increased α more levels are involved in (c) and the shape of the corresponding classical distribution in (d) starts to appear, even though there are a lot of interference structures. This behavior is consistent with the value ofhk/ p = 1.58, which is of the order of unity. A broader initial momentum distribution washes out the discrete levels, but still the Wigner functions for α = 1/3 in (f) withhk/ p = 0.87, and in (j) withhk/ p = 0.43 are quite different from their classical counterparts in (h) and in (l), respectively. For large α = 10 and ℘ = 1, that is,hk/ p = 0.16, there are only a few interference fringes left in (g). An even wider initial distribution with ℘ = 2, corresponding tohk/ p = 0.08, makes them nearly invisible such that Wigner function in (k) and classical distribution in (l) resemble each other. the higher-order derivatives only for some time. Consequently, every initial distribution that evolves coherently sooner or later becomes nonclassical. Decoherence effects like spontaneous emission or space charge may impede or even remove the appearance of the quantum features in the distribution and hence might lead to a classical evolution [37].
Even though all existing x-ray FELs operate in the highgain regime, we rely on their parameters for an estimation of . For a small time τ = π/12, the interference fringes have a width proportional to the initial distribution, but a low amplitude. In the course of time, the amplitude grows (note the different scaling of the colormap for each plot) due to the increasing magnitude of the higher-order derivatives of the Wigner function. Those, in turn, are a consequence of the narrow width in momentum space that appears first in the proximity to the separatrix due to dispersion. The appearing interference fringes are also narrow, leading to a self-reinforcing effect. However, the amplitude of the deviation remains in the depicted case about one order of magnitude smaller than the maximal values of W in Fig. 2(k). the decoherence timescales, since they are to date the only FEL devices with large recoils that approach the quantum regime. We emphasize that the high-gain regime is strictly speaking not covered by our single-electron model.
For the European XFEL [44], we obtainhk/ p = 0.014, a value that is still much smaller than unity but at least corrections in the Wigner function due to the discrete momentum steps might become conceivable. The interference pattern in Fig. 3(c) emerges at τ ∼ = π , that is the typical time for which the FEL gain saturates [45]. For the European XFEL the saturation length amounts to L sat = 133 m corresponding to a saturation time T sat = 440 ns.
The corresponding timescales of possible decoherence mechanisms [38] are calculated [46] to T se = 0.7 ns (corresponding to the length L se = 0.2 m) for spontaneous emission and T p = 1 μs (corresponding to the length L p = 300 m) for space charge, respectively. The comparison of these three timescales reveals that, while we can neglect space charge effects, spontaneous emission occurs long before an interference pattern can emerge. Hence, we expect interference effects to be suppressed in state-of-the-art machines.
In order to quantify the deviation of quantum and classical phase-space evolution we introduce the quantity which defines a Hilbert-Schmidt-like distance [47] of the Wigner function to the classical distribution, normalized to values between zero and unity. The smaller the value of d cl , the closer is the Wigner function to the classical distribution. Due to the integration over the whole phase space, d cl constitutes a global measure in contrast to the local modulations that we observe in Fig. 3.
The introduction of d cl gives us a convenient tool to study the influence of the different parameters on the similarity of quantum and classical evolution. In Fig. 4(a), we plot the time evolution d cl = d cl (τ ) for different values of α and ℘ used in Fig. 2. We observe that the distance d cl increases rapidly for larger values ofhk/ p, which matches the results from the perturbative treatment. Moreover, we notice that for small values ofhk/ p, the distance is very small for short times, see the curves corresponding to ℘ = 2. This behavior fulfills our expectations of a nearly classical dynamics. However, after some time the quantum features in the Wigner function become more prominent and the distance suddenly increases before it saturates. The saturation value is still about an order of magnitude below the possible maximum.
In Figs. 4(b) and 4(c), we plot d cl as a function of α and ℘ for two different times τ . In both panels we observe the behavior that small values of α and ℘, see the lower left corner, lead to a large distance while it decreases when these parameters are increased, corresponding to a more classical evolution in phase space. The contour lines approximately correspond to hyperbolae of constanthk/ p. At the earlier time τ = π/2, shown in Fig. 4(b), the distance quickly falls off when approaching the upper right corner. For the later time τ = π depicted in Fig. 4(c) the distance from a classical evolution declines much slower, since the small at short times for small values ofhk/ p. After a longer time, the distance suddenly increases even for smallhk/ p but saturates afterwards. In (b) and (c), we draw d cl as a function of the initial momentum width ℘ and of the quantum parameter α at times τ = π/2 and τ = π , respectively, corresponding to the grey vertical lines in plot (a). For small values of α and ℘, that is the lower left corner, d cl is maximized and close to unity, while it decreases for large values of α and ℘ towards the upper right corner, that is, the evolution becomes more classical. At the earlier time τ = π/2, the decrease is steep, while it is more gradual at the later time τ = π . Even if the distance is very small at an earlier time, this does not mean that it will be small also for later times.
Wigner function and the classical distribution function differ after some time, as seen in Fig. 4(a).

III. FEL GAIN
So far, we have only discussed the motion of the electron. However, in an experiment we are mainly interested in the dynamics of the laser field, which is why we investigate in the following quantum effects of the FEL gain.

A. Small-signal limit
Inserting the expression from Eq. (B5) for W (1) into the equation of motion (5) for the laser field, and averaging over the phase θ yields a linear differential equation for ε which can be straightforwardly solved. This solution reads [48] where we have introduced the gain of the laser field in the small-signal limit.
To derive analytical expressions we restrict ourselves to the extreme cases of a cold and a warm electron beam [45] defined by the limits ℘τ 1 and ℘τ 1, respectively. Making the approximation ρ(℘) ∼ = δ(℘ −℘) leads us from Eq. (10) to the FEL gain for a cold beam, while the assumption sinc 2 (℘τ /2) ∼ = 2π δ(℘)/τ yields the corresponding quantity in the warm-beam limit.
In contrast to the Wigner function and the warm-beam case, the quantum corrections for a cold beam do not scale with powers ofhk/ p but instead with even powers of the recoil parameter, defined as ω r t = τ / √ 4α. We note that the relationh connects the three relevant parameters. For a cold beam with ℘τ 1, the recoil parameter ω r t is always smaller than the ratio of recoil and momentum spread, that is ω r t hk/ p. Hence, quantum effects apparent in the Wigner function are washed out due to the averaging for the FEL gain. In the opposing limit of a warm beam, Eq. (13) predicts thathk/ p in turn is smaller than the recoil parameter ω r t dominating the cold-beam limit, that is hk/ p ω r t due to ℘τ 1. In Fig. 5, we compare the classical gain function to the one including quantum corrections for a cold and a warm electron beam, respectively. We observe the same qualitative behavior in both cases, that is a slight shift of the positions for maximum and minimum gain and a decrease of its magnitude for increasing values of the recoil parameter ω r t. However, while moderate quantum corrections in the cold-beam case emerge already for ω r t = 1, it is not until ω r t = 3 that they appear for a warm beam and they remain small even for ω r t = 7. 023027-7 FIG. 5. Quantum corrections to the FEL gain: we have drawn the small-signal gain for a cold (above), Eq. (11), and a warm electron beam (below), Eq. (12), as functions of℘τ = 2kpt /m and℘/ ℘ = p/ p, respectively. In the cold-beam case, ℘τ 1 we compare the classical gain (black line) to gain curves including the lowestorder quantum corrections [m = 1 in Eq. (11)] for the values ω r t = 1 (blue, dotted line) and ω r t = 2 (red, dashed line), respectively, of the recoil parameter ω r t. We observe that for increasing values of ω r t the positions of minimum and maximum are slightly shifted to the left and to the right, respectively, while their magnitude is decreased. The behavior in the warm beam case with ℘τ 1, is qualitatively similar. However, for our-already moderate-choice of ℘τ = 10 visible quantum effects [again in lowest order, that is m = 1 in Eq. (12)] occur not until ω r t = 3 (blue, dotted line) and are still quite small for ω r t = 7 (red, dashed line). Hence, quantum corrections in the warm beam case are small in comparison to a cold beam in accordance to Eq. (13). We conclude that a large momentum spread suppresses quantum effects in the FEL gain.
Hence, quantum effects are suppressed for a warm beam in comparison to a cold beam in accordance with Eq. (13).
We conclude that the averaging over all momenta as well as a large momentum spread prevents the appearance of quantum effects in the FEL gain, at least in the small-signal limit.
The expression in Eq. (11) is not a new result. It is in principle the same formula derived by Madey and many others [4][5][6][7][8][9][10]. For m = 0, we indeed recover from Eq. (11) the classical result-Madey's famous gain formula. In contrast to these earlier works, our approach is not restricted to singlephoton processes but instead generalizes the procedure [5] for the classical Boltzmann equation (3).

B. Longer times
It is possible to relate the gain of the laser field to the change of the mean momentum of the electrons. Inserting the quantum Liouville equation (1) into the time derivative of the expectation value of the electron momentum and making use of Eq. (5), we find after assuming ε(τ ) ∼ = 1. Identifying the gain as the relative change of the field in accordance with Eq. (9), we arrive at the expression Having already calculated the Wigner function, this expression for the gain can easily be evaluated. The classical gain can also be inferred from Eq. (14) with the help of the classical distribution function f cl . In Fig. 6, we plot the time evolution of the FEL gain, obtained from Eq. (14) and the numerical solution of the Wigner function. Because the gain vanishes for℘ = 0, we chose a nonzero initial mean momentum (℘ = 1.6 in the plots).
For longer times, the gain does not increase monotonically, but exhibits oscillations in a saturation regime. The classical time evolution of the phase-space distribution yields a rotation inside the separatrix, leading to an oscillating change of the mean momentum and consequently also to the oscillations of the gain. This behavior was for example discussed in Ref. [45].
In Fig. 6, we do not present the case α 1, since here the gain shows a behavior completely different from the classical description [49]. With increasing α the frequency of oscillation approaches the classical one such that the quantitative differences between quantum and classical evolution become less prominent and only appear at relatively large times, such that both theories lead to a similar behavior of the gain. This convergence becomes apparent in Figs. 6(a) and 6(b). A larger ℘ further reduces the difference between the two gain curves.
In contrast to the phase-space distance d cl , which increases rapidly in the beginning for an initially small momentum spread, the classical and quantum gain resemble each other very well for short times. We have seen this behavior already in the section above, where the quantum corrections to the gain scaled asymptotically with ω r t = τ / √ 4α. The good agreement between quantum and classical theory arises from the fact that the gain depends on the change of mean momentum, that is a quantity where the differences between classical and Wigner distribution are averaged out, even if they are as substantial as in Fig. 2(c), while the distance takes all deviations into account. Further, we observe that the quantum corrections become smaller for decreasing values of 1/(2 √ α ℘) =hk/ p.

IV. CONCLUSIONS
By comparing the time evolution of the Wigner function and of the classical phase-space distribution in the lowgain regime, we have analyzed the effects of the initial momentum width ℘ and the quantum parameter α. The asymptotic treatment valid for the small-signal regime, that is short times, allowed us to identify the critical quantity √ α ℘ governing the transition from quantum to classical. The numerical solution for longer interaction times confirmed that neither the variation of ℘ nor α by themselves are sufficient to obtain matching dynamics of the Wigner function and the classical distribution function, but only the combination of both. Furthermore, we observe always that, after sufficient time, both distributions differ due to nonlinear effects imposed by the anharmonic potential.
For the FEL gain, the quantum effects are much less prominent and hence hard to observe in comparison to the Wigner function. At the highly relativistic electron energies of an FEL, we cannot determine the Wigner function sufficiently precisely from experiment. However, our theory is not limited to this specific energy range and can also be applied to electrons at lower energies. Other experimental situations, for example Kapitza-Dirac [50], provide a more favorable basis for the obvervation of the quantum to classical transition, since here the momentum level structure and the oscillations between the levels are experimentally resolvable [51][52][53]. There is even the possibility to reconstruct the Wigner function from measurement data [54].
To extend our model to the high-gain regime [55], a many-electron description including a varying laser field becomes necessary. However, the quantum limit of the high-gain regime also reduces to a system of two momenta for each electron, such that the process is dominated by single-photon transitions. The classical regime, for low-as well as for highgain FELs, differs from the quantum limit by the emission of multiple photons per electron. In the high-gain regime, many electrons simultaneously interact with the light fields, and thus they may become entangled [3]. However, if we neglect these quantum correlations between the individual electrons, the high-gain regime mainly differs through an rapidly increasing laser field from the low-gain limit. Hence, we expect that the parameterhk/ p plays a similar role in both regimes.

ACKNOWLEDGMENTS
We thank P. Anisimov, A. Debus, M. Efremov, M. Freyberger, A. Gover, Y. Pan, and K. Steiniger for many fruitful discussions. W.P.S. is grateful to Texas A&M University for a Faculty Fellowship at the Hagler Institute of Advanced Study at Texas A&M University, and Texas A&M AgriLife Research for the support of his work. The Research of IQ ST is financially supported by the Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg.

APPENDIX A: DERIVATION OF MODEL
In this Appendix, we derive our model for the FEL dynamics in terms of the Wigner function. While we describe the motion of the electron quantum mechanically, the laser field is treated as a classical quantity in analogy to semiclassical laser theory [56].

Pendulum Hamiltonian
The one-dimensional motion of an electron with mass m and charge e in the FEL is dictated by the Hamiltonian [35] H =p withẑ andp being the position and the momentum operators along the wiggler axis. We assume that the vector potentials A L = A L e x and A W = A W e x of the laser and the wiggler field, respectively, are linearly polarized perpendicular to the longitudinal axis. Moreover, we model their amplitudes as plane waves [33] A L (z, t ) = E L (t ) ω sin (ωt + φ L (t ) − kz) and (A2) where ω and k = ω/c denote the frequency and wave number, respectively, of the two counterpropagating fields, while c gives the velocity of light. We consider the reference 023027-9 frame [29] where the frequencies of the two modes coincide and thus the electron interacts with a standing light field. For resonant interaction the electron motion in this frame is always nonrelativistic [2] justifying the Hamiltonian in Eq. (A1). The amplitude and phase corresponding to the electric field of the laser mode, E L = E L (t ) and φ L = φ L (t ), respectively, are slowly varying with time. In contrast, we assume that the corresponding quantities for the strong magnetic field of the wiggler, B 0 and φ W ≡ 0, are constant.
By inserting the fields from Eqs. (A2) and (A3) into Eq. (A1) and neglecting rapid oscillations with 2ω in a rotating-wave-like approximation [45], we arrive at the pendulum Hamiltonian [6,29] H =p 2 2m for the electron dynamics in the FEL. Here we have introduced the amplitude of the potential as well as the dimensionless electric field ε(t ) ≡ E L (t )/E 0 normalized with respect to the initial field amplitude E 0 ≡ E L (0) before the electrons enter the wiggler.
In the low-gain regime of FEL operation the value of ε is always close to unity.

Laser field couples to electron current
In contrast to the quantized electron, we describe the laser field as a classical quantity that evolves according to Maxwell's equations. Hereby we strongly follow the lines of semiclassical laser theory [56]. The classical wave equation [57] couples A L to the x component j el of the electron current with μ 0 denoting the vacuum permeability. Inserting Eq. (A2) into Eq. (A5), projecting on the laser mode, performing the slowly-varying phase and amplitude approximation [35,56], and taking the imaginary part leads to the equation of motion for the laser amplitude that includes the vacuum permittivity ε 0 = 1/(c 2 μ 0 ). We observe that only the Fourier component of the current that corresponds to the laser mode appears in the equation of motion, with L z denoting the longitudinal extend of the quantization volume. In the course of the slowly-varying phase and amplitude approximation we have neglected terms withË L ,φ L ,φ LĖL ,φ 2 L ,φ L J L , andJ L in accordance with Ref. [56].
Since the gradient of the electron wave function ψ = ψ (z, t ) in x direction vanishes, we deduce that is the x component of a current that satisfies a continuity equation following from the Schrödinger equation of a charged particle in the electromagnetic field. Here, N gives the number of electrons in the bunch while L x and L y describe normalization lengths. We note that we have neglected a contribution with A L A W . By inserting Eq. (A8) into Eq. (A7) and neglecting a rapidly oscillating term, we finally arrive from Eq. (A6) at the equation of motion for the field amplitude coupling to the expectation value of cos 2kẑ. Here we have introduced the initial density n el ≡ N/V of the electrons after identifying the quantization volume V ≡ L x L y L z with the volume of the electron bunch.

Formulation with Wigner function
With the pendulum Hamiltonian for the dynamics of the electron and the equation of motion for the laser amplitude we have obtained the two important relations of our semiclassical theory for the FEL. In order to illuminate the transition to a purely classical theory [5], we introduce the Wigner representation [27] for the density operatorρ of an electron. This function depends on the position z and its conjugate momentum p which together form the Wigner phase space. It is convenient to use the dimensionless variables which correspond to position, momentum, and time, respectively. Accordingly, we write the Wigner function W = W(θ,℘; τ ) in its dimensionless form, which, for the correct normalization, includes the factor for the Wigner function also known as quantum Liouville equation [27]. While the left-hand side of this equation corresponds to the free motion of the electron, the right-hand side emerges due to the periodic potential.
In the Wigner description of quantum mechanics, the expectation value of an observable is formed by weighting the Weyl representation of the corresponding operator with the Wigner function and integrating over phase space [27,28]. For the expectation value in Eq. (A9) this procedure yields the relation d dτ ε(τ ) = −χ dθ d℘ W(θ,℘; τ ) sin θ for the dynamics of the dimensionless field amplitude ε with characterizing the coupling between field and electrons.

APPENDIX B: PERTURBATION THEORY
This Appendix is devoted to the solution of the quantum Liouville equation by means of perturbation theory [25] which is valid in the small-signal limit of FEL operation.

Structure of equation
We first cast Eq. (A11) into the form where we have defined the operators for the free motion and the potential, respectively. By performing a Taylor expansion of the shifted Wigner functions in Eq. (A11) in powers of 1/(2 √ α) we identify L (1) as differential operator containing infinitely many derivatives with respect to momentum.
We interpret the right-hand side of Eq. (B1) as a perturbation to the free motion and make the asymptotic expansion W ∼ = W (0) + W (1) + W (2) + · · · in analogy to Ref. [25]. By inserting this expansion into Eq. (B1) we obtain the equations for zeroth and higher orders (n > 0), respectively, which we solve iteratively.

Zeroth order
The equation in lowest order corresponds to a free particle and is solved by In our case, the lowest-order contribution of the asymptotic expansion equals the initial distribution because it is independent of θ according to Eq. (4).

First order
We now turn to the first-order calculations. The formal solution of Eq. (B2) for n = 1 reads where we have recalled from Ref. [25] the Green's function for the unperturbed quantum Liouville equation.
Inserting the lowest-order solution W (0) from Eq. (B3) into Eq. (B4) and evaluating the operators and integrals yields the result for the first-order term W (1) , which still contains infinitely many derivatives of ρ. We note that we have treated the relative field change ε as a constant in this procedure, that is ε(τ ) ∼ = ε(τ ), in accordance with our low-gain approach.
The assumption that the initial momentum distribution ρ is Gaussian with mean value℘ and standard deviation ℘ leads to a series of Hermite polynomials of odd order. From this series, we finally arrive at the closed expression for W (1) after using a generating function for Hermite polynomials [58] as well as 4α ℘ 2 −1/2 =hk/ p.

Connection to classical theory
We briefly compare our results for the Wigner function to classical FEL theory [5] in terms of the classical distribution function f cl = f cl (θ,℘; τ ) in phase space. In the equation of motion for f cl , we simply replace L (1) in the quantum Liouville equation by L (1) cl ≡ −ε(τ ) sin θ ∂ ∂℘ and arrive at the colissonless Boltzmann equation (3).
Consequently, the procedure for finding a perturbative solution of f cl is analogous to the approach for W in the preceding section. For the same initial distribution, we thus obtain the first-order contribution [59] f (1) cl (θ,℘; τ ) = −ε(τ ) τ cos (θ − ℘τ ) − cos θ ℘τ for the classical distribution function with H 1 denoting the Hermite polynomial of first order evaluated at ξ ≡ (℘ −℘)/ √ 2 ℘ . By comparing Eqs. (B6) to (B7), we realize that the Wigner function and the classical distribution just differ by a momentumdependent contribution Q. These "quantum corrections" are given by the series of Hermite polyniomials H 2m+1 of odd order that arise from the derivatives with respect to ℘ in Eq. (B5).

APPENDIX C: EXACT TIME EVOLUTION OF THE WIGNER FUNCTION
Since the infinitely many derivatives in the quantum Liouville equation (1) impede the solution in a closed form, we resort to solving the Schrödinger equation. From the resulting solutions, we can construct the time-evolution operator. We then use this operator to obtain the time evolution of the density operator from which the Wigner function is computed.

Schrödinger equation and Mathieu functions
At first, we determine the eigenfunctions and energy eigenvalues of the Hamiltonian (A4). They directly lead to the time evolution of an arbitrary quantum state subject to the periodic potential U 0 cos θ .
We write the stationary Schrödinger equation in position space with the dimensionless position θ and the quantum parameter α and arrive at where we have introduced the dimensionless energy E scaled in units ofhω r . This equation is known as Mathieu equation [60,61]. Due to the periodicity of the potential, we restrict the range of θ to the interval [0, 2π ].
We use here only a subset of the solutions of Eq. (C1), namely, the bounded Mathieu functions u(α, θ ) ≡ me ν (α, θ ), since the other solutions diverge for θ → ±∞ and hence do not describe a physical quantum state. The functions me ν are, in general, not expressible in terms of standard functions and are 2π/ν-pseudoperiodic. That is, they obey the relation The parameter ν ∈ R is hence determined by the particular choice of periodic boundary conditions [62]. Connected to each solution is the characteristic value E ν (α), that is, the eigenenergy of that state. By solving the time-dependent Schrödinger equation, we find describing the time evolution of an energy eigenstate.
In order to attribute a physical meaning to the parameter ν, we recall that the periodicity of a momentum eigenstate is given by the De Brogli wavelength λ DB = 2πh/p, which translates in the dimensionless coordinates to 2π √ α/℘ = 2π/ν. Hence, we identify the parameter ν of the Mathieu functions as the initial momentum. Thus, the initial momentum by itself determines the boundary conditions we have to impose.
In order to explicitly find the solutions of Eq. (C1), we exploit the periodicity of the potential with the help of Bloch theory. Accordingly, the solutions can be represented in a special form of Fourier series, given by where the c ν r (α) denote the Fourier coefficients. Here the series itself has the periodicity of the potential, while the pseudoperiodicity is achieved only through the prefactor e iνθ . For the sake of readability we will from now on suppress the parameter α.
By inserting Eq. (C3) into the Schrödinger equation, Eq. (C1), we obtain which is an infinite set of coupled algebraic equations for the coefficients c ν r . This representation is equivalent to an infinitedimensional matrix eigenvalue problem with the eigenvalues E ν and eigenvectors c ν . After truncation to finite dimension, we solve the problem by numerical means [63]. Recursive insertion of Eq. (C4) into itself leads to a continued fraction expansion for the coefficients [60,64]. Not only can this be used for the numerical computation, but it also is particularly useful to derive asymptotic expressions for the limit α 1. Furthermore, the set {me ν+r (α, θ )} ∞ r=−∞ forms a complete basis of the Hilbert space. This allows us to introduce the concise Dirac notation |r ν , such that we have in position representation. It holds the orthonormality relation r ν |l ν = δ rl , where δ rl denotes the Kronecker delta, and the completeness relation ∞ r=−∞ |r ν r ν | = 1. These properties allow us to expand other functions in terms of Mathieu functions.
In order to calculate the time evolution of a momentum eigenstate, that is a plane wave with initial momentum ν 0 = ℘ 0 / √ α, we expand such a state in terms of the Mathieu functions. In this basis, we read off the time evolution from Eq. (C2). A transformation back into the momentum representation yields where we have introduced the scattering amplitude This representation highlights that an initial momentum p 0 = 2hkν 0 is coupled only to other momenta p 0 + s(2hk), which are separated by multiples of the recoil momentum 2hk. In the limit α 1, most coefficients (except c ν 0 0 and c ν 0 −n for ν 0 ≈ n 2 ∈ Z) vanish and the sums in Eqs. (C5) and (C6) can be reduced to only a few terms. This procedure leads to Rabi oscillations between momentum levels that dominate in the quantum regime of the FEL [2].
In the next step we generalize Eq. (C5) by integrating over all possible momenta as initial states. Thus we obtain the timeevolution operator in momentum representation which can be applied to an arbitrary initial state.

Wigner function
We are now in the position to calculate the time evolution of a quantum state in the Wigner representation. Here we restrict ourselves to an initial state which is initially uniform in space and hence fully determined by the initial momentum distribution ρ(℘), see Eq. (4).
A similar approach to calculate the Wigner function has been persued in Ref. [20], but not for mixed states. Other attempts to use the Mathieu functions to obtain the Wigner function of an electron in a periodic potential [65] do not cover dynamics nor arbitrary momenta. In the context of the FEL, Mathieu functions have been used to derive asymptotic expressions for the FEL gain [13], but not in Wigner phase space.