Accelerator beam dynamics at the European X-ray Free Electron Laser

We consider the collective beam dynamics in the accelerator of the European X-ray Free Electron Laser for electron bunches with charges of 250 and 500 pC. The dynamics has been studied numerically with different computer programs, and we struggle to include properly the whole beam physics with collective effects. Several scenarios of the longitudinal dynamics are considered for the peak currents of 5 and 10 kA. A procedure for the choice of the machine parameters is described, and the tolerances of radio frequency parameters of accelerating modules along the accelerator are calculated. The bunch properties at the end of the accelerator before the undulator line are analyzed. A comparison between measurements and simulations for the beam energy loss due to wakefields is presented.


I. INTRODUCTION
A free electron laser (FEL) in the x-ray wavelength requires well-conditioned electron bunches with a small emittance, small energy spread, and high peak current. In order to prepare such bunches, the modern FEL facilities include very complicated accelerators based on radio frequency (rf) technology with several stages of beam compression.
The choice of the magnetic lattice layout and the rf module specifications are done during the design time. The beam dynamics studies at the design period of the European X-ray Free Electron Laser (XFEL) [1] are published, for example, in Refs. [2][3][4]. Now the European XFEL is constructed, and the space of beam dynamics parameters is fixed by technical constraints of installed hardware elements.
In this paper, we suggest a procedure for choosing "working points" in the space of parameters of the longitudinal beam dynamics under existing technical constraints. We show how to map these parameters into rf voltages and phases of the accelerating modules. In our studies, we combine analytical estimations with complex and time-consuming numerical iterative solutions in order to take into account self-consistent beam dynamics with collective effects.
Our studies are carried out for the current state of the European XFEL facility. We show with beam dynamics simulations that there are working points which produce the electron bunches with the design parameters at the entrance of the undulator lines and which will provide a high level of self-amplified spontaneous emission (SASE) radiation at hard x-ray wavelengths.
In comparison with our earlier publications [2][3][4], we use now three-dimensional treating of space charge forces along the whole accelerator. Additionally, we have included all sources of the impedances available as Green functions in the European XFEL database [5]. The working points suggested in this paper provide 30%-50% better rf tolerances as compared to ones suggested earlier in Refs. [3,4].
We start in Sec. II from the layout of the European XFEL and the technical constraints. Then, in Sec. III, the simulations in a rf gun are presented and analyzed. Section IV describes a simple model of longitudinal beam dynamics which is used in Sec. V to calculate rf tolerances of the accelerator modules. A procedure for the choice of working points of the longitudinal beam dynamics is developed in Sec. V as well. This procedure is applied in order to choose the working points for electron bunches with charges of 250 and 500 pC. An iterative algorithm [6] to convert the working point parameters in the rf parameters of the accelerating modules is described in Sec. VI and Appendix A. The code Ocelot [7] used in our studies is described shortly in Appendix B. Section VII describes the results of numerical modeling with an analysis of bunch properties at the end of the accelerator part of the European XFEL. Finally, Sec. VIII compares the obtained results with other simulations and measurements.

II. LAYOUT AND TECHNICAL CONSTRAINTS OF THE MACHINE
The European XFEL produces hard x rays up to the angstrom wavelength. It consists of a 17.5 GeV linear accelerator (linac) and several undulator lines for SASE radiation: SASE1, SASE2, and SASE3. The layout of the facility with the SASE1 branch is shown in Fig. 1. In order to compress the beam to a high peak current, the European XFEL electron beam line incorporates three vertical bunch compressors of C type.
In this paper, we are discussing the choice of the machine parameters in the linac for 250 and 500 pC bunch charges. The parameters are chosen based on the stability of the compression and the properties of the electron bunch after the compression.
The maximal accelerating voltages at different sections of the accelerator are listed in Table I. These values are taken from the current state of rf stations. The ranges for the compression compaction factors are listed in Table II [8].
In our numerical studies, we use the magnetic lattice parameters of the European XFEL facility from September 2018. The linear optics is shown in Fig. 2. It starts after the rf gun at position z ¼ 3.2 m from the gun cathode. The position z ¼ 3.2 m corresponds to the beginning of the booster rf module with eight Tesla superconducting cavities (see Fig. 1). The optics has a special phase advance between the last two bunch compressors in order to reduce effects due to coherent synchrotron radiation (CSR) kicks at them. It can be seen that the lattice has several additional dispersive elements which have to be taken into account when looking for working points with a desired global compression. In the current design, the total momentum compaction factor r inj 56 of the laser heater and the injector dogleg is equal to −35 mm. The second-order optics gives the second-order momentum compaction factor t inj 566 of these two elements equal to 94 mm. We will use these values in order to correct the momentum compaction factors of the first compression stage in the simple analytical model in Sec. IV.

III. SIMULATION OF RADIO FREQUENCY GUN
The electron bunch is produced by a shaped laser pulse in the rf gun. In our previous studies [2][3][4] of beam dynamics, we had assumed a flattop longitudinal shape of the laser pulse. In this work, we are using the Gaussian shape which corresponds to the current setup of the laser in the machine. The parameters used in the gun simulations for two charges are listed in Table III. The simulations are done with code Krack3 [9] and cross-checked with code ASTRA [10]. Both codes use the same mathematical model. In our simulations, we use 2 × 10 6 macroparticles, and the results have been crosschecked with simulations which use 2 × 10 5 macroparticles (see Sec. VIII). The slice parameters and the phase space projections at the distance 3.2 m from the cathode are shown in Fig. 3 for a bunch charge of 500 pC and in Fig. 4 for a bunch charge of 250 pC. In order to calculate the slice parameters, we have used 5 × 10 3 particles per slice. The "emittance" means in this paper the rms normalized emittance. For example, in the horizontal plane the projected normalized rms emittance is calculated from the total particle distribution by the formula ϵ proj where p x is the particle momentum and hi defines the second central moment of the particle distribution. We will use these particle distributions in the particle tracking and analysis described in the next sections. Let us note here that the projected emittance at this distance from the cathode is relatively large. The emittance will reduce in the booster considerably according with the emittance compensation process [11].

IV. ANALYTICAL MODEL OF MULTISTAGE BUNCH COMPRESSION
In this section, we introduce a one-dimensional model of multistage bunch compression in order to formulate in Sec. V a simple problem of longitudinal beam dynamics which allows an analytical solution: the conversion of parameters of longitudinal beam dynamics into rf parameters of accelerating modules. This solution will be used later on to build a "preconditioner" [12] in iterative tracking with full three-dimensional beam dynamics.
In order to describe the longitudinal beam dynamics inside the bunch during the compression and acceleration, let us introduce several coordinate systems. As a starting point, we consider a bunch of electrons after the rf gun at position z ¼ 3.2 m. The relative coordinate along the bunch will be noted as s. It has an origin at the position of the reference particle and increases in the direction of the bunch motion: The head of the bunch has a positive value of s, the tail has a negative one, and the reference particle is somewhere in the middle of the bunch (a mean position of all electrons). After the rf gun, the reference particle has energy equal to the slice energy at this position: Here and in the following, the low index i presents a reference point along the multistage bunch compression system (after the rf gun i ¼ 0), and the upper index 0 shows that the quantity is taken at the position of the reference particle.
Let us introduce the relative energy deviation coordinate In the longitudinal phase space ðs; δ 0 Þ, the energy distribution from the gun can be approximated by a third-order polynomial: where we use a suggestion that the reference particle after the rf gun has a design energy of δ 0 ð0Þ ¼ 0. Here and in the following, the prime denotes a partial derivative with respect to the particle position s (inside the bunch after the rf gun). The values of the derivatives in Eq. (1) for two bunches in Sec. III are listed in Table IV. They are obtained from the numerical fit to the particle distributions at the booster entrance, z ¼ 3.2 m from the gun cathode. The corresponding curves are shown by dots in Figs. 3 and 4. The reference particle has coordinate "0" at these plots. It can be seen from Figs. 3 and 4 that the reference particle position is not only the mean position of all macroparticles, but it corresponds approximately to the position of the peak current as well. If these two positions differ considerably, a better choice for the reference particle position could be the exact position of the peak current.
Let us consider the transformation of the longitudinal phase space distribution in a multistage bunch compression and accelerating system shown in Fig. 1. The system has three bunch compressors fBC 1 ; BC 2 ; BC 3 g and several accelerating sections fL 1 ; L 2 ; L 3 g. The injector section L 1 includes the booster and the third harmonic module.
In order to describe the longitudinal beam dynamics, we introduce several additional reference points to the already described one (after the rf gun). The longitudinal coordinate after bunch compressor number i will be denoted as  s i , and the energy coordinate at the position immediately after the bunch compressor will be denoted as δ i . The reference particle is always at the origin of the coordinate system. Note that throughout the paper the coordinate s (position in the bunch after the rf gun) will be used as an independent coordinate. All other functions depend on it.
For example, the function s i ðsÞ means that the particle with the initial position s (in the bunch after the rf gun) has the position s i after bunch compressor BC i . In the following, we omit the dependence on coordinate s in the notation. A phase space dilution due to rf phase effects occurs, because different electrons experience different rf phases as they pass through the rf modules. Thus, they gain different amounts of energy. For relativistic electrons interacting with a sinusoidally time-varying field, the energy gain of the electron is proportional to the cosine of the phase angle between its position and the position of maximum energy gain. Hence, the energy changes in the accelerating modules can be approximated as where e is the electron charge, φ 11 and V 11 are a phase and an on-crest voltage of the booster, respectively, φ 13 and V 13 are a phase and an on-crest voltage of the third harmonic module, respectively, φ i and V i are a phase and an on-crest voltage of accelerating section L i , respectively, and k is a wave number. The relative energy deviations in the reference points after the bunch compressors read The transformation of the longitudinal coordinate in compressor number i can be approximated by the expression Equations (1)-(4) present a simple nonlinear model of a multistage bunch compression system. For the fixed values of rf parameters and momentum compaction factors, we define the compression functions in each bunch compressor: The global compression function C 3 ðsÞ presents the compression after compressor BC 3 , which is obtained for the particles in the neighborhood of position s (the position in the bunch after the rf gun). For example, if we would like to increase the peak current by a factor of 50 at the position of the reference particle, then C 3 ð0Þ ¼ 50. In other words, function C 3 ðsÞ describes the increase of the current in the slice with initial position s. Above, we have introduced complementary functions fZ i g which we will call the inverse compression functions and which will be used in the analysis of rf tolerances in the next section. Finally, we would like to stress that the simple onedimensional model of the longitudinal beam dynamics described in this section does not include any collective effects. Additionally, it neglects the velocity bunching in straight sections due to the energy chirp in the bunch. The peak current of the bunch affects only the choice of the bunch compressor ratios fC i g as described in the next section.

V. RF TOLERANCES, SPACE CHARGE DEFOCUSING, AND CHOICE OF WORKING POINTS
The final bunch length and the peak current are sensitive to the energy chirp and, thus, to the precise values of the rf parameters. Let us calculate a change of the compression due to a change of the rf parameters. Unlike other publications on rf tolerances, we do not consider the voltage and the phase tolerances separately. Such onedimensional treatment is replaced in Ref. [6] by twodimensional consideration, where the phase and the voltage can change simultaneously and we estimate the tolerance from the worst possible direction of the change.
In the following, we are going to consider rf tolerances only in the booster. This tolerance is the tightest one (see Table VI). To simplify the notation, we define new variables for rf parameters through the complex number notation: In order to obtain a stable compression, we require that the relative change of the global compression C 3 at s ¼ 0 is smaller than Θ: where v 11 ¼ ðX 11 ; Y 11 Þ T is a vector of rf parameters and 0 stays for the design values. Following Ref. [6], we define the rf tolerance of the booster as θ 11 : where the gradient of the global inverse compression function Þ T can be found through iterative relations: Here ∂ means a derivative in X 11 or Y 11 , and Z 0 ≡ 1. In the same manner, we can define voltage and phase tolerances by the relations Such introduced tolerances are related by the simple relation The European XFEL has three bunch compressors, and we have to define eight rf parameters ðX 11 ; Y 11 ; X 13 ; Y 13 ; X 2 ; Y 2 ; X 3 ; Y 3 Þ in the machine. In order to define them, we have to fix 15 independent parameters of the longitudinal beam dynamics: 3 -compression factors after each bunch compressor and the first and the second derivatives of the global compression.
Instead of deflection radii in bunch compressors, we will use in some cases the compaction factors r 56i , i ¼ 1, 2, 3. And instead of derivatives of the global compression function, we will use derivatives of the inverse global compression function: Z 0 3 ; Z 00 3 . The conversion of these 15 parameters in the eight rf parameters is described in Appendix A. In the following, we will discuss an optimal choice of the parameters.
The initial conditions come from the gun simulations and are listed in Table IV. The energies at the beam compressors are fixed by design studies and are listed in Table V. The global compression C 3 usually comes from requirements on the peak current from FEL simulations. At the design report of the European XFEL [1], the peak current is fixed to 5 kA.
The first and the second derivatives of the global compression are special parameters which allow one to tune the flatness and the symmetry of the current profile. Initially, we will put them to zero. It means that we would like to have a flat current profile at the vicinity of the reference position s ¼ 0. Hence, at the end we need to find optimal values only for five parameters: deflecting radii in the bunch compressors r 1 , r 2 , and r 3 and compression ratios in the bunch compressors C 1 and C 2 .
Let us consider the choice of these parameters for the case of a bunch with a charge of 500 pC compressed to the peak current of 5 kA. In order to choose these parameters, we will do a scan of the parameter space. Let us describe the procedure. In our considerations, we will first define the maximal phases in sections L 1 and L 2 : where we use a safety margin S < 1 and V max 2 and V max 3 are maximal accelerating voltages from Table I in sections L 2 and L 3 correspondingly. In the following, we set parameter S to 0.9 in order to be able to compensate the self-field effects with additional shifts in voltages and phases later on.
The scan will be done numerically in three dimensions ðr 561 ; C 1 ; C 2 Þ, where under r 561 we understand the total momentum compaction factor from the gun (z ¼ 3.2 m) to the end of the compressor BC 1 : r 561 ¼ r inj 56 þ r BC 1 56 . Here r BC 1 56 is the momentum compaction factor of compressor BC 1 alone. We fix arbitrary the momentum compaction factor r 0 561 . After this, we do a two-dimensional scan of parameters C 1 and C 2 . For each pair C 0 1 ; C 0 2 , we calculate the maximal energy chirp before compressor BC 2 : The momentum compaction factor in compressor BC 2 is chosen as where we took into account the technical constraint in Table II. The maximal energy chirp before compressor BC 3 and the momentum compaction factor in BC 3 are found aŝ For the set of parameters ðC 0 1 ; C 0 2 ; r 0 561 ; r 0 562 ; r 0 563 Þ, we calculate analytically (as described in Appendix A) eight rf parameters ðX 0 1;1 ; Y 0 1;1 ; X 0 1;3 ; Y 0 1;3 ; X 0 2 ; Y 0 2 ; X 0 3 ; Y 0 3 Þ and use these rf parameters to calculate the rf tolerance θ 1;1 in the booster. In the following, we always show the tolerances for Θ ¼ 0.1, which means that the total compression deviation is restricted by 10%.
The two-dimensional plots of the rf tolerance in the booster are shown in Fig. 5 for two different compaction factors of compressor BC 1 : r BC 1 56 ¼ −40 mm and r BC 1 56 ¼ −60 mm. The dotted lines on the plots present the limits of the available rf voltage (see Table I). It can be seen that the smaller absolute value of the compaction factor r 561 results in larger tolerances, but the area between the dotted lines due to the technical constraints on the rf voltage is reduced. We will set r 561 to −40 mm, as yet smaller absolute values will be difficult to optimize in an even smaller area between the dotted lines. The working point is set to C 1 ¼ 3 and C 2 ¼ 21 as shown by the red dot in the plots. This choice is based on two arguments. From one side, we would like to have a relatively high compressed beam already before compressor BC 3 in order to use a larger deflecting radius in BC 3 and to reduce in such a way the CSR effects in it. Simultaneously, we would like to have a low compression in compressor BC 1 to avoid the strong space charge forces at the relatively low beam energy between compressors BC 1 and BC 2 .
In order to quantify the transverse space charge forces, we introduce the normalized space charge parameters: where ϵ x and ϵ y are normalized beam emittances, I 0 is the peak current, I A is the Alfvén current, γ is the Lorentz factor of the reference particle, and β x and β y are the beta functions of the linear optics (without space charge) shown in Fig. 2. These space charge parameters relate the space charge forces to the averaged focusing in the magnetic lattice (see [13]). Figure 6 presents the space charge parameters along the first (low-energy) part of the accelerator. We see that, for the chosen working point ðC 1 ¼ 3; C 2 ¼ 21Þ, the defocusing due to the space charge forces is approximately equal to or less than the focusing of the quadrupoles. In order to compress the bunch with a charge of 500 pC to the higher peak current of 10 kA, we change only the deflecting radius in the last bunch compressor BC 3 . As expected, the rf tolerance in the booster will be 2 times smaller, and the defocusing due to the space charge after the last bunch compressor will be 2 times stronger. The rf tolerance for this case is shown in Fig. 7, left. The space parameters [Eq. (8)] remain near or below 1 along the whole accelerator.
The working point for a bunch with a charge of 250 pC compressed to 5 kA peak current is shown in Fig. 7, right. This booster tolerance plot is quite similar to one on the left. Indeed, it is expected, as the total compression in both cases is approximately the same. The space charge defocusing is similar to the previous case, too, as the ratio of the peak currents to the emittances is approximately the same in both cases.
We summarize the parameters of the longitudinal beam dynamics for the three cases in Table V. The last row contains the parameters of the second derivative of the inverse global compression function as used in the simulations in the following sections.
The tolerances of the rf parameters of different accelerator modules for the three working points considered above are presented in Table VI. In order to estimate these tolerances, we have used the analytical solution of the onedimensional problem (see Appendix A). In comparison with our earlier publications [2][3][4], the suggested in this paper working points provide 30%-50% better rf tolerances. It is shown in Ref. [14] that the stability of the rf  parameters at the European XFEL is better than those required in Table VI. Before proceeding to the simulations, let us summarize the methodology to choose the working point at the threestage bunch compression system of the European XFEL at the existing technical constraints. As a first step, the energies in the bunch compressors E 0 1 , E 0 2 , and E 0 3 are fixed to the values in Table V. This choice is due to the maximal accelerating voltages available in the facility (see Table I). The total compression factor C 3 is chosen from the requirement to reach the peak current of 5 or 10 kA after the final stage of compression in BC 3 . The first derivative of the inverse compression function Z 0 3 is set to zero in order to have a "symmetric" current profile near to the peak current. We set the phases in sections L2 and L3 near to their maximal values [see Eq. (5)]. It allows us to produce the maximal energy chirps before bunch compressors BC 2 and BC 3 and, therefore, to use the maximal possible deflecting radii in them to reduce the strength of the CSR fields generated by the bunches. The deflecting radius r 1 in the first bunch compressor BC 1 is fixed to the maximal value allowed by the available voltage in the booster and the high harmonic module. It is shown in Fig. 5 that this choice results in better rf tolerances in the booster. After this, we scan the rf tolerance in the booster for different compression ratios C 1 and C 2 in the bunch compressors BC 1 and BC 2 . In the scan, we simultaneously set the deflecting radii in BC 2 and BC 3 to the maximal possible values [see Eqs. (6) and (7)]. Again, it is done to reduce CSR fields in the bunch compressors. The final values of compression ratios C 1 and C 2 are taken near to the maximum rf tolerance of the booster (this tolerance is most critical as can be seen from Table VI). The shifts of the working points (from the area of the maximal rf tolerance) are done in order to avoid strong transverse space charge forces after BC 1 (see Fig. 6) and to have an already high compression ratio after BC 2 . Finally, we choose some minimal positive value of the third derivative of the inverse global compression function Z 00 3 in order to avoid a strong spike in the current at the bunch head.

VI. TRACKING WITH COLLECTIVE EFFECTS IN THE ACCELERATOR
For the fixed initial beam parameters E 0 0 , δ 0 0 , δ 00 0 , and δ 000 0 , the fixed magnetic lattice (including the deflecting radii r 1 , r 2 , and r 3 in the bunch compressors) and some choice of rf parameters, we can introduce a map of rf parameters into the remaining eight parameters of the longitudinal beam dynamics: This map exists in the real machine, but in the following we will consider only a case when this map is realized by some tracking program. In order to describe the real machine as closely as possible, the tracking algorithm is quite complicated and includes the self-fields self-consistently. For a special case of the simple model in Sec. IV, the operator A can be inverted analytically as described in Appendix A. Let us denote this operator as A 0 and write the solution formally: The rf parameters x 0 will not produce the required compression: Aðx 0 Þ ≠ f . In order to find the solution of Eq. (9), we suggested in Ref. [6] an iterative algorithm, which we write here as We introduce in Eq. (10) an additional parameter λ; 0 < λ < 1, as suggested in Ref. [15]. In our simulations, we use mainly λ ¼ 1, as it provides the fastest convergence. Smaller values of λ should be used if a divergence appears or a smoother convergence is desired. In order to implement the particle tracking Aðx n−1 Þ, we use the code Ocelot. Because of the usage of the analytical solution A −1 0 , the algorithm converges very fast. Usually, 5-10 iterations are sufficient in order to find the rf parameters x which produce the desired compression scenario f .
The physical models and the numerical algorithms of Ocelot are described shortly in Appendix B. The numerical modeling of the accelerator beam dynamics presented below includes all sources of the impedances available as Green functions in the form of Eq. (B1) in the European XFEL database [5].

VII. NUMERICAL MODELING OF COLLECTIVE BEAM DYNAMICS AT THE EUROPEAN XFEL ACCELERATOR
In this section, we present the results of the beam dynamics simulations of the working points introduced in Sec. V.
In the first scenario, we compress the bunch charge of 500 pC to the peak current of 5 kA. The parameters of the longitudinal beam dynamics are listed in Table V. As was discussed in Ref. [6], the positive values of the second derivative of the inverse compression function Z 00 3 allow us to reduce the compression in the head and the tail of the electron bunch in order to avoid strong current spikes. In order to avoid a microbunching (which we are not able to model properly with the small number of macroparticles used), the energy spread in the bunch was increased in the laser heater, which is installed in the machine before the injector dogleg (see Fig. 1). The power of the laser is chosen to produce the energy modulation amplitude on axis of 8 keV and to have after the compression the rms slice energy spread near to 1 MeV. The rf parameters obtained with the iterative algorithm in Sec. V are listed in Table VII. This table contains three rows. In the first row, the rf parameters obtained analytically are given. These analytical parameters do not produce the required compression even without self-fields, as the model in Sec. IV does not take into account the velocity bunching at low energies. Hence, initially we track the particle distributions in Ocelot without any self-fields and use 3-5 iterations to correct the rf parameters. The corrected rf parameters which take into account the velocity bunching are listed in the second row in Table VII. Finally, we include all self-fields models available in the code (space charge, CSR, and wakes) and use again 3-5 iterations to correct the rf parameters. The final set of rf parameters which should be used in the real machine is listed in the last row in the table.
The results of the simulations with the rf parameters found above are presented in Fig. 8. The bunch is accelerated to the energy of 17.5 GeV in the main linac. We show the bunch shape and the slice parameters at the entrance of the SASE1 undulator line. The bunch has the peak current of 5 kA with a low slice energy spread (below 1 MeV) and small slice emittances (0.6 μm in the bunch core). The projected transverse emittances are listed in Fig. 8 as well. The larger projected emittance in the y plane can be explained by the vertical orientation of the bunch compressors: The CSR fields change the energies and the trajectories of the particles in the deflection plane.
In the second scenario, we compress the same bunch to the peak current of 10 kA. The power of the laser is chosen to produce the energy modulation amplitude on axis of 4 keV. The longitudinal dynamics is very close to the previous case. Only in the last bunch compressor, we use a smaller deflecting radius r 3 and a higher compression rate C 3 . The rf parameters obtained with the iterative algorithm in Sec. V are listed in Table VIII. The results of the simulation for the bunch accelerated to the energy of 17.5 GeV in the main linac are shown in Fig. 9. The bunch has the peak current of 10 kA with approximately the TABLE VII. The rf parameters for a charge of 500 pC compressed to 5 kA peak current.  8. The properties of the bunch with a charge of 500 pC compressed to 5 kA peak current at the SASE1 line undulator entrance. same horizontal projected and slice emittances as in the previous case. However, coherent synchrotron radiation in the last bunch compressor spoils the vertical emittances considerably.
Finally, we consider compression of the bunch with a charge of 250 pC to the peak current of 5 kA. The power of the laser is chosen to produce the energy modulation amplitude on axis of 4 keV and to have after the compression the slice energy spread again near to 1 MeV. The longitudinal dynamics is very close to the case of a bunch with a charge of 500 pC compressed to 10 kA peak current. The rf parameters obtained with the iterative algorithm in Sec. V are listed in Table IX. The results of the simulation for the bunch accelerated to the energy of 17.5 GeV in the main linac are shown in Fig. 10. The bunch has the peak current of 5 kA with a better slice emittance (0.4 μm in the bunch core) compared to both cases with a charge of 500 pC. Figure 11 presents the snapshots of the longitudinal phase space along the accelerator. A considerable energy chirp is used in the last compressor BC 2 . The chirp is compensated partially by the longitudinal wakefields in the linac. At the same time, the space charge creates in the main linac a chirp of alternative sign at the position of the peak current [see the snapshot in Fig. 11(e)]. The wake in the collimator section is strong and of a resistive nature (the form of the wake is proportional to the bunch shape). It can be seen in the last snapshot that the profile of the current is imprinted in the longitudinal phase space after the collimator.

VIII. COMPARISON WITH OTHER SIMULATIONS AND MEASUREMENTS
The simulations presented above solve the complicated nonlinear problem of self-consistent dynamics of charged FIG. 9. The properties of the bunch with a charge of 500 pC compressed to 10 kA peak current at the SASE1 line undulator entrance. IX. The rf parameters for a charge of 250 pC compressed to 5 kA peak current.  VIII. The rf parameters for charge 500 pC compressed to 10 kA peak current. particles in self-induced and external electromagnetic fields. In our studies, we have used the code Ocelot [7]. The physical models of this code are described in Appendix B.
We have spent a considerable amount of time in order to cross-check Ocelot with other codes. As a first step, we have tested that the beam dynamics without collective effects is correct. The first-and second-order optical functions from Ocelot have been compared with those obtained by codes MAD [16] and Elegant [17]. Then we have used the direct particle tracking in the code CSRtrack [18] to check that the chromatic effects are modeled correctly in Ocelot. For example, we have reproduced the results for the European XFEL injector dogleg optics with sextupoles described in Ref. [19].
The space charge forces are modeled in Ocelot with the help of a three-dimensional Poisson equation in the bunch frame. In order to use the iterative algorithm in Sec. VI, we should start from a particle distribution at the low energy of 6.5 MeV after the rf gun at position z ¼ 3.2 m from the cathode. We have used particle-in-cell code ASTRA [10] to check that the space charge forces and the rf focusing [20] in the accelerating modules are correctly modeled in Ocelot. The code ASTRA uses the same space charge model (Poisson equation in the bunch frame), but it solves the equations of motion by the Runge-Kutta method. It uses a discrete representation of the rf field in a cavity. The code Ocelot uses transport matrices to track the particles, and the rf field is described analytically. In Fig. 12, we show the projected horizontal emittance and the horizontal β function β x ¼ hx 2 iγ=ϵ proj  Fig. 13. The minor differences between ASTRA and Ocelot results are mainly due to different models of the rf fields used in these codes: ASTRA uses a rf field obtained from a numerical simulation of a Tesla cavity by a special electromagnetic code; Ocelot uses an analytical description of the rf curvature. These results confirm that we can use Ocelot already at the low-energetic part of the injector.
The projected and the slice emittances for a charge of 500 pC have been measured in the European XFEL injector [21,22] and agree with those obtained in the simulations presented in Figs. 12 and 13.
The modeling of dispersive sections with strong coherent synchrotron radiation effects has been tested by a comparison with the code CSRtrack. We have obtained a perfect agreement with the one-dimensional projected model in both codes [7]. Additionally, we have implemented in Ocelot a possibility to model curved trajectories with deflections in both transverse planes simultaneously. We need this possibility to model the extraction arc to SASE2, for example (see Fig. 15 and Ref. [23]).
The results presented in this paper are checked by varying of sampling parameters and by changing the number of macroparticles. Figure 14  by the deviation of a small number of particles available in the overcompressed branch in the head of the bunch (see Fig. 8).
Recently, it was verified [24] that the wakefield model used in Ocelot agrees with the experimental data. The experiment was done at the European XFEL for a bunch with a charge of 250 pC accelerated to an energy of 14 GeV. The bunch length was adjusted with help of the phase and the voltage in section L2. It was measured with a transverse deflecting structure after bunch compressor BC 3 . The energy of the bunch was measured at two positions shown in Fig. 15: after the undulator line SASE1 (the blue arrow) and after the undulator line SASE3 (the red arrow). The right plot in Fig. 15 shows the energy loss of the beam ΔE for different bunch lengths σ z . The simulation results are shown by solid lines. The measured data are presented by points. The results after SASE1 are shown in blue, and the results after SASE3 are shown in red. We see reasonable agreement between the measured and the simulated data.
Finally, our estimations of the slice emittances agree well with the results published earlier in Ref. [4]. The projected emittances in the simulations of this paper are smaller than those obtained in Ref. [4]. It is due to the usage of the new optics with special phase advances between the bunch compressors to compensate partially the kicks from CSR fields.

IX. SUMMARY
We have discussed a procedure to choose a working point of the longitudinal beam dynamics at the European XFEL. Following this procedure, we have considered three working points for bunch charges of 250 and 500 pC compressed to several kiloamperes. Beam dynamics simulations which include collective effects are used to find the rf parameters and estimate the bunch parameters before SASE undulators. In comparison with our earlier publications [2][3][4], the working points suggested in this paper provide 30%-50% better rf tolerances. Additionally, for the first time, our simulations include three-dimensional space charge forces and wakefields along the whole accelerator.
In this paper, we have presented only a few measurements and comparisons. In our future work, we will try to carry out a more detailed comparison of the theoretical simulations with measurements in the real facility.

ACKNOWLEDGMENTS
The authors thank W. Decking and E. Schneidmiller for helpful discussions.

APPENDIX A: ANALYTICAL SOLUTION FOR THE THREE-STAGE BUNCH COMPRESSION SYSTEM OF THE EUROPEAN XFEL
The rf parameters in the accelerator of the European XFEL can be found analytically if we neglect the selfforces and the velocity bunching. The rf parameters in the accelerating sections L 2 and L 3 read The rf parameters in the booster and in the third harmonic module can be found as where α 2 ≡ δ 00 1 ð0Þ and α 3 ≡ δ 000 1 ð0Þ are unknown yet derivatives of the energy curve after the first bunch compressor BC 1 . For an arbitrary derivative of the global inverse compression function Z 0 3 , the parameter α 2 ≡ δ 00 1 ð0Þ can be found by the relations ;ȳ i ¼ȳ i − kY ixi−1 ; i¼ 2; 3; Finally, for an arbitrary second-order derivative of the global inverse compression function Z 00 3 , the last parameter α 3 ≡ δ 000 1 ð0Þ can be found as Z 0 2 ¼ Z 0 1 − r 562 δ 00 2 − 2t 562 δ 0 2 ; Z 0 1 ¼ −r 561 α 2 − 2t 561 α 1 ; x 1 ¼ −6u 561 ðα 1 Þ 3 − 6t 56i α 1 α 2 ;ŷ 1 ¼ 0:

APPENDIX B: PHYSICAL MODELS AND NUMERICAL METHODS USED IN OCELOT
The tracking of particles in Ocelot is done in the same way as, for example, in Elegant [17]. Quadrupoles, dipoles, sextupoles, rf cavities, and other lattice elements are modeled by linear and second-order maps. The focusing effect of rf cavities is taken into account according to the Rosenzweig-Serafini model [20].
The space charge forces are calculated by solving the Poisson equation in the bunch frame. Then the Lorentztransformed electromagnetic field is applied as a kick in the laboratory frame. For the solution of the Poisson equation, we use an integral representation of the electrostatic potential by convolution of the free-space Green's function with the charge distribution. The convolution equation is solved with the help of the fast Fourier transform. The same algorithm for the solution of the three-dimensional Poisson equation is used, for example, in ASTRA.
The CSR module uses a fast projected one-dimensional method and follows the approach presented in Refs. [18,25]. The particle tracking uses matrices up to the second order. The CSR wake is calculated continuously through the beam line of arbitrary flat geometry. The transverse self-forces are neglected completely. The method calculates the longitudinal self-field of a onedimensional beam that is obtained by a projection of the real three-dimensional beam onto a reference trajectory. A smooth one-dimensional charge density is calculated by binning and filtering, which is crucial for the stability and accuracy of the simulation, since the instability is sensitive to high-frequency components in the charge density.
The European XFEL contains hundreds of sources of the coupled impedances. In order to obtain the wake functions of different elements, we have used analytical and numerical methods. The wake functions of relativistic charge have usually singularities and can be described only in terms of distributions (generalized functions). An approach to tabulate such functions and use them later to obtain wake potentials for different bunch shapes was introduced in Refs. [5,26,27].
The longitudinal wake function near the reference trajectory r a is presented through the second-order Taylor expansion w z ðr; sÞ ¼ w z ðr a ; sÞ þ h∇w z ðr a ; sÞ; Δri þ 1 2 h∇ 2 w z ðr a ; sÞΔr; Δri þ OðΔr 3 Þ; where we have incorporated in one vector the transverse coordinates of the source and the witness particles, r ¼ ðx 0 ; y 0 ; x; yÞ T , Δr ¼ r − r a , and s is a distance between these particles. For arbitrary geometry without any symmetry, the Hessian matrix ∇ 2 w z ðr a ; sÞ contains eight different elements: where we have taken into account the harmonicity of the wake function in coordinates of the source and the witness particles [28]. Hence, in the general case we use 13 onedimensional functions to represent the longitudinal component of the wake function for arbitrary offsets of the source and the witness particles near to the reference axis. For each of these coefficients, we use the representation [5] hðsÞ ¼ w 0 ðsÞ þ 1 C þ RcδðsÞ þ c ∂ ∂s ½LcδðsÞ þ w 1 ðsÞ; ðB1Þ where w 0 , w 1 are nonsingular functions, which can be tabulated easily, and constants R, L, and C have the meaning of resistivity, inductance, and capacitance, correspondingly. The wake potential for arbitrary bunch shape λðsÞ can be found by the formula where c is the light velocity in a vacuum, λ 0 is a derivative of charge density λ, and the asterisk means the convolution.