Variable polarization states in free-electron lasers

Free-electron lasers (FELs) can emit light with different optical polarizations including linear, elliptic and circular polarizations corresponding to the characteristics of the undulators used. X-ray FELs depend upon long undulator lines consisting of a sequence of short undulators. Linearly polarized undulators are most commonly used; hence the optical output is linearly polarized. Alternately, APPLE-II, Delta undulator designs, or a sequence of linearly polarized undulators with alternating orientations can be used to produce undulating magnetic fields with arbitrary polarizations. We present a three-dimensional, time-dependent formulation that self-consistently includes two optical orientations and, therefore, treats any given sequence or combination of undulator including undulator imperfections and degradation There are two principal characteristics of the formulation that underpin this capability. First, particles are tracked using the full Newton–Lorentz force equations with analytic models of the undulators fields. This permits an accurate model of the interaction of the electrons with a large variety of undulator fields and orientations. Second, the electrons can couple simultaneously to two independent electromagnetic polarizations and, therefore, the optical polarization evolves self-consistently along the undulator line. We present the numerical model and give some examples using prevailing undulator configurations.


Introduction
X-ray free-electron lasers (FELs) are an active area of research and development [1][2][3][4][5], resulting in a variety of configurations that generate ever shorter wavelengths. In particular, there is interest in the generation of a variety of different optical polarizations. These FELs use long sequences of relatively short undulators separated by gaps that may contain focusing magnets (typically quadrupoles), phase shifters (dipoles) and various diagnostics. The most common undulator type in use is linearly polarized and this means that the optical output from this configuration is also linearly polarized. However, undulator imperfections and possible undulator degradation may result in a polarization state different from the purely linear. Furthermore, there is interest in generating a variety of different polarizations. The most straightforward way to accomplish this is to use a line of linearly polarized undulators where the orientations of the undulator magnetic fields vary in some fashion along the undulator line. For example, the undulator line might be configured consisting of a first stage with a given undulator polarization followed by a second stage with a different undulator polarization that functions as an afterburner [6][7][8][9]. Since the interaction is close to saturation in the first stage, the electron beam is preconditioned for strong excitation in the afterburner and this can generate elliptically polarized output. A configuration like this has been implemented on the Linac Coherent Light Source (LCLS) using the so-called Delta undulator [10] as an afterburner [11]. Alternately, the APPLE-II [12] undulator design is configurable to produce any type of undulator field from linearly polarized through elliptic polarization to a helically polarized field and this can be used for the complete undulator line [13] or as an afterburner. In view of the growing interest in novel undulator configurations, it is important to develop a numerical formulation that can selfconsistently model the evolution of the polarization of the light through an arbitrary sequence of undulators.
The bulk of existent FEL codes rely on the slowly-varying envelope approximation (SVEA) [14][15][16][17][18][19][20] and use field models with fixed polarizations to describe the interaction with undulators with linear, helical or elliptical polarizations. An exception to this is the particle-in cell code PUFFIN which can treat arbitrary polarizations [21].
Here, we describe a time dependent, three-dimensional nonlinear formulation that is also based on the SVEA and which self-consistently models arbitrary optical polarizations due to varying undulator configurations, imperfections in the magnetic field, possible undulator degradation, or any combination of these effects. We apply this formulation to the study of variable polarizations produced by an APPLE-II undulator line. The full Lorentz force equations are used to propagate the electrons through the optical and magnetostatic fields without applying the wiggler-orbit approximation. The description of the magnetostatic fields includes three-dimensional models for linear, helical, and elliptical undulators and includes the fringing fields associated with entry and exit tapers. Analytic models for quadrupoles and dipoles are also included. This permits the treatment of particle dynamics in the native fields of any arbitrary configuration. The formulation uses a superposition of Gaussian modes to model the electromagnetic field, which typically starts from noise, followed by an exponential growth before reaching a nonlinear post-saturation state.
As such, the present formulation follows that described previously for the magnetostatic fields and particle tracking [17,18]. The difference comes in with the treatment of the electromagnetic field. Instead of prescribing a fixed phasor polarization for the optical field, the new formulation introduces two perpendicular field components (x, y) with the z-direction along the axis of symmetry of the undulator line, that simultaneously interacts with the electrons. The formulation allows the polarization of the light to dynamically evolve upon propagation along the undulator line. Both formulations use the SVEA where the dynamic equations for the optical field are time-averaged to remove the rapid field oscillations such that these equations reduce to equations for the slowly-varying amplitude and phase. In the new formulation there are now twice as many field equations to be integrated as previously. Fortunately, this presents a relatively small additional computation load on modern computers where vectorization and multiple processors can be employed.
The organization of the paper is as follows. The formulation is described in section 2. Examples are given in section 3. A summary and discussion is given in section 4.

The numerical formulation
To summarize, both the old [18,19] and new formulations are time-dependent and use three independent spatial dimensions to describe the particles and fields. As before, we still use MINERVA as the name for the new formulation. Magnetic field elements are used to create a transport line and can be placed in any arbitrary sequence and typically contain single or multiple undulator segments. To create a FODO lattice, quadrupoles can be placed at any location and are allowed to overlap with undulator segments. In addition, for modeling high-gain harmonic generation (HGHG), optical klystron configurations or phase shifters, dipole chicanes can be placed between undulator segments. The magnetic field models and a statement of the dynamical orbit equations have been described previously [18,19].
In the earlier formulation [18,19], the optical field was decomposed into Gauss-Hermite modes when using planar undulators, while Gauss-Laguerre modes were used for helical or elliptical undulators. The present model uses only the Gauss-Hermite modes; however, comparison of helical wiggler simulations with both formulations shows near-identical results using the two Gaussian mode representations; hence, there is no loss of generality using only the Gauss-Hermite modes. As the optical field is included in the Newton-Lorentz force equation, harmonics are included self-consistently.
A fourth-order Runge-Kutta integrator is used for all parts of the transport line, i.e., within the undulator(s) as well as in the gaps, dipoles, and quadrupoles. As such, the optical phase is determined self-consistently for all parts of the transport line. Likewise, the slippage between optical field and electrons is also included selfconsistently along the transport line.
The representation of the Gauss-Hermite modes in present use are given by where (l,n) denote the transverse mode numbers, h is the harmonic number, and dA l n h , , are the mode amplitudes. Here  where j=x,y and H l and H n are the Hermite polynomials, ( ) are the waist and spot sizes of the hth harmonic respectively. The phase is ( ) a h j denotes the curvature of the phase front of the hth harmonic. The mode amplitudes, the spot sizes and the curvature are assumed to be a slowly-varying function of z and t.
The total power in each mode, P l,n,h , is the sum of the powers in each polarization direction over the entire cross section which use the source-dependent expansion (SDE) [22], where ( ) S l n h i j , , , are the source terms and i=(1, 2) and j=(x,y) denoting the x-and y-directions, d/dz is the convective derivative, where ω b is the beam plasma frequency, <(K)> denotes an average over the initial particle distribution [17,18] and υ j ( = j x y , ) are the components of the electron velocity. In the old formulation there was one spot size and one curvature associated with the fixed polarization phasor. In the new formulation, there is a separate spot size and curvature associated with each polarization direction and the SDE equations are We remark that when no electron beam is present ( ) and vacuum diffraction is recovered.

The number of equations in the simulation is
where N slices is the number of slices in the simulation, and N particles is the number of particles in each slice, N modes is the total number of optical modes in all the harmonics, and N harmonics is the number of harmonics. Note, the number of modes is not necessarily the same for each of the harmonics. The Runge-Kutta algorithm allows the step size to change so different steps can be used in each magnetic element or in the drift spaces so this imposes no limitation on the placement of components along the electron beam path.
It should be remarked that the particle trajectories are integrated using the complete, three-dimensional Lorentz force equations using the complete field representations for the electromagnetic and magnetostatic fields. Hence, harmonic interactions are implicitly and self-consistently included in the formulation.
We rely on the Stokes parameters (s 0 , s 1 , s 2 , s 3 ) to characterize the polarization of the optical field [23]. Further, we determine the eccentricity of the polarization ellipse as well as the orientation of the semi-major axis of the ellipse with respect to the x-axis [23] in case of elliptically polarized light, as illustrated in figure 1. If we denote the angle between the semi-major axis of the ellipse and the x-axis as θ, then this angle is determined by The eccentricity, ε, of the ellipse is given by The eccentricity of a purely circularly polarized optical field is zero while that for a purely linearly polarized optical field is unity. Elliptical polarizations are characterized by eccentricities between these extrema. Measurements of the polarization in FELs have been characterized by the Stokes parameters [24]. At the present time, we have implemented a preliminary diagnostic for the Stokes parameters and the eccentricity based upon the on-axis field of the fundamental TEM 00 mode. A more complete diagnostic is under development. Furthermore, since the optical field is composed of multiple temporal slices each of which may behave differently, we provide the on-axis eccentricity ε and angle θ as a power-weighted average over the slices.

Numerical simulations
In this section we consider three different configurations. The first is a comparison with the SPARC experiment which is a planar undulator-based SASE FEL [25] and which serves as validation of the numerical formulation. This experiment has been amply studied with a variety of simulation codes [18,19,25,26]. The fundamental parameters of the SPARC experiment are used for all the subsequent configurations, but with variations. The planar undulators and the quadrupole focusing lattice used in the experiment are replaced by a helical undulator in the second configuration. In this case, the electrons are matched into the weak focusing helical undulator.
Here we compare the new formulation with Guass-Hermite modes and the dynamic evolution of the polarization with the old formulation that uses a prescribed circular polatization with Gauss-Laguerre modes. In the third configuration, we replace the planar undulators used in the experiment with APPLE-II undulators which allows us to vary the polarization of the undulator line. Again, we compare the new formulation with the old formulation which uses a fixed elliptical polarization with Gauss-Laguerre modes.

The SPARC experiment
In this section, simulations using the old and new formulations are compared using the parameters of the SPARC experiment for which there are experimental measurements to anchor the two formulations. The experimental parameters are summarized in table 1. It was observed that the bunch charge was insufficient to reach saturation. In the simulation, the electron bunch was modeled using a parabolic profile with a peak current of 53 A and a width of 12.67 ps at the base. It was assumed in simulation that each undulator segment had one period entrance and exit taper and that the electron beam was matched to the undulator line. The simulation time window consists of 1600 temporal slices. For the parameters of table 1, the resonant wavelength is at 491.5 nm.  As a first validation, we used the experimental data to compare with both formulations. For this case, both formulations are essentially identical except that old formulation has a prescribed optical polarization and the new formulation allows self-consistent evolution of the optical polarization when the electron bunch propagates through the undulator line. A comparison of the simulated pulse energy as found in the old (shown by black exes) and the new (blue) formulations and from the experiment (red markers) is shown in figure 2, where the simulation results are averages taken over a total of 20 simulation runs, each with a different noise seed. This yields convergence to better than 5%. The experimental data is courtesy of L. Giannessi. Agreement between the two formulations is very good, although not identical, and agreement between the simulations and the measured performance is excellent. Energy conservation in both the old and new formulations is maintained to within better than one part in 10 3 .
Measurements of the 3rd harmonic radiation at the exit from the undulator line were also reported [25]. The measured pulse energy (red triangles) is shown in figure 2 along with the average pulse energy (blue circle) found over 20 simulation runs with different noise seeds. The error bars denote the range of fluctuations from shot-toshot or run-to-run. Note that while the measurements and simulations correspond to the pulse energy after 15 meters of the undulator line, we have offset the data points by +/− 20 centimeters respectively to facilitate the comparison. The substantial overlap in the 3rd harmonic pulse energies found in the experiment and in simulation indicate substantial agreement with the experiment was obtained using the new formulation.
It should be remarked that the large range of fluctuations in the 3rd harmonic pulse energies arise because the undulator line was too short to achieve saturation in either the fundamental or 3rd harmonic. Since the 3rd harmonic is generated through the nonlinear harmonic generation [17] mechanism which depends upon high intensities at the fundamental and the growth rate is three times that of the fundamental, small changes in the fundamental intensity can result in relatively larger fluctuations in the 3rd harmonic pulse energy prior to saturation.
The variation of the relative linewidth with propagation distance as determined from simulations with the old formulation (black exes), the new formulation (blue) and by measurement (red squares, data courtesy of L. Giannessi) is shown in figure 3, where the two results of the two formulations are very close. We observe that the simulated relative bandwidths are in excellent agreement with each other and that the measured linewidth is in substantial agreement with the simulations (to within 35% at a propagation distance of 15 m).
The planar undulators used in the simulation were oriented in such a way that the wiggle-motion was aligned with the x-axis. As such, we expect that the optical field generated in simulation will be linearly polarized along the x-axis as well. For this case, the old formulation assumes the light to be polarized in the x-direction. In contrast, the polarization of the optical field will evolve self-consistently in the new formulation and the on-axis Stokes parameters at the end of the undulator line yield θ=0 with an eccentricity ε=1 as expected.

Helical undulator simulation
The second case we consider is a bifilar helical undulator. For this case, the old formulation uses Gauss-Laguerre modes with a prescribed circular optical polarization. In contrast, the new formulation uses Gauss-Hermite modes and an undefined initial polarization. Since both the Gauss-Laguerre and Gauss-Hermite modes constitute complete sets, this difference should not be important as long as a sufficient number of modes is included.
For the simulations, we employ the same electron beam parameters used for the simulation of the SPARC experiment but replace the undulator/quadrupole line with a long undulator with the same period as used in the SPARC experiment that produces a field with a fixed helical polarization [27]. In order to preserve the resonance at 491.5 nm, the peak on-axis field of the helical undulator is reduced to 5.717 kG. In addition, the Twiss parameters are chosen to match the beam into the weak-focusing helical undulator [18,19]. As in the case of the  SPARC experiment, we simulate SASE for this configuration but only use a single (and identical) noise seed for both formulations which means that the initial phase spaces in each simulation will be identical.
Since there is no experiment which we can use to validate the formulation, we compare the present formulation with the results of simulation with the preceding model used in MINERVA [18,19], and we remark that the earlier formulation used a circularly polarized optical field representation based upon a superposition of Gauss-Laguerre modes in conjunction with a helical undulator. In contrast, the present formulation allows for the polarization of the optical field to vary based upon the electron motion induced by the undulator and is based upon a superposition of Gauss-Hermite modes. Since both the Gauss-Laguerre and Gauss-Hermite modes constitute complete sets, this difference should not be important as long as a sufficient number of modes is included.
Simulations of this configuration have been conducted with the earlier (old) formulation using 15 Gauss-Laguerre modes and with the present formulation using 12 Gauss-Hermite modes with 1600 temporal slices in each case. The growth of the pulse energy with propagation distance is shown in figure 4, and we observe that excellent agreement is found between the two simulations. As in the simulations for the SPARC experiment, energy conservation is maintained to within better than one part in 10 3 .
The spectra produced by the two formulations is also in good agreement. This is evident in figure 5 where we plot the spectra at the undulator exit as found by the old (a) and new (b) formulations.
These results show that, as long as sufficient modes are included, either Gauss-Hermite or Gauss-Laguerre modes can be used to model the interaction of the electrons with the optical field. Furthermore, the polarization evolved self-consistently within the new formulation into a circular polarized state with the handiness corresponding to the helical undulator used in the simulation.

The APPLE-II undulator
The third and final case we consider here is the APPLE-II undulator [12] that can be described using an analytic undulator model [18,19] and which consists of the superposition of two flat-pole-faced planar undulators oriented perpendicularly to each other and can be displaced with respect to each other. This configuration allows for adjustable polarization of the magnetic field, which is characterized by an ellipticity parameter u e [18] ranging from linear ( ) . Note that the APPLE-II undulator model [27] in the helical configuration differs from the pure helical undulator model in its off-axis field. The APPLE-II undulator consists in the superposition of two planar undulators which are shifted in relative phase. Hence, the field approximates the purely (bifilar) helical undulator only near the axis of symmetry and diverges from the pure helical undulator due to the different dependence on transverse coordinates of the two field models. For most high energy electron beams generated using radio frequency linacs, the crosssectional area of the beams is relatively small and this is usually not an issue.
In the simulations we replaced the planar undulators in the SPARC experiment with APPLE-II undulators of the same lengths and periods, and consider five different choices for the ellipticity: u e =0.00, 0.25, 0.50, 0.75, and 1.00. In order to ensure that the resonance remains at 491.5 nm, the on-axis field amplitudes for these cases were found to be 7.8796 kG, 7.6600 kG, 7.0577 kG, 6.3150 kG, and 5.5817 kG respectively as shown in figure 6 and as determined by Henderson et al [28]. In order to study the interactions out to saturation, and since the bunch charge was insufficient to reach saturation after 15 meters of undulators, we extend the undulator line to reach 28 meters in total length. We also note that the quadrupole positions and field gradients were held fixed as well as the initial Twiss parameters of the electron beam. This means that optimal focusing of the electrons through the undulator line may not be achieved in all three cases; however, the transmission of the electron beam is sufficient to illustrate the variations in the polarizations of the optical field.
We consider SASE using a single, identical noise seed for each of these five choices for the ellipticity. This ensures that the injected electron phase space is also identical for each of these cases. Furthermore, the simulation used 1600 temporal slices and the optical field was described using 15 Gauss-Hermite modes in the  The growth of the ensemble average of the pulse energies with propagation distance as found using the new formulation is shown in figure 7 for three choices of the ellipticity: u e =0.0, 0.50, and 1.00. Observe that the saturation distance increases (i.e., the growth rate decreases) as the ellipticity decreases from unity (corresponding to a circularly polarized undulator) to zero (corresponding to a linearly polarized undulator). This follows the decrease in the coupling coefficient associated with the increase in the oscillation in the magnitude of the transverse velocity. This is described by the generalized JJ-factor [28]. As in the cases for the planar and helical undulators, energy conservation is maintained to within better than one part in 10 3 .
An important validation measure is a comparison between the exponentiation (gain) length (L G ) found in both the old and new formulations with a theoretical estimate. A theoretical expression for the exponentiation length for elliptic undulators was given in [28] as a generalization of the parameterization described by Ming Xie [29,30] using the generalized JJ-factor. Using this generalization, a comparison between the gain lengths found in both simulations and theory for the five ellipticities under consideration is shown in table 2 which shows good agreement.
Using this present diagnostic for the Stokes parameters, the variation in the optical eccentricity versus the ellipticity of the undulator field is shown in figure 8 in blue while the prescribed polarization used in the fixedphasor representation in the old formulation is indicated by the red line. Since the optical field is composed of multiple temporal slices, each of which may behave somewhat differently, the average ellipticity shown in the figure is a weighted average over the power in each slice. As is evident in the figure, the eccentricity is unity and θ corresponds to the direction of wiggle-motion when the Apple-II undulator is configured for linear polarization, as expected, and decreases approximately linearly with increasing ellipticity but more slowly than would be expected from the fixed-phasor representation. We do not find that the eccentricity is zero when the Apple-II undulator is configured for an ellipticity of unity which is close to that of a helical field. It is important to remark here that the Apple-II undulator model is not the same as the helical undulator field used in Sec. III.B and, in general, diverges from the helical field slowly away from the axis of symmetry. Thus, it may not be possible to achieve pure circular polarization with an Apple-II undulator. The seeded FEL at Fermi Light Source in Trieste uses APPLE-II undulators and, while this may be due in part to the optical beam transport line from the end of the undulator line, they have indeed observed that the optical polarization is not purely circular [13] at the location of measurement when the ellipticity of the undulators is set to unity.
As in the aforementioned cases, this undulator line presents a fixed polarization over the entire distance; hence, a fixed phasor representation of the optical field may be appropriate. In order to test this, we consider a single polarization corresponding to a helical undulator with an ellipticity of unity. In addition, we consider a single noise seed which generates identical initial phase spaces for both formulations.
A comparison plot showing the pulse energy over the course of the undulator line obtained from the old and new formulations is shown in figure 9. The rms discrepancy over the entire undulator line is only about 9%. In addition, the rms spot size found using the two formulations shows a discrepancy of only about 0.8%. The polarization used in the old formulation is purely circular and, while these discrepancies are not large, these comparisons reveal the limitations of a fixed phasor representation when simulating APPLE-II undulators which diverge from a purely helical undulator off the axis of symmetry.

Summary and conclusion
Our purpose in this paper is twofold. First to present a numerical formalism that allows, for the first time among FEL simulation codes based upon the slowly-varying envelope approximation, for the self-consistent evolution of the optical polarization when the electrons propagate through an arbitrary undulator line. This formalism not only allows for the simulation of arbitrary combinations/sequences of undulators with varying magnetic polarizations, but also for the simulation of how imperfections in the magnetic orientation affect the FEL interaction; particularly on the purity of the polarization state when appropriate models for such imperfections are implemented. This can also be extended to study the effects of magnet degradation that may occur over time for permanent magnet-based undulators. Second, to study the utility of the concept and the functionality of the formalism. To this end, we described validations of the formalism by comparison with the experimental results of the SPARC experiment which employed planar undulators and by comparison with an earlier formalism for a helical undulator configuration, and found excellent agreement for both of these cases. We found good to excellent agreement for various properties of the optical field, such as pulse energy, gain length, and mode size and spectrum, when compared with the old formulation with prescribed optical polarizations. We also showed that energy conservation was maintained to within better than one part in 10 3 . However, we found that, with the new formulation, the polarization state starts to deviate from the fixed-phasor expectation when the ellipticity of the APPLE-II undulator increases from zero to one, i.e., the light is not circularly polarized when the ellipticity is unity.
We conclude that this formulation is able to faithfully simulate the control of the polarization of the output optical field for a variety of undulator combinations/sequences.
We have not addressed harmonic generation in this paper; however, this may be an important feature for many future light sources. In this regard, the MINERVA formulation includes both linear and nonlinear harmonic generation self-consistently because the electron trajectories are integrated using the complete Lorentz force equations including the transverse optical field elements; hence, harmonic generation is treated implicitly. Analytical treatments of harmonic generation in elliptical undulators have appeared in the literature [31][32][33][34][35]. Harmonic simulations with MINERVA are underway, which are complementary to the analytical models, and will be reported in a future paper.