General method for analyzing three-dimensional effects in free-electron laser amplifiers

Free-electron laser (FEL) configurations in which the parameters of the electron beam vary along the undulator become relevant when considering new aspects of existing FELs or when exploring novel concepts. This paper describes a fully three-dimensional, analytical method suitable for studying such systems. As an example, we consider a seeded FEL driven by a beam with varying transverse sizes. In the context of the Vlasov-Maxwell formalism, a self-consistent equation governing the evolution of the radiation field amplitude is derived. An approximate solution to this equation is then obtained by employing an orthogonal expansion technique. This approach yields accurate estimates for both the amplified power and the radiation beam size. Specific numerical results are presented for two different sets of x-ray FEL parameters.


I. INTRODUCTION
In recent years, there has been renewed interest in freeelectron lasers (FELs) [1], mainly because of their great potential as tunable sources of intense, coherent radiation.In view of the success of facilities such as FLASH [2] and the Linac Coherent Light Source (LCLS) [3], special emphasis is placed nowadays on high-gain, single pass FELs, in which the radiation is generated as a relativistic electron beam passes once through a suitable undulator system.For the majority of existing or planned machines of this type, the wavelength of the output radiation lies in the UV or x-ray region.For such short wavelength FELs, threedimensional (3D) effects, i.e., those due to the transverse size and emittance of the electron beam, become significant and must be taken into account.A comprehensive theory of high-gain FELs, including 3D effects, has been developed in the past and successfully applied in the design of many FEL projects [4][5][6].Strictly speaking, however, this theory is valid only when the various parameters of the undulator and the electron beam do not depend upon the longitudinal position z along the undulator axis.In fact, FEL schemes for which this simplifying assumption is not satisfied can become important when studying new aspects of conventional FEL designs or when considering novel concepts.The main goal of this paper is to present a systematic, fully 3D analytical method suitable for studying such generalized FEL configurations, focusing on those with varying electron beam parameters.
Following the standard perturbative approach, we derive the coupled, frequency-domain, linearized Vlasov-Maxwell equations for the amplitudes of the radiation field and the perturbation to the distribution function of the electron beam [7].These relations, along with the equation governing the evolution of the unperturbed distribution f 0 , form a closed set which accurately describes the operation of a high-gain FEL in the linear regime, up to the onset of saturation effects.In general, existing FEL theory mostly restricts itself to cases in which the background distribution f 0 is z independent, a typical example being that of a beam that is matched to a focusing channel with a constant transverse gradient.For such systems, it is possible to represent the general solution of the initial value problem in terms of the self-similar, guided eigenmodes of the FEL [8].Each of these modes is characterized by a constant growth rate and a fixed (z invariant) mode profile, which in turn satisfies a properly formulated eigenmode equation.Additionally, if the variation of the electron beam parameters is sufficiently slow, one can construct a generalization of the standard guided mode theory and describe the system in terms of a set of approximate FEL eigenmodes with adiabatically varying growth rates and mode profiles [9].This approach, though insightful and effective, is not the one we choose to follow in this paper.Instead, we endeavor to directly solve the initial value problem without making any assumption about the degree of adiabaticity of the system.
As an example, we consider an FEL that is driven by an initially unmodulated electron beam with varying transverse sizes.This model can represent a mismatched beam traveling through a constant-focusing undulator or an unfocused beam.After some manipulation of the linearized FEL equations, we obtain a single, self-consistent equation for the radiation field amplitude.Our strategy is to expand this quantity in terms of a complete set of orthogonal transverse modes.In general, the latter contain a number of free parameters (usually one or two) that may be complex valued and z dependent.The actual type of basis modes that is to be used depends on the symmetry of the problem at hand.In our case, it is assumed that the system has axial symmetry (round electron/radiation beams and symmetric focusing) so a Gauss-Laguerre expansion is employed.Furthermore, by specifying the single free parameter, one can choose between conventional vacuum modes and waveguide-like guided modes [10] or set up a self-consistent scheme in which the basis parameter evolves in correlation with the complex coefficients of the orthogonal expansion [11].The latter option is typically the preferred one because it leads to highly accurate results with a minimum number of expansion modes.Projecting the amplitude equation onto a suitable basis of our choice leads to an infinite set of coupled evolution equations for the expansion coefficients.The numerical solution of a properly truncated version of this set yields accurate estimates of important quantities such as the FEL gain and the radiation beam size as functions of z.Another notable advantage of such an expansion method is that it is equally robust and efficient for both the initial, transient stage of the interaction and the subsequent exponential gain region.This is not the case for an eigenmode formalism (exact or approximate), which is typically much easier to use in the high-gain regime, where a single FEL mode usually suffices to approximate the radiation field.
This paper is organized as follows: Section II covers the main theoretical development, including the derivation of the equation for the amplitude of the radiation field and the formulation of the expansion technique.In Sec.III, we present numerical results obtained for two different sets of x-ray FEL parameters.Several topics have been explored, including a demonstration of gain guiding, seeding with a non-Gaussian field, mismatch effects and the influence of the absence of focusing upon the FEL performance.The data are shown to be in good agreement both with simulation and with the analytical prediction for the special case of a beam with constant transverse size.The results and conclusions are summarized in Sec.IV.

II. THEORY A. Equation for the amplitude of the radiation field
To begin with, let us assume that the FEL radiation is generated as a relativistic, bunched electron beam moves along the z direction through a planar, parabolic-pole-face undulator.The vector potential A u of the undulator field is given by where k u ¼ 2= u is the undulator wave number, " , and B 0 is the peak on-axis field.In our case, it is assumed that the undulator poles are shaped in such a way that an equal degree of natural focusing is available in both transverse directions.Whenever necessary, additional external focusing is to be provided by thin-lens quadrupoles inserted between the undulator segments in a FODO pattern (i.e., a sequence of equally spaced quadrupoles with focusing strengths of equal magnitude and alternating sign).Moreover, the electric field E r of the radiation is taken to be linearly polarized in the x direction, i.e.E r ¼ E r x, where E r can be expressed as dE ðx; zÞe ik r ðzÀctÞ : Here, x ¼ ðx; yÞ is the transverse position, c is the speed of light, E is the Fourier amplitude of the radiation field, ¼ !=! r is the scaled frequency, and k r ¼ !r =c is the resonant wave number.The latter quantity is defined by , where 0 ) 1 is the average electron Lorentz factor and K ¼ eB 0 =ðk u m 0 cÞ is the familiar undulator parameter-e and m 0 are the electron charge and mass.It is assumed that E is a slowly varying function of z, i.e., that it changes little over an undulator period u .It is useful to point out that E r being real implies that . Since the emission of radiation occurs primarily near the resonant frequency !r , we anticipate that E will be appreciably different from zero only within a very narrow region around ¼ AE1 (excluding higher harmonics).We also note that this analysis does not include space charge effects.
The interaction between the electron beam and the radiation it emits in the presence of the undulator can be studied in a rigorous and self-consistent manner through the Vlasov-Maxwell equations for the combined system.When the random fluctuations arising from the discreteness of the electrons can be disregarded, the beam can be described in terms of a smooth, ensemble-averaged distribution function f ¼ fð; ; x; p; zÞ.Here, ¼ðk u þk r ÞzÀ !r tþ½k r K 2 =ð8k u 2 0 Þsinð2k u zÞ is the electron phase, ¼ ð À 0 Þ= 0 is the relative energy deviation, and p ¼ dx=dz.For the startup and the exponential gain regimes of the FEL process, a perturbation formalism is applicable, according to which the distribution function can be expressed as f ¼ f 0 þ f 1 , where f 0 is the background distribution of the beam and f 1 represents the perturbation due to the coherent modulations.We will only concern ourselves with cases where f 0 is independent, which is a good approximation for a sufficiently long electron bunch.In the context of this approach, the dynamics of the FEL interaction is governed by the following set of linearized, frequency-domain Vlasov-Maxwell equations (for > 0): while k is the combined focusing strength of the undulator.The latter quantity takes into account both the undulator natural focusing and the contribution of the external focusing lattice, which we assume can be treated in the smooth approximation.Furthermore, the unperturbed distribution f 0 evolves according to the zeroth-order Vlasov equation, and its normalization is given by where l b and N b are, respectively, the length of the bunch and the number of electrons it contains.A detailed derivation of Eqs. ( 3)-( 6) from first principles is given in Appendices A and B. We now introduce the interim variables, where z 0 z À z e and z e is a constant offset.From a particle dynamics point of view, x À and p À represent the transverse position of an electron and the slope of its trajectory at z ¼ z e when their respective values at a given z are x and p (we note that . In terms of these new variables, Eqs. ( 5) and (3) become @ f0 @z ¼ 0 and where a tilde over a quantity denotes that it has been expressed in terms of x À and p À instead of x and p. Equation ( 9) implies that f0 does not have an explicit dependence on z, so f0 ¼ gð; x À ; p À Þ or, equivalently, , where g is an arbitrary function.In particular, we choose an unperturbed distribution that corresponds to a mismatched beam with a round, Gaussian profile: or, equivalently, where is the rms energy spread and À¼ 02 =ð 2 k 2 ÞÀ1 is the mismatch parameter, defined in terms of and 0 , which are-respectively-the rms size and angular divergence of the electron beam when sinðk z 0 Þ ¼ 0. Since the transverse distribution is upright at the aforementioned z locations, the emittance is still equal to 0 .The electron beam size e ðzÞ is given by The special case with À ¼ 0-for which = 0 ¼ m 1=k -corresponds to a Gaussian beam that is matched to the focusing system and has a constant transverse size.The general solution of Eq. ( 10) is e ÀiÁk u e i 0 ðÀzÞ ; where 0 À z e and 0 In this paper, we only consider FEL configurations for which the electron beam is initially unmodulated, so f ðz ¼ 0Þ ¼ 0. Switching back to the original transverse variables yields where À z and " The use of x À and p À as intermediate variables is an example of the method of integration along the unperturbed trajectories [12].Substituting Eq. ( 14) into Eq.( 4) and performing the integration over , we obtain a single integrodifferential equation for the radiation field amplitude E : The integral kernel Ç is given by GENERAL METHOD FOR ANALYZING THREE-. . .Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-3 where is the Pierce parameter [13], expressed in terms of the peak current I ¼ eN b c=l b and the Alfven current We also point out that, in the process of deriving Eqs. ( 15) and ( 16), we have used the fact that % 1 in order to simplify those terms that contain the scaled frequency.The sole exception is the detuning term, which cannot be eliminated since k u z ) 1 for the better part of the interaction.Changing the momentum integration variable from p to " x, we obtain the driven paraxial wave equation, where (recall that z 0 ¼ z À z e , 0 ¼ À z e , and

B. Orthogonal expansion method
Our goal is to obtain a solution to Eq. ( 18) that is consistent with a specified initial amplitude E ðx; 0Þ.The method we adopt is based on expanding the field amplitude in terms of a complete set of orthogonal basis functions.The transverse basis we employ consists of generalized Gauss-Laguerre (GL pq ) modes where ðr; Þ are polar coordinates in the transverse plane, p and q are integers with p ! 0, and L jqj p are the associated Laguerre polynomials.Here, is the fundamental basis mode, defined by means of a complex-valued function and are, respectively, the spot size and the Gouy phase associated with it.From a physical point of view, 1 ðzÞ-which is always assumed to be positive-represents the local Rayleigh length of the mode while À 2 ðzÞ corresponds to the local waist position.The basis elements described above satisfy the orthonormality condition, and reduce to the standard vacuum modes of paraxial optics when is a constant.On the other hand, when the complex beam parameter z À i ¼ z þ 2 À i 1 is chosen as constant, one obtains a set of modes with fixed transverse profiles, analogous to the guided modes encountered in waveguide theory.For the moment, we will assume that ðzÞ is an arbitrary complex function.The expansion for E is C pq ðzÞc pq ðx; zÞ; where the coefficients C pq are chosen to be dimensionless and " " is a constant that can be related to the initial amplitude [see Eq. ( 49)].Inserting the first equality of Eq. (24) into Eq.( 18) yields the operator relation, where Using Eqs. ( 20)-( 23), along with properties of the associated Laguerre polynomials, we obtain the relation which, in conjunction with orthogonality, yields Furthermore, after a lengthy calculation, we obtain with and In the relations given above, 2 F 1 is a Gaussian hypergeometric function, is shorthand for ðÞ and 1 1 ðÞ ¼ Re½ .As we shall see later on, we will mostly treat the basis parameter not as a variable to be specified a priori but rather as an additional unknown quantity to be determined during the solution of the initial value problem.Anticipating this, we have made explicit the dependence of Ã nm pm on .Moreover, and GENERAL METHOD FOR ANALYZING THREE-. . .Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-5 An interesting limiting case arises when k !0, or-more precisely-when 0 =ðk Þ ) 1.The mismatch parameter is then approximated by À % 02 =ð 2 k 2 Þ and the electron beam size is expressed by e ðzÞ ¼ ð 2 þ 02 z 2 0 Þ 1=2 , a result which corresponds to an unfocused, coasting beam with a single waist at z ¼ z e .In this limit, Eqs. ( 30), ( 34), (37), and (38) still hold but Eqs. ( 31)-( 33), (35), and (36) are, respectively, replaced by and Combining Eqs. ( 25) and ( 27)-( 29), we obtain a set of coupled integrodifferential equations for the expansion coefficients: where n ¼ 0; 1; 2; . . .and m ¼ 0; AE1; AE2; . . . .The above equations describe the FEL interaction in terms of the evolution of the expansion coefficients C nm ðzÞ along the undulator.An important preliminary observation is that coefficients with different azimuthal indices m are always uncoupled.This reflects the fact that there exist separable solutions of the form g m ðr; zÞe im for the amplitude of the radiation field, a property which is due to the axial symmetry of the system under study.As a result, we can decompose the general initial value problem into an infinite number of independent component problems, each of which is characterized by a particular value of m.In what follows, we will concentrate on such a reduced seeded problem, the solution of which yields a field amplitude that can be expressed as C nm ðzÞc nm ðx; zÞ: We also point out that the function ðzÞ which is used in the definition of the basis elements is as yet unspecified.
Since the full, infinite-dimensional system of Eq. ( 44) is obviously not amenable to analytical or numerical treatment, truncating the expansion for the field amplitude becomes a necessity.In general, if we keep a large number of modes in the right-hand side (rhs) of Eq. ( 45), the error introduced by such a truncation is small regardless of the choice of .In practice, however, the need for fast numerical calculations, coupled with the desire to obtain useful analytical results whenever possible, limit the number of modes that can be retained.Thus, an additional objective emerges: to determine ðzÞ in such a way that, when only a few modes are used, the resulting truncated system of evolution equations captures as much as possible of the physics contained in its full counterpart.This task is twofold.First, we have to specify the initial value ð0Þ.
Given the external seed E ðx; 0Þ $ e im , the initial conditions for C nm are Taking into account the fact that we obtain the convenient normalization condition Generally speaking, it is desirable to choose ð0Þ so that the magnitude of the coefficients C nm ð0Þ falls off as rapidly as possible with increasing n.In this way, it can be expected that the number of modes that are excited at z ¼ 0 will be kept to a minimum.Here, we shall only consider the special case where the input field consists of a finite BAXEVANIS, RUTH, AND HUANG Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-6 number of vacuum Gauss-Laguerre modes with the same azimuthal index m and-constant-beta parameter in .By choosing ð0Þ ¼ in , we ensure that C nm ð0Þ ¼ 0 for all n > M, where M is the maximum radial index of the seed modes.The next step is to determine how evolves with z.
Let us assume that our goal is to construct an approximate solution using the minimum number of modes advisable for this problem, i.e., the first M þ 1 terms in the expansion of the field amplitude (0 n M).In order to facilitate the optimization of the basis function, we temporarily include the next-order mode in our calculations and truncate at n ¼ M þ 1.With this in mind, the ''equation of motion'' for the coefficient of the aforementioned mode becomes it is evident that the evolution of C Mþ1;m is no longer correlated with that of the lower-order coefficients.In addition, the equation satisfied by the highest-order expansion coefficient becomes homogeneous and admits the unique solution C Mþ1;m ðzÞ ¼ 0. In other words, C Mþ1;m vanishes identically when its initial value is equal to zero, as is the case in our example.The evolution equations for the remaining coefficients are where n ¼ 0; 1; 2; . . .; M while s nM ¼ 0 for n ¼ M and 1 for n < M. To recapitulate, starting from an approximation using M þ 2 modes and an unspecified basis function, we were able-through a suitable choice of ðzÞ-to eliminate the contribution of the highest-order mode, thus enhancing the accuracy of the resulting M þ 1-mode approximation scheme.To include additional expansion modes in our calculation, we choose to retain Eq. (50) (whose form is determined by the seed modes) and simply use Eq. ( 51) with M replaced by is the total number of modes.In terms of methodology, the expansion technique discussed in this section can be viewed as a generalization of the source-dependent expansion introduced in Ref. [11].The latter approach was developed using a differently defined and parametrized set of Gauss-Laguerre modes and only considered seeding with a single (Gaussian) mode.Before we move on, it is essential to mention two very useful results about the radiation power and beam size.When the field amplitude has an e im angular dependence, the general expression for the radiation power derived in Appendix C-Eq.(C6)-can be rewritten as where ¼ l b =c and we have utilized Eq. ( 45) and orthonormality.For a monochromatic signal, Eq. ( 52) reduces to where P 0 ¼ ½" 0 =ðk r Þ R 1 0 dj " " j 2 is the input power and I ¼ P 0 j " E j 2 is the intensity of the radiation.The radiation beam size r ðzÞ is defined by where we have made use of the fact that the intensity profile is radially symmetric.Taking into account Eqs. ( 20), ( 22), (53), and (54), we ultimately obtain

A. Demonstration of gain guiding
To illustrate the expansion method, we have used it to study various aspects of the operation of a high-gain FEL in the linear regime.We have considered two different FEL parameter sets, both of which correspond to hard x-ray machines (Table I).Set 1 roughly describes the design parameters of the LCLS while Set 2 refers to a machine with lower beam energy and more ambitious undulator specifications [14,15].To begin with, we seek to make a connection with standard theory [7] by considering an FEL with the parameters of Set 1, driven by a matched electron beam with a constant size m ¼ 23:14 m.For such a configuration, the amplitude of the radiation field can be expressed as a superposition of the guided FEL eigenmodes, i.e., where nm and R nm ðrÞ are, respectively, the complex growth rate and the radial profile of the E nm FEL eigenmode (c nm are constant coupling coefficients and m is the azimuthal index).Eigenmodes with different m are orthogonal, though the same is not true of modes with the same azimuthal but different radial indices.The fundamental (E 00 ) eigenmode of the FEL is characterized by a Gaussian-like profile and usually has the growth rate with the largest imaginary part (in absolute value).Thus, if c 00 Þ 0, it dominates all the other modes in the high-gain regime, which allows us to approximate the radiation amplitude by E ðx; zÞ % c 00 e i 00 z R 00 ðrÞ; (57) for z ) L 3D , where L 3D ¼ À1=ð2Im½ 00 Þ is the 3D power gain length.Here, we would like to verify this well-known theoretical result (namely, the emergence of a single, Gaussian-like guided mode in the high-gain region of the linear regime) by directly solving the initial value problem with our expansion technique.Assuming a Gaussian external seed for the FEL, we first attempt to obtain a solution employing only a single (Gaussian) expansion mode, so that E ðx; zÞ / C 00 ðzÞc 00 ðx; zÞ.Recalling Eq. ( 21), we note that the GL 00 expansion mode can be conveniently rewritten in terms of the radiation complex beam parameter q r ¼ z À i q 2 À iq 1 (q To further facilitate comparison between our results and those of the eigenmode theory, we can also define a z-dependent, local complex growth rate ðzÞ through the relation C 00 ðzÞ ¼ exp½i R z 0 ð ẑÞd ẑ or, equivalently, For calculations involving only one expansion mode, the (negative) imaginary part of the above quantity is equal to half the local power growth rate ð1=PÞdP=dz ¼ dG=dz, where G ¼ log½P=P 0 is the FEL gain.To improve the accuracy of our single-mode approximation, we want the transverse profile of the input radiation field to match as much as possible that of the fundamental FEL mode.To achieve this, some quantitative knowledge of the properties of the latter is required.A very useful technique for calculating the FEL growth rates and mode parameters is the variational method [6,[16][17][18].Specifically, to obtain the variational data that are included in this paper, we have relied on the algorithms described in Refs.[6] (for the E 00 mode) and [17] (for higher-order eigenmodes).In the case of the fundamental FEL mode, a trial solution of the form expðÀAr 2 Þ is assumed for the radial profile, yielding an highly accurate approximation to 00 and a fairly reliable estimate of the mode shape.For our discussion, it is more meaningful to parametrize the trial solution in terms of the betalike quantity 0 ¼ iq 0 ¼ k r =ð2AÞ, where q 0 is the complex beam parameter of the mode.Using the variational method, one can also optimize the fundamental growth rate, i.e., maximize its imaginary part with respect to the detuning.For the LCLS parameters, the variational approach yields an optimum normalized growth rate of 00 ¼ 00 =ð2 m k u Þ ¼ 0:247 À i0:717 and a corresponding mode parameter given by 0 = m ¼ 0:384 þ i0:206 (i.e., a mode size 57.5% the size of the electron beam) for a scaled detuning ¼ Á=ð2 m Þ ¼ À0:38, which in turn corresponds to a wavelength of 1.500 62 A ˚. We select the latter as the seed wavelength and also choose ð0Þ ¼ 0 [the remaining initial condition is simply C 00 ð0Þ ¼ 1].Using the approximation scheme defined by Eqs. ( 50) and (51) for M ¼ 0, we obtain the single-mode results included in Figs.1-4.In particular, in Fig. 1-which shows the evolution of the FEL gain-we observe that an initial lethargy stage is followed by a rapid power buildup, eventually leading to a region in which the gain increases linearly with z.Moreover, from Fig. 2, it is evident that the local complex growth rate asymptotically approaches a constant value in the high-gain limit.This is also the case    for the complex beam parameter q r ¼ z þ 2 À i 1 , as we can see in Fig. 3.The last two results confirm that the radiation field eventually evolves into a guided mode with a constant size.Additional verification of this gain guiding effect is provided by Fig. 4, where the radiation beam size, after an initial transient, is also shown to approach a constant value in the high-gain limit.The obtained asymptotic values for (in units of 2 m k u ), q r (in units of m ), and r (in units of m ) are, respectively, 0:250 À i0:717, 0:196 À i0:425, and 58.7%.We note that these results are in good agreement with the variational prediction, especially in the case of the growth rate (which is, in any case, the quantity most accurately calculated by the variational method).A systematic analysis of the asymptotic solutions, including their generalization in case the size of the electron beam is not constant, requires the study of dispersion relations similar to those obtained in the context of the eigenmode formalism.We will defer a detailed discussion of this important topic to a subsequent publication.As a closing remark, we note that Figs. 1 and 4 also include data obtained using multiple (up to five) expansion modes.In fact, Fig. 4 shows the gradual convergence of the radiation size estimates as the number of modes increases.
For the radiation power PðzÞ ¼ P 0 e GðzÞ and the radiation beam size, our calculations show that the 1-mode results differ from their 5-mode counterparts by about 3% and 1%, respectively (in an rms sense).This suggests the potential usefulness of the single-mode approximation for deriving reliable results with a minimum of computational effort.

B. Higher-order mode (HOM) seeding
Next, we would like to explore the degree to which seeding with a non-Gaussian input field affects the FEL performance.We focus on the parameters of Set 2 and again assume that the electron beam is matched to the undulator, which only relies on natural focusing.However, we now consider two distinct cases: seeding with a Gaussian and a GL 01 mode.For both cases, the input parameters are chosen to be ð0Þ= m ¼ 0:92 þ i0:91 and ¼ Á=ð2 m Þ ¼ À0:94 (these parameters roughly match the Gaussian seed to the optimized fundamental FEL eigenmode).Figures 5 and 6 summarize the resulting single-mode data for the two seeding schemes.From Fig. 5, it follows that, after traversing an undulator length L u ¼ 34:45 m, the total FEL gain for the GL 01 seed is less than what it would be for the GL 00 seed by an amount ÁG 01 % À3.In other words, if the input power is the same in both cases, the amplified power ratio is P 01 =P 00 ¼ e ÁG 01 % 5%, a result which clearly reflects the less than optimum overlap of the higher-order mode with the electron beam.Moreover, the two gain curves have constant but different slopes in the high-gain region.This suggests that different eigenmodes prevail in each  FIG. 5. Gain functions for an FEL driven by a matched beam (Set 2 parameters) and seeded with a GL 0m mode: data for m ¼ 0 (blue) and m ¼ 1 (red) are shown.Figure 6 also refers to this configuration.Both 1-mode (solid curves) and 5-mode results (symbols) are included.In terms of the radiation power, the rms deviation of the latter from the former is about 1% (2%) for the GL 00 (GL 01 ) seed.

BAXEVANIS, RUTH, AND HUANG
Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-10 case.In Fig. 6, we compare the z-dependent, local growth rates with the growth rates of the dominant FEL modes.In the case of the GL 01 seed, we need to point out the following: (a) The local growth rate is now defined by ðzÞ ¼ À½i=C 10 ðzÞdC 10 ðzÞ=dz.(b) Of the FEL modes with nonzero coupling to the input field (i.e.those with m ¼ 1), the dominant one is the E 01 eigenmode, which has a transverse profile similar to that of the seed.The variational method yields 00 ¼ 0:259 À i0:514 and 00 ¼ 0:227 À i0:438 while the corresponding asymptotic values of the local growth rates are 0:261 À i0:514 and 0:232 À i0:438 (in units of 2 m k u ).Good agreement is again observed between the two quantities.
Another aspect of this analysis involves seeding with a GL n0 mode, where n Þ 0. In particular, we assume that the previously studied FEL configuration (matched beam/Set 2 parameters) is now seeded with a GL 10 mode, for which ð0Þ= m ¼ 0:230 þ i0:227 and ¼ À0:94.In this example, we use at least two expansion modes to construct a solution so we set M ¼ 1 in Eqs. ( 50) and (51).The numerical results are presented in Figs.7-9, which also include data for the case of the Gaussian seed.From Fig. 7, we deduce that the difference in total gain (in favor of the GL 00 seed) is ÁG 10 % À0:9, which yields a power ratio P 10 =P 00 ¼ e ÁG 10 % 40% at the exit of the undulator.Unlike the previous case, the two gain curves appear to have approximately the same constant slope in the highgain regime, the only difference being in the extent of the lethargy region.This is verified by an inspection of Fig. 8, in which the power growth rate ð1=PÞdP=dz (in units of 4 m k u ) is plotted as a function of z.This property is due to the fact that both the GL 00 and the GL 10 seed modes have nonzero coupling to the E 00 eigenmode so the latter is the dominant FEL mode in both cases.Finally, in Fig. 9 we plot the ratios (60) (for n ¼ 0; 1) versus z.Bearing in mind that (in the context of the 2-mode approximation) the amplified power is given by P ¼ P 0 ½jC 00 j 2 þ jC 10 j 2 , it becomes clear that these quantities represent the square root of the fraction of the power due to each expansion mode.The transition of the system from an initial state defined by the GL 10 seed to an eventual state dominated by the Gaussian-like E 00 FEL eigenmode is evident.

C. Mismatched/unfocused electron beam
In all the cases studied so far, the size of the electron beam was constant.The remaining part of the numerical  FIG. 9. Amplitude ratios r 00 and r 10 -defined by Eq. ( 60)versus z.The plots were generated using the 2-mode results for the case of the GL 10 seed.results deals primarily with z-dependent schemes.Let us consider again the example studied earlier in the context of demonstrating gain guiding (matched beam/Gaussian seed/ LCLS parameters).We now expand this example by including (for the same undulator and input field) a case for which the electron beam is underfocused, with z e ¼ 0, = m ¼ ffiffiffiffiffiffi ffi 2:5 p , and À ¼ 4 m = 4 À 1 ¼ À0:84.In Figs. 10 and 11, we show the comparison between the results obtained with our technique and simulation data given by GENESIS [19] (for the latter, an input power P 0 ¼ 2kW has been assumed).In the linear regime, the theoretical results are in good agreement with simulation, even though only a single Gaussian mode has been used in obtaining the former.Including higher-order modes in our calculation makes the comparison with simulation even more favorable.As we can see in Fig. 10(b), the power growth rate for the mismatched beam lags behind that for the matched beam in the earlier stages of the interaction but this situation is eventually reversed.As a result, the total gain is almost the same for the two cases [see Fig. 10(a)-interestingly, this observation appears to be valid both for the gain derived from the linearized solution and for the actual saturation gain given by GENESIS].This is due to the fact that the oscillation of the electron beam size in the case of the underfocused beam [Fig.11(b)] allows that variable to attain values smaller than the matched size of 23:14 m and closer to the optimum beam size opt % 0:47 m (which corresponds to a beta function of 6.72 m).This optimized value can be determined by minimizing the 3D gain length with respect to the matched beta using Ming Xie's fitting formula [6], though a value larger than the optimum is usually chosen.In our case, the electron beam size is minimum and equal to min ¼ ð1 þ ÀÞ 1=2 % 0:63 m when z ¼ m =2 % 47 m.In the neighborhood of this minimum, the electron beam size is approximately constant and as close as possible to the optimum size, which is why the power growth rate becomes maximum in the same region.From Fig. 11(a), we also observe that the radiation size roughly follows the evolution of the electron beam size in  the high-gain regime, as has been noted in the literature (Ref.[20]).
To conclude this section, we switch back to the parameters of Set 2 and consider an FEL based on an undulator with the listed values of K and u but with virtually no focusing (natural or external).This option can refer to a planar rf undulator (Ref.[21]), where the transverse defocusing effect is typically very weak.As we have mentioned before, this case describes an FEL that is driven by an electron beam with a waist at z ¼ z e .We choose the waist beta function e = 0 to be numerically equal to the natural focusing value of m ¼ 13:78 m, which is fairly close to the optimum beta for these parameters (10.6 m).Thus, we expect the neighborhood of the waist to be the optimum region for such a beam, as far as the FEL gain is concerned.In Figs. 12 and 13, we compare 1-mode (solid curves) with 5-mode results (markers) for three beams with z e ¼ 0=0:5L u =L u , where L u ¼ 2:5 e ¼ 34:45 m is the undulator length.For these, as well as all remaining calculations, we have assumed the same Gaussian seed that was used for deriving part of the data for Figs.5-9.Overall, the rms deviation between the 1-mode and the 5-mode estimates for the amplified power and the radiation size is of the order of 1% or less.From Fig. 12(a), we note that the highest total gain is registered by the configuration with z e ¼ L u =2, for which the overlap of the optimum region of the beam with the undulator region appears to be best [Fig.13(b)].Moreover, an inspection of Fig. 12(b) verifies that the position at which the power growth rate becomes maximum is roughly identical to the waist position of each beam [though this is less clear for the blue curve (z e ¼ 0), for which the optimum region around z % 0 lies mostly in the transient regime].The corresponding maximum growth rate (in the high-gain region) seems to be independent of the position of the waist.Figure 13(a) again establishes that the evolution of the size of the radiation beam is closely correlated with that of the electron beam size in the exponential gain regime.Finally, having shown that the 1-mode formalism is reliable in this case, we can use it to study in more detail the influence of the waist position upon the total FEL gain (for the same e and L u as before).
The results are summarized in Fig. 14, where we also include the total gain for a beam that is matched to the conventional version of the undulator.From this parameter scan, we conclude that one can indeed recover as much as 93% of the gain for the matched beam when z e ¼ 1:25 e ¼ 0:5L u , i.e., by placing the waist of the unfocused beam in the middle of the undulator segment.

IV. CONCLUSIONS
An expansion method has been developed for solving the initial value problem of an FEL with variable electron beam parameters, taking full account of 3D effects.We have used this technique in our study of a high-gain FEL that is driven by a beam with varying transverse sizes.In particular, we have considered a round, Gaussian electron beam with no initial modulation that, in general, is not properly matched to a constant-focusing undulator system.This model includes the constant-size beam as a special case and can also deal with the case of an unfocused, coasting beam by taking the limit of vanishing undulator focusing strength.Starting from the linearized Vlasov-Maxwell equations for the FEL, we obtained a driven paraxial wave equation for the slowly varying amplitude of the radiation field.By expanding the radiation amplitude in terms of a set of Gauss-Laguerre transverse modes, we eventually derived a truncated set of coupled, integrodifferential equations for the mode coefficients and the variable basis parameter.The latter is selected so as to approximately match the retained basis modes to the evolving radiation profile.Given the input field, a straightforward numerical solution of this set leads to accurate estimates of the radiation power and beam size.Several topics were covered during this numerical study.For a matched beam, these included a direct demonstration of the gain guiding effect and seeding with a higher-order mode (for which less gain was registered-compared to the case of a Gaussian seed-due to the less than optimum coupling to the electron beam).For the cases with an explicit z dependence, we explored the effect of betatron mismatch on the FEL output and also showed how the parameters of an unfocused beam can be optimized so as to yield the maximum possible FEL gain.Overall, the results obtained are in good agreement both with GENESIS simulations and with the predictions of the standard guided mode theory (the latter, of course, refers to the special case of a matched electron beam).The ability to describe the radiation amplitude with a small number of expansion modes makes this method potentially useful for parameter studies.Moreover, this approach can be thought of as a possible alternative (or benchmark) to simulation when one seeks to solve the initial value problem in the linear regime in order to validate other theoretical techniques, such as the adiabatic version of the eigenmode formalism presented in Ref. [9].Finally, we would like to stress that the broad applicability of this method enables us to use it in the analysis of a wide range of FEL configurations and parameter regimes.

APPENDIX A: SINGLE PARTICLE EQUATIONS OF MOTION
The purpose of this Appendix is to derive the equations of motion for a single electron of charge Àe and mass m 0 in the combined field of the undulator system and the radiation.In a gauge where the scalar potential is equal to zero, the radiation field is derivable from a vector potential A r ¼ A r x ¼ À R E r dt.Thus, the total vector zÞ represents the longitudinal vector potential from which the external focusing fields are derived.For the linear focusing channel under consideration, A z is given by where the focusing function kðzÞ is periodic with period equal to the length L of the basic FODO cell.The relativistic, single particle Hamiltonian is where E is the electron energy and P i (i ¼ x; y; z) are the canonical momenta.Next, we would like to use the longitudinal position z as the independent variable and also eliminate the time t by introducing the phase c ¼ ðk u þ k r Þz À !r t as a new coordinate.These tasks can easily be accomplished simultaneously by noting that the action integral can-by virtue of the phase definition-be rewritten as BAXEVANIS, RUTH, AND HUANG Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-14 where P c ¼ H 1 =! r ¼ E=! r is the new conjugate momentum and is the new Hamiltonian.Solving Eq. (A2) for P z and substituting into Eq.(A3), we obtain where we have used the small angle approximation, i.e., the fact that the total mechanical momentum . The next step is to introduce the fractional energy deviation, where E 0 ¼ 0 m 0 c 2 is the average energy of the electrons, as the new canonical momentum conjugate to the phase c .
To this purpose, we scale H 2 -along with the remaining canonical momenta-by the appropriate factor and expand p m in terms of .Up to second order in , this manipulation yields where P x;y ¼ !r P x;y =E 0 , Ãx ¼ Ãux þ Ãr , and Ãj ¼ ð! r e=E 0 ÞA j (j ¼ ux; uy; z; r).Assuming that the radiation vector potential is a small quantity-that is, jA r j ( jA u j-the term in the brackets in Eq. (A6) can be written as where we have omitted the quadratic Ã2 r term in the expansion of Ã2 x .For small transverse displacements from the undulator axis, i.e., when k u jxj and k u jyj ( 1, Eq. ( 1) yields In view of the above relations and after some rearranging and truncation, Eq. (A6) becomes Here, the fourth-order term Ã2 uy $ x 2 y 2 has been omitted and k n ¼ Kk u =ð2 0 Þ ¼ K=ð 0 u Þ is the undulator natural focusing strength.As is well known (see, for instance, Ref. [7]), the horizontal position x of an electron can be expressed as x ¼ x w þ x , where x w ¼ À½K=ð 0 k u Þ Â sinðk u zÞ is the term due to the fast wiggling motion and x is the contribution of the slow betatron motion.We note that the expression for x w does not include a small quadratic correction due to the transverse dependence of the undulator field.Furthermore, in a planar undulator, the longitudinal velocity is not constant-as would be the case for a helical geometry-but oscillates about an average value at twice the transverse wiggle frequency ck u .As a result, c can be decomposed into a fast oscillating part c w ¼ À½k r K 2 =ð8k u 2 0 Þ sinð2k u zÞ and a slowly varying phase .To switch to the slowly varying coordinates (, x ), we employ the second-type canonical transformation defined through the generating function, F 2 ¼ F 2 ðc ; x; y; ; P 1 ; P 2 ; zÞ Finally, we wish to retain only those terms in the Hamiltonian that vary slowly on the scale of the undulator period u .To this end, we first pay closer attention to the terms that involve Ãr .Using Eqs. ( 2) and (A10), the latter quantity can be written as Þ and we have also used the identity We also need the expansion E ðx; zÞ ¼ E ðx ; zÞ À @E @x r w sinðk u zÞ þ Á Á Á ; where r w ¼ K=ð 0 k u Þ is the amplitude of the wiggling motion.In what follows, we assume that r w is always much smaller than the radiation/electron beam sizes.In that case, we may keep only the zeroth-order term in the expansion given above.In view of this approximation and bearing in mind that (a) % AE1, and (b) P x , P y , and are slowly varying quantities, it becomes clear that Ãr P x only has fast oscillating components.This is not the case for the Ãux Ãr =k r term, which (ignoring the transverse dependence of Ãux ) can be written as Using the previously quoted identity-Eq.(A18)-and properties of the Bessel functions, we conclude that the slowly varying part of the cross product term is where 1 ¼ eKJJ=ð4 2 0 m 0 c 2 Þ and JJ ¼ J 0 ðQÞ À J 1 ðQÞ while a is defined by a E e ÀiðÀ1Þk u z ; > 0 e Àiðþ1Þk u z ; < 0 (A20) for Þ 0 and is equal to zero for ¼ 0 (note that a À ¼ a Ã ).The averaging of the remaining terms on the rhs of Eq. (A16) proceeds along similar lines.The final, wiggleaveraged Hamiltonian is given by In one important special case, the above Hamiltonian can be further simplified.If k n L ( 1, the undulator segments can be treated as drift spaces as far as determining the betatron functions is concerned.From the standard results for a FODO cell, it is known that the lattice beta functions x and y oscillate between a maximum and a minimum value given by AE ¼ L½1 AE sinðÈ=2Þ= sinÈ, where È is the phase advance per cell.The latter quantity is in turn expressed by the relation sinðÈ=2Þ ¼ L=ð4f q Þ, where f q is the magnitude of the quadrupole focal lengths.Furthermore, the alpha functions a x and a y are zero at the quads and equal to AE1= cosðÈ=2Þ at the midpoints of the drift spaces.In the limit of small phase advance (È ( 1), the beta functions are approximately constant and equal to an average value " % L=È % 2f q while the evolution of the alpha functions can be approximated by a periodic step function with values equal to AE1.To rigorously derive the simplified Hamiltonian in the framework of this so-called constant focusing or smooth approximation, we perform the canonical transformation given by which can be obtained from the second-type generating function, " F 2 ¼ " F 2 ð; x ; y ; " ; Px ; Py ; zÞ The new Hamiltonian H ¼ H þ @ " F 2 =@z is given by BAXEVANIS, RUTH, AND HUANG Phys.Rev. ST Accel.Beams 16, 010705 (2013) 010705-16 where a prime denotes differentiation with respect to z.The transverse conjugate variables can be related to the actionangle variables ðJ x ; J y ; c x ; c y Þ of the linear betatron motion via Given that (a) the actions J x , J y are exact invariants of the motion in the absence of radiation, and (b) the beta functions x , y are almost constant while c x ðz þ LÞ À c x ðzÞ ¼ c y ðz þ LÞ À c y ðzÞ ¼ È ( 1, we conclude that x ¼ ðx ; y Þ and P ¼ ð Px ; Py Þ change little over the length L of a FODO cell.Moreover, a x , a y and the focusing function kðzÞ average to zero over that same length.Using the above observations and the periodicity of the betatron functions, it becomes evident that the last three terms in the Hamiltonian do not contribute to an averaging over the cell length.Taking into account the fact that x ; y % " and a 2 x ; a 2 y % 1, we obtain the final form of the averaged Hamiltonian: is the total focusing strength.Introducing the new momentum variable p ¼ P=k r , taking into account the fact that jj ( 1, and neglecting those terms which describe the influence of radiation on the betatron motion, the single particle equations of motion become Equations (A27) and (A28) describe the betatron oscillations in the smooth approximation while Eqs.(A29) and (A30) constitute the 3D generalization of the 1D FEL pendulum equations, including the emittance effect.

APPENDIX B: VLASOV-MAXWELL FORMALISM
On a microscopic level, the phase space density of the electron beam is given by the Klimontovich distribution, where N b is the total number of electrons in the bunch.In cases where random effects such as shot noise can be neglected, we may instead deal with a smoothed distribution function f ¼ hf K i, where the brackets denote an ensemble averaging.The evolution of f along the undulator is governed by the Vlasov equation, while its normalization is given by where z b ¼ =k r is the longitudinal position with respect to the centroid of the bunch.In the startup and the exponential gain regimes of the FEL interaction, one can employ perturbation theory in order to obtain a linearized version of the Vlasov equation.We set f ¼ f 0 þ f 1 , where f 0 is the background distribution and f 1 is the perturbation caused by the microbunching effect (jf 1 j ( f 0 ).Here, we will only consider cases with @f 0 =@ ¼ 0, which corresponds to a bunch with a uniform longitudinal profile.This approximation is valid when the bunch length l b is much larger than the slippage length N u r , where N u is the number of undulator periods.Taking into account Eqs.(A27)-(A30) and the fact that the radiation amplitude is to be treated as a first-order quantity, we find that the unperturbed distribution f 0 ð; x ; p; zÞ satisfies the equation and is normalized according to df 0 ð; x ; p; zÞ ¼ N b (B5) while the linearized Vlasov equation for the perturbation f 1 ð; ; x ; p; zÞ is da ðx ; zÞe i ; (B6) where 0 ¼ d=dz is given by Eq. (A29).Making use of the Fourier representation f 1 ¼ R 1 À1 df e i -where f ¼ R 1 À1 df 1 e Ài =ð2Þ)-Eq.(B6) is carried entirely into the frequency domain, since it can be written as The above result is valid for all .As expected, different frequencies remain uncoupled in the linear regime.Furthermore, we recall that E r and f 1 are real quantities, which implies that E À ¼ E Ã and f À ¼ f Ã .This property allows us to restrict the remainder of our analysis to the positive frequency domain, in which the Vlasov equation has the form @f @z þ p @f @x À k 2 x @f @p þ i 0 f ¼ À 1 E ðx ; zÞe ÀiÁk u z @f 0 @ ; (B8) where we have introduced the detuning Á ¼ À 1.
We now turn our attention to the electrodynamics aspect of our derivation.The radiation E r satisfies the wave equation, where e (j x ) is the charge (current) density.Substituting Eq. ( 2) in the rhs of the relation given above yields where we have dropped the @ 2 E =@z 2 term in the parentheses as a consequence of the slowly varying approximation (SVA) for the field amplitude E .Inverting the Fourier transform, we obtain dt @ e @x þ 1 c 2 @j x @t e Àik r ðzÀctÞ ¼ e iÁk u z " 0 X 1 l¼À1 J l ðQÞe ið1þ2lÞk u z Â Z 1

À1
d @ e @x À i k r c j x e Ài : (B11) Here, we have performed an integration by parts in order to transform the @j x =@t term and have also utilized the definition of and Eq.(A18).Next, we calculate the charge and current densities in terms of the distribution function f.This can be done in a systematic way by first considering their microscopic counterparts.We have m e ¼ Àen m , where the (microscopic) electron number density n m is given by Both multiple integrals on the rhs of Eq. (B13) are slowly varying functions of z as the SVA is equally valid for the distribution function f.After Eq. (B13) is substituted into Eq.(B11) and the familiar averaging over an undulator period is performed, we obtain @ @z þ r Once more, we substitute f ¼ f 0 þ f 1 (with @f 0 =@ ¼ 0) and recall that E is a first-order quantity.The nontrivial result is a driven paraxial wave equation [Eq.( 4) in the main text] @ @z þ r where 2 ¼ eKJJ=ð2" 0 0 Þ.

APPENDIX C: RADIATION POWER
In this section, we derive an analytical expression for the radiation power, working from first principles.The magnetic field of the radiation is expressed by B r ¼ r Â A r ¼ r Â A r x ¼ rA r Â x.The corresponding Poynting vector is where 0 is the vacuum permeability, and its longitudinal component is given by S z ¼ ðE r = 0 Þð@A r =@zÞ.From Eq. (2), we obtain As a result of the paraxial approximation, we can disregard the first term in the parenthesis.Thus, we obtain @A r @z % 1 2c dd 0 E E 0 e iðþ 0 Þk r z e Àiðþ 0 Þk r ct : (C4) The total radiation energy from a single electron bunch is dd 0 E E 0 e iðþ 0 Þk r z Z 1

À1
dte Àiðþ 0 Þk r ct ¼ " 0 k r (C5) In the last step, we have utilized the integral representation of the delta function as well as the relation E À ¼ E Ã .Taking into account the bunched structure of the electron beam, the average radiation power is given by where ¼ l b =c is the bunch duration.

FIG. 3 .
FIG. 3. Single-mode data for (a) the local Rayleigh length and (b) the real part of the radiation complex beam parameter in units of m (blue curves).The dashed lines are the corresponding variational values for the E 00 eigenmode.

FIG. 1 .
FIG.1.Gain function GðzÞ ¼ log½PðzÞ=P 0 for an FEL driven by a matched beam and seeded with a Gaussian field (Set 1 parameters).Figures2-4also refer to this configuration.Both 1mode (blue curve) and 5-mode results (symbols) are shown.

FIG. 4 .
FIG.4.Radiation beam size in units of m , using various numbers of expansion modes.Also shown is the variational estimate of the fundamental mode size (dashed line).

FIG. 6 .
FIG. 6.(a) Real and (b) negative imaginary parts of the scaled, local growth rate 0 ðzÞ ¼ ðzÞ=ð2 m k u Þ for a GL 0m seed mode (solid lines, based on single-mode data) compared to the variational values for the growth rate of the E 0m eigenmode (dashed lines).Results for m ¼ 0 (blue) and m ¼ 1 (red) are shown.

FIG. 8 .
FIG.8.Local power growth rate ð1=PÞdP=dz (in units of 4 m k u ) for a Gaussian (blue) and a GL 10 seed mode (red).The dashed line is the growth rate of the fundamental FEL mode, obtained from the variational method.

FIG. 7 .
FIG.7.Gain functions for an FEL driven by a matched beam (Set 2 parameters) and seeded with a Gaussian (blue solid line) and a GL 10 mode (red solid line).Figures8 and 9also refer to this configuration.These results have been obtained using the minimum number of expansion modes suitable for each case.Also included are 7-mode data for the HOM seed (red markers), which are in good agreement with the 2-mode results.

FIG. 11 .
FIG. 11.(a) Radiation and (b) electron beam sizes in units ofm for a matched and a mismatched electron beam.The legend is the same as in Fig.10.
FIG. 10. (a)FEL gain GðzÞ ¼ log½PðzÞ=P 0 and (b) normalized local power growth rate ð4 m k u Þ À1 dGðzÞ=dz for a matched (blue curves) and a mismatched beam (red curves)-LCLS parameters.Both analytical results (solid curves) and GENESIS simulation data (dashed curves) are shown.

FIG. 13 .
FIG. 13.(a) Radiation and (b) electron beam sizes in units of m for three different unfocused beams.The legend is the same as in Fig. 12.