Efficient model for low-energy transverse beam dynamics in a nine-cell 1.3 GHz cavity

FLASH and the European XFEL are SASE-FEL user facilities, at which superconducting TESLA cavities are operated in a pulsed mode to accelerate long bunch-trains. Several cavities are powered by one klystron. While the low-level rf system is able to stabilize the vector sum of the accelerating gradient of one rf station sufficiently, the rf parameters of individual cavities vary within the bunch-train. In correlation with misalignments, intrabunch-train trajectory variations are induced. An efficient model is developed to describe the effect at low beam energy, using numerically adjusted transfer matrices and discrete coupler kick coefficients, respectively. Comparison with start-to-end tracking and dedicated experiments at the FLASH injector will be shown. The short computation time of the derived model allows for comprehensive numerical studies on the impact of misalignments and variable rf parameters on the transverse intra-bunch-train beam stability at the injector module. Results from both, statistical multibunch performance studies and the deduction of misalignments from multibunch experiments are presented.


I. INTRODUCTION
At FLASH (Free-Electron Laser in Hamburg) [1,2] and European XFEL (European X-Ray Free-Electron Laser) [3,4], superconducting 9-cell TESLA (TeV-Energy Superconducting Linear Accelerator) cavities [5] accelerate the electron bunches in pulsed operation. Due to the high achievable duty cycle, thus long radio-frequency (rf) pulse structure, bunch-trains containing up to 800 and 27000 bunches can be provided at FLASH and European XFEL, respectively. Several cavities with individual operational limits [6] are supplied by one rf power source. Within the bunch-train, the low-level-rf system [7,8] is able to restrict the variation of the vector sum of the accelerating gradient of one rf station sufficiently [9]. However, individual cavities have an intrinsic variation of rf parameters within one bunch-train, caused by the effects of beam loading and Lorentz force detuning [10]. Misaligned cavities in combination with variable rf parameters induce intrabunch-train trajectory variations which decrease the multibunch FEL performance significantly [10].
Understanding the complexity of intrabunch-train trajectory variations is therefore vitally important for a successful multibunch FEL operation at FLASH and European XFEL. In this paper we focus on the transverse dynamics in the injector module. Considering the low beam energy, thus beam sensitivity to off-axis fields, the injector module is of great importance in limiting intrabunch-train trajectory variations. There are several ways for getting a proper description of the transverse beam dynamics in a rf accelerating structure, e.g., using tracking algorithms [11] or simplified analytic models [12,13]. Assuming knowledge of the electromagnetic field distribution and the initial conditions of the particles, tracking provides accurate solutions, even for very low particle energies. Since the track step has to be small compared to a cell length, many steps are required, which needs considerable computation time for simulations with high dimensional parameter scans. Furthermore most tracking codes are not optimized for dealing with different bunches with different rf parameters. Established simplified analytic models on the other hand may calculate the beam transport by few matrix multiplications. However, they are based on assumptions, most importantly ultrarelativistic beams, which do not apply at most particle injectors. Thus, a major challenge is to set up a model for low particle energies γ ¼ ½10…200, which is simple enough in order to calculate its output within milliseconds, yet able to reproduce key features of rf dynamics such as rf focussing and coupler kicks. Our approach uses a combination of numerically calculated axially symmetrical beam transport matrices and discretized coupler kicks, coefficients of which are derived via a Runge-Kutta tracking algorithm using a high precision 3D field map of the TESLA cavity. The final model uses the matrix formalism to calculate the beam transport through an accelerating module consisting out of eight cavities in the order of ms for 400 bunches. Its output will be benchmarked against start-to-end tracking and experimental data. The short computation time of the model function allows for comprehensive numerical studies on the impact of misaligned accelerating structures and rf parameters on the intrabunch-train transverse dynamics in the injector module. It will be shown that a distinct determination of the misalignments from fitting experimental data is prevented by experimental constraints. However, the model will be used in order to give statistically meaningful predictions about the multibunch performance of the injector module regarding differential limits of rf parameters and misalignments. Figure 1 shows a schematic drawing of the TESLA cavity. It is a 9-cell standing wave structure of about 1 m length whose lowest TM mode resonates at 1.3 GHz. Two field maps are available for describing the electromagnetic field configuration.

II. MODEL DEVELOPMENT
The cavity fundamental mode is nominally described by its axial electric field E z ðzÞ ¼ E z ðr ¼ 0; zÞ. The corresponding field map is obtained from a 2D ½r; z simulation of the TESLA cavity and can be found in Ref. [14]. The field configuration Eðr ≠ 0; zÞ follows from the Maxwell equations. The bold letter indicates a vector. This axially symmetrical (RZ) field map describes the accelerating mode without geometric disturbances.
The 3D field map [15] describes this mode including the fields induced by both higher order mode (HOM) and power coupler and can be found in Ref. [16]. It is given as a table of sine and cosine like amplitudes, E cos r ðrÞ and E sin r ðrÞ, respectively, for a purely reflected wave, thus a wave traveling from the cavity into the waveguide and r ¼ ½x; y; z. The electric field of the reflected wave, E r ðr; tÞ, has the following time-space-dependency E r ðr; tÞ ¼ E cos r ðrÞ cosðωtÞ þ E sin r ðrÞ sinðωtÞ; ð1Þ with ω being the angular frequency. Using the Maxwell equations, the electric filed for the forward wave E f ðr; tÞ can be calculated by reversing time and follows as E f ðr; tÞ ¼ E cos r ðrÞ cosðωtÞ − E sin r ðrÞ sinðωtÞ: ð2Þ Let A ½f=r and ϕ ½f=r being the amplitude and phase of the forward and reflected wave to/from the power coupler, respectively. The overall electric field component for the general case with given accelerating voltage V 0 and phase ϕ with respect to the beam can then be calculated with from the 3D field map for the pure decay mode, thus no incoming wave.V r normalizes the field to the Eigenmodesolution of the field map. The voltage standing wave ratio describes the ratio between the difference of the forwarded and reflected wave in respect to the overall accelerating field. The magnetic component behaves analogously, using similar symmetry properties of the field components. Please note that Γ is defined at planes in the waveguide at which the standing wave mode has maximum field amplitude, e.g., at the end of the coupler antenna. Note also that beam loading induces a reflected wave which is not included in Eq. (3).
A. Beam transport in axially symmetrical cavities Bunches are described as single particles, hence only referring to their centroid dynamics. Space charge effects and intra bunch wakefields are not considered. The change of transverse coordinates of a particle induced by an axially symmetrical cavity can be written in terms of a matrix formalism as with u and u 0 holding the particle transverse input and output coordinates u ¼ ½x; x 0 ; y; y 0 , respectively, and M RZ being the beam transport matrix of the cavity. The derivation of the analytic beam transport matrix in Ref. [12] assumes the ultrarelativistic limit. The typical initial beam energy in the injector module is about 5 MeV. At this energy the solution of the beam transport equation (5) using the analytical matrix shows significant disagreement with tracking, as can be seen on the left-hand side of Fig. 2. Plotted is the difference of horizontal offset Δx in respect to an ASTRA [11] tracking for different initial beam energies. The accelerating gradient is 20 MV m −1 For each energy 5000 randomly distributed beam initial trajectories in the range x 0 ¼ AE½4 mm; 4 mrad are evaluated. In the energy range of the injector module the disagreement is of up to almost 500 μm at the end of one cavity. Taking into account the beam size of about 1 mm this is not acceptable for modeling the data. In order to rely on the matrix formalism for describing the beam transport, numerical adjustments on the transfer matrix have to be made. Using the Maxwell equations, a quasi-3D field map can be calculated from Ref. [14]. A Runge-Kutta algorithm is used to solve the equation of motion for one cavity for an ensemble of initial particles, entering the cavity at different offsets and angles. The calculation of the beam transport matrix then becomes a linear regression problem of the form arg min where u 0=1;i are the transverse coordinates of the beam at the entrance and at the end of the cavity. However, explicit information of the dependencies of the matrix elements on the initial beam energy E 0 , the amplitude V 0 and phase ϕ of the accelerating field is lost in the numerical process. It is therefore necessary to calculate the numerical beam transport matrix at a sufficiently fine grid of model parameters.
A cubic spline interpolation of each matrix element is used to obtain the transfer matrix at each point in the parameter space ½E 0 × V 0 × ϕ. The variation of the matrix elements with respect to the model parameters is moderate and smooth. The interpolated solution converges reasonably fast with increased grid points. For this work 14×7×11 ¼ 1309 grid points within E 0 ¼ ½5…150 MeV, V 0 ¼ ½13…30 MV m −1 and ϕ ¼ AE30°are chosen. The righthand side of Fig. 2 shows the difference between the solution of the beam transport equation (5) using the numerical transfer matrix and the previously described ASTRA reference. The agreement in the energy range of the injector module is in the order of 0.1 μm. It is well below the required accuracy, considering the corresponding beam size of about 1 mm. Therefore, the transverse dynamics related to axially symmetrical rf fields can reasonably well be described by the numerically derived beam transport matrix. The energy gain ΔE of a particle in the TESLA cavity is determined by the accelerating mode and is to a very good approximation independent of the coupler fields. It's dependency on the beam energy E 0 , accelerating phase ϕ and gradient V 0 is found to be well described via for the previously mentioned parameter range with fitted coefficients a i .

B. Discrete coupler kicks
HOM and power coupler break the axial symmetry of the cavity and influence the transverse beam dynamics considerably, as illustrated in Fig. 3. Plotted are two particle trajectories through the injector module, where for both trajectories the equation of motion is solved independently for both, the RZ and the 3D field map. The difference between the trajectories obtained with different field maps sums up significantly, especially on the horizontal plane. An appropriate description of the transverse motion therefore has to incorporate the field disturbances caused by the couplers. The transverse rf kick is the total beam transverse momentum change along the trajectory normalized by the longitudinal momentum of the beam. The integrated transverse field strength experienced by an ultrarelativistic paraxial particle does not only depend on the absolute distance from axis, r, but from x and y independently if the rotational symmetry is broken by couplers. The integrated transverse field induced by the couplers can be separated from the axially symmetrical rf focussing part of the field, V RZ , using the 3D field map of the TESLA cavity and extracting the monopole part. The real part of V coupler corresponds to a net deflection of the bunch centroid, whereas the imaginary part represents a kick which depends on the phase between the oscillating field and each particle. This time-dependent kick induced by the couplers distorts the longitudinal slices of the beam by a different amount, which results in an increase of the projected emittance [17]. The normalized kick k ¼½x 0 ;y 0 on a bunch centroid induced by a coupler can be calculated as [18] kðx; yÞ with E 0 being the beam energy and x and y the beam's transverse position at the coupler position. V 0 is the amplitude of the accelerating field and e the electron charge.Ṽ ⊥ is the normalized complex kick factor, defined asṼ with V ∥ ¼ R dz e z Eð0; 0; zÞ e i ωz c .Ṽ ⊥ holds the information of the axially asymmetrical field disturbances induced by the couplers. Its real part is plotted in Fig. 4 for both, the upstream and the downstream coupler. Different Γ, thus modes of cavity operation are evaluated. The kick induced by the upstream HOM coupler does not depend on Γ. This is due to the fact that the electromagnetic field away from the power input coupler is to a very good approximation described by a standing wave and is not affected by the ratio of the forward and reflected traveling wave. The static part of the downstream kick in respect to different Γ relates to the downstream HOM coupler. The Γ-dependent part relates to the power coupler, which primarily acts horizontally. Coupler kicks can therefore respond independently from the resonating accelerating field to variations of the forward and reflected traveling wave.
In order to describe the transverse beam dynamics properly we use the axial-symmetric beam transport matrices and insert linearized discrete kicks at a certain location. The kick k on a bunch's centroid induced by a coupler can therefore be expressed as The normalized complex kick coefficients V ij describe the normalized transverse deflection induced by the coupler and have to be found numerically, since the paraxial assumption of Eq. (8) is not fulfilled at low beam energy. The full beam transport equation of one cavity becomes where M RZ i are the axially symmetrical beam transport matrices calculated according to Eq. (6) between the corresponding reference planes. M RZ up describes the beam transport between the entrance of the cavity and the first coupler, M RZ center between the upstream and downstream coupler and M RZ down between the downstream coupler and the exit of the cavity. k up ðuÞ and k down ðuÞ evaluate the normalized upstream and downstream coupler kick, respectively, at the transverse coordinate u ¼ ½x; x 0 ; y; y 0 , such that kðuÞ ¼ ½x; x 0 þ k x ðx; yÞ; y; y 0 þ k y ðx; yÞ.
For low beam energy the V ij depend not only on the mode of cavity operation, thus the real and imaginary part of the voltage standing wave ratio Γ, but also implicitly on the particle's initial energy E 0 and on both amplitude V 0 and phase ϕ of the accelerating gradient, spanning a parameter space E 0 × V 0 × ϕ × ℜΓ × ℑΓ. The reason for this is the beam trajectory dependence on these parameters, since the ultrarelativistic limit is not reached. The parameter fit is done as follows: At every point in this parameter space a particle's centroid distribution is created at the entrance of the cavity. The particle distribution at the exit of the cavity is obtained via tracking using the 3D-field map. In addition, the particle distribution in the center of the cavity is recorded. This gives three reference distributions, before and after each coupler region, for which the coupler induced field disturbances are taken into account. Next, the axially symmetrical beam transport matrices are determined: the tracking is redone with the same initial particle distribution and the same rf parameters using the RZ-field map. This time, the particle distribution is recorded additionally at the coupler positions. Between each of these 5 points the axially symmetrical beam transport matrices FIG. 4. Real part of the normalized complex kick factor for the upstream (left) and downstream couplers (right) as a function of the transverse coordinates x and y. All vectors are scaled by the same amount in order to assure a quantitative comparison. The three colors in the right correspond to the case of pure filling (blue) of the cavity, e.g., no reflected wave and Γ ¼ −1, standingwave operation (red, Γ ¼ 0) and pure decay mode (yellow, Γ ¼ 1), where there is no incoming wave. Note that the kick induced by the downstream HOM coupler does not change for different Γ. The net effect of the power coupler is primarily horizontal.
are calculated according to Eq. (6). For both couplers the reference distributions are compared with the output calculated with the linear beam transport using Eq. (13). A fitting routine was used to then find the V ij which best describe the coupler kick of Eq. (12) using the ultrarelativistic limit [18] as a starting point. The global variation of the V ij were found to be described via with a n ðE 0 ; ϕÞ ¼ A n ðE 0 Þ · cos ϕ þ B n ðE 0 Þ · sin ϕ þ C n ðE 0 Þ; The ½w; y A n ; z A n ; y B n ; z B n ; y C n ; z C n ij are 49 constants for each coefficient V ij and were found with a fitting routine.

III. MODEL VALIDATION
The developed model is compared to the results of a start-to-end tracking using ASTRA [11] and to experimentally derived data at FLASH. The evaluation limits are

A. Comparison with ASTRA
The rf and beam input parameters are randomly created within the limits. The rms difference of the transverse position u for one cavity as a function of beam input energy is shown in Fig. 5 using different models for calculating the beam transport. At energies above 100 MeV the ultrarelativistic limit is by a very good approximation reached and the beam transport can be calculated according to Ref. [12] including coupler kicks according to Ref. [18]. Especially in the first cavities, however, it is important to use the fitted solutions for both the transfer matrices and the coupler kick coefficients V ij . Finally the beam transport through the whole injector module including eight cavities is studied. The initial beam energy is set to E 0 ¼ 5.6 MeV. Rf-and the beam input parameter are randomly created within the range of modeling as defined above. The rms difference of the output of the model function using the previously mentioned high energy approach is hΔu rms i ¼ 3 mm, the mixed approach results in hΔu rms i ¼ 1.3 mm. The difference of the low energy approach to the ASTRA reference is hΔu rms i ¼ 56 μm. Taking the beam size of about 1 mm into account this is a satisfactory result.

B. Comparison with dedicated experiments
This section gives a comparison of the model function with dedicated experiments at the injector section at FLASH. The experimental setup is illustrated in Fig. 6. Two beam position monitors (BPM) are located in the drift space between the solenoid magnet of the rf gun and the injector module. Therefore, the beam input parameter u 0 ¼ ½x 0 ; x 0 0 ; y 0 ; y 0 0 at z 0 ¼ 1.31 m can be calculated. The measured position at the BPM downstream the last cavity gives the reference The model function is evaluated for each bunch independently. The rf signals of the forward and reflected traveling wave are measured inside the waveguides at the circulators. The rf parameters required for the model in order to calculate the beam transport matrices and coupler kick coefficients are ½V refl ; ϕ refl ; V forw ; ϕ forw , with V refl=forw being the amplitude and ϕ refl=forw the phase of the reflected and forward traveling wave, respectively. The rf data recorded in the data acquisition system [19] is manually recalibrated FIG. 5. Rms difference between the output of the model function and an ASTRA tracking hΔu rms i as a function of initial beam energy E 0 , evaluated for one cavity and averaged over both transverse planes. The beam transport matrices are calculated with different models. Plotted is the developed model (blue, Eq. (13)), numerical transfer matrices with ultra-relativistic V ij (red, cf. [18]) and the analytic model (purple, cf. [12]). according to Ref. [20] in order to remove cross-coupling effects between the forward and reflected signals. The energy at the entrance of the module is calculated with phase and amplitude of the gun rf, ½V gun ; ϕ gun , assuming the gun calibration to accelerate the beam to 5.6 MeV=c at V gun ¼ 54 MV m −1 and ϕ gun ¼ −1.85°. In order to verify the model function, an experimental setup is designed in which the impact of unknown parameters, e.g., the offsets of BPMs, is supposed to cancel out. This is done by comparing different sets of rf parameters and BPM readings, as will be described later. Each measurement of BPM-and rf data is an average over approximately 300 consecutive bunch-trains to deal with short-term jitter. The bunch spacing is 1 μs.
If the accelerating gradient is changed, the transverse offaxis fields change as well as the on-axis accelerating field of the cavity. If the variation of the klystron power is slow enough to ensure a steady-state condition, coupler kicks should vary only in strength, not in direction and in inferior order. A reference is measured. The klystron power is modulated with f mod ¼ 3 kHz, assuring the accelerating field to be in resonance. The difference of BPM readouts in respect to the reference measurement for the horizontal and vertical plane, Δx and Δy, respectively, is plotted in Fig. 7 in black. The output of the model function is calculated and subtracted correspondingly and plotted in blue and red in Fig. 7. This, so to speak, partial derivative of the offset on the vector sum of the accelerating field shows a reasonable agreement with the BPM difference reading, considering that misalignments are not taken into account in the calculation.
As a second step, the implementation of coupler kicks is focused on. Caused by the limited bandwidth of the cavity, an increase of the modulation frequency of the klystron power will lead to a smaller amplitude of the modulation of the accelerating field, thus a higher reflected power. The ratio of the forward and reflected wave should vary and therefore, as pointed out in Fig. 4, the mainly horizontal forces induced by the power input coupler. Thus the impact of the power coupler compared to the overall transverse dynamics should increase with higher modulation frequencies. Especially misalignments should, for the most part, cancel out. Figure 8 shows the results of four measurements. The modulation frequency of the klystron power is increased subsequently from 3 kHz to 5 kHz, 50 kHz, and 100 kHz. Shown are the differences of the horizontal BPM readout between the modulated setup and the reference setup and the difference of the corresponding output of the model function in black and blue, respectively. The left column is calculated without coupler kicks, the right column shows the calculation including coupler kicks. Comparison between the columns in the last two rows points out that the beam dynamics above a modulation frequency of several 10 kHz are dominated by coupler kicks. These transverse dynamics are well described by the developed model.
It can be concluded that the implemented model function is both qualitatively and quantitatively able to reproduce the experimentally generated transverse trajectory features at the injector module. bunch-trains can be performed efficiently. In this section the deduction of misalignments from experimentally derived multibunch rf and BPM data is considered. It will be shown that the numbers of degrees of freedom, the limited amount of available data and the considered errors of the BPMs prevent reasonable fitting of the cavity and module misalignments.
Misalignments of involved structures are modeled by coordinate system transformations. The nomenclature used in this work is as follows: the beam parameters are Δu for the trajectory offset and Δu 0 for trajectory tilt angle with respect to the design axis, where u stands for the transverse planes x and y. The offset of a structure, for example the cavity or the module, is defined at its longitudinal center. The tilt of a cavity is evaluated around its longitudinal center, whereas the tilt of the module is evaluated around its entrance. The reference axis of the cavities is the axis of the module, while its reference is the design axis of the gun section. The misalignments of cavities and modules are Δu cav=mod and Δu 0 cav=mod for the offset and tilt angle, respectively. An algorithm was implemented, that allows simultaneous fitting of these parameters to an arbitrary amount of experimental data sets. This allows to increase the amount of experimental data, including the possibility that a certain set of rf and beam input parameters is insensitive to a particular misalignment, see Ref. [21] for more details.
Dedicated experiments [22] at the injector module at FLASH with 400 bunches are made to increase the sensitivity of single data sets to misalignments. The detuning of individual cavities and their loaded quality factors are subsequently changed. This results in a distinct change of rf parameters. Additionally, the beam input trajectory is varied. A total amount of 28 data sets are used for fitting. Each data set is averaged over 300 consecutive bunch-trains to deal with short-term jitter.
Despite the variation of rf parameters it was not possible to change the sensitivity of the model function with respect to misalignments of downstream cavities significantly. Apart from this, the amount of fit parameters, thus degrees of freedom, is large compared to the observable parameters. The consequential under-determination and the limitation of the fitting algorithm is evaluated first. Rf and beam input parameters from the above described data sets are used as input of the model function. Random misalignments are created. The output of the model function at the end of the module reflects a pseudo BPM reading and is used as reference for the fitting routine. Thus, it can be tested if the fit algorithm is able to identify the correct values of the misalignments.
In a first simplified example only cavity misalignments are studied. BPMs are supposed to work accurately and the module is on axis. A total set of 10 4 random misalignments according to the specification limits of Δu cav ¼ 0.5 mm and Δu 0 cav ¼ 0.5 mrad are evaluated. The initial values of the fitting routine are created randomly within the upper and lower bounds Δu bound ¼ AE0.5 mm and Δu 0 bound ¼ AE0.5 mrad, respectively. The rms difference between the simulated misalignments and the fitted misalignments is shown in Fig. 9 for the horizontal plane. The difference is significant. Note that the rms value of a real number which is randomly distributed between AE0.5 is about 0.3. It can therefore be concluded that based on the experimental data available it is not possible to fit the misalignments of individual cavities to a reasonable accuracy.
As a second step the prediction accuracy of the fit algorithm regarding the module misalignment is estimated. A total set of 10 4 random misalignments according to the uncertainties of Δu mod ¼ 5 mm and Δu 0 mod ¼ 0.5 mrad are evaluated. Analogously to the above mentioned simulation-and fit-procedure, the rms difference between the fitted and the simulated misalignments is calculated. Cavity misalignments and BPM errors are not considered. The rms error is found to be Δu mod;rms ¼ 0.51 mm and Δu 0 mod;rms ¼ 0.29 mrad. However, errors of the calibration factor δu BPM and offset Δu BPM of the involved BPMs can not be neglected. The prediction accuracy of the fit algorithm is expected to depend significantly on the accuracy of the BPMs. Its particular value is estimated as follows. The misalignment of the module is fitted for different bounds of the BPM error factor Σ. It incorporates the error of the calibration factor and the offset of each BPM and will be defined such that Σ ¼ n corresponds to δu BPM ¼ n·30% and Δu BPM ¼ n·0.5 mm. For example, Σ ¼ 1 reflects the possibility that all involved BPMs have a maximum calibration error of AE30% and maximum offset error of AE0.5 mm. The ratio between these values is reasonable for FLASH [23]. At each step Σ i , the prediction accuracy of the fitted module misalignment is estimated analog to the above described method and used for classification of the upcoming fit. Finally the unmodified recorded data sets are used for fitting. Fit parameters are the misalignment of cavities, the misalignment of the module in respect to the gun section and the calibration factors and offset errors of the BPMs. The boundary values for the misalignments are set according to the previously mentioned limits. The boundary values for the BPMs are incrementally increased according to the previously described ratio Σ. At each step Σ i , 500 different sets of randomly created starting points for the fit algorithm are independently evaluated. Solutions with χ 2 k;i < 2 · χ 2 best;i are accounted, where χ 2 best;i describes the lowest found value for each Σ i .
The results are shown in Fig. 10. The left-hand side shows the fitted module offset Δu mod in respect to the gun section as a function of the BPM error factor Σ. The horizontal and vertical plane are plotted in blue and red, respectively. The right-hand side shows the fitted values for the module tilt angle Δu 0 mod . The crosses are fitted solutions, while the dashed lines indicate the estimated prediction accuracy of the fit algorithm based on the previously described simulations. In the lower row χ 2 k;i is plotted for each solution. The initial χ 2 0 is in the order of 10 5 .
For perfectly calibrated BPMs, thus Σ ¼ 0, the fitted offset spans a range of about 2 mm, whereas the range of the fitted tilt angle is 800 μrad. Note that the longitudinal size of the module is 12 m. A tilt angle of 800 μrad therefore corresponds to an offset of about 5 mm on both ends in respect to the design axis of the accelerator. The range of the fitted values is unreasonably large. Subsequent studies on the accuracy of the BPMs located in the injector section indicate that Σ > 3. Furthermore long term drifts, most likely caused by technical issues in the readout electronics were identified. The actual BPM calibration during the multibunch measurements can therefore not be determined beyond reasonable doubt. It can be concluded that the recorded data sets are not distinct and comprehensive enough to deal with the large amount of degrees of freedom and the parameter limits of the regarded system. If the module misalignment should be determined sufficiently by multibunch-based misalignment measurements, several tasks would have to be accomplished. At FLASH, a quadrupole-triplet is located between the first and the second BPM downstream from the injector module. The difference between the theoretical and actual transfer matrix of this magnet is known to be significant [24]. If this problem is solved, both downstream BPMs could be used for the analysis. This results in full phase space information at the reference plane. In addition, dedicated studies on the calibration of the involved BPMs are advised.

V. STATISTICAL INVESTIGATION ON INTRABUNCH-TRAIN TRAJECTORY VARIATIONS
Statistical studies are made in order to generally describe the influence of misalignments in correlation with variable rf parameters on the intrabunch-train transverse dynamics. As a figure of merit it is reasonable to introduce the multibunch emittance blow-up in order to describe the ratio of the multibunch emittance ϵ MB in respect to the emittance of a single bunch ϵ SB . The multibunch emittance is calculated according to Refs. [25,26], assuming that only centroid dynamics differ within one bunch-train. The transverse single bunch emittance is assumed to be 1 μm and the Courant-Snyder parameters at the reference plane are assumed as α x ¼−0.076, α y ¼−0.049, β x ¼12.15 m, and β y ¼12.18 m. As described in this paper, the model requires the voltage standing wave ratio Γ as input. However, its variability within one bunch-train is not a common specification. The following workflow was developed in order to generate artificial rf data sets, while obeying the limitations set by the actual low level rf setup. The amplitudeV 0 and phaseφ of the vector sum of the accelerating field have to be chosen as well as the maximum slope of the amplitude ΔV 0 and phase Δϕ of the accelerating field of the individual cavities and their maximum detuning Δf within one bunch-train. The mean accelerating fields of the individual cavities are derived according to the waveguide-setup at the injector module at FLASH. The difference of the individual mean amplitude to the one eighth of the vector sum is calculated. The individual cavity accelerating field slope is randomly created, while keeping the vector sum constant. For each cavity, the signals of the forward and reflected wave have to be calculated. A detailed discussion on rf cavities can be found in Ref. [7]. We assume the special case of a superconducting cavity operating close to the steady state condition, nearly on crest, with beam loading and a detuning small compared to the resonance frequency (see Sec. 3.3.2 in Ref. [7]). The cavity voltage V cav , the forward-and reflected wave V for and V ref and the beam current I b are then related according to where the bold letters indicate complex numbers, for The ratio of the shunt impedance R and the unloaded quality factor Q 0 depends only on the geometry of a cavity, where the loaded quality factor Q L in this special case depends only on the coupling between the cavity and the waveguide system. The detuning angle ψ is defined as with f 0 being the resonance frequency of the cavity and f the operating frequency, leading to the detuning Δf ¼ f 0 − f. Using Eqs. (18) and (17) and assuming a constant charge and repetition rate of the bunches, the forward and reflected wave can be expressed as a function of the phase difference ϕ cav between the bunch and the cavity voltage with an amplitude of V cav and the detuning Δf. The voltage standing wave ratio Γ follows. With these rf parameters the transfer matrices and coupler kick coefficients of the model function are calculated for each bunch and each cavity individually. In order to avoid the loss of generality by analyzing only a coincidental case of parameters, a Monte Carlo filtering based method, also referred to as regional sensitivity analysis [27], is being followed. The goal is to analyze a multidimensional stochastic output statistically by conditioning the input space. The fundamental idea should be outlined on the basis of a simple example first. Consider the value of interest τ is a function of several parameters x, y. Consider the projection of τ on the τ−x-plane to be distributed as shown in Fig. 11 in the upper row on the left. In the upper right its correspondent cumulative distribution function (CDF) of τðxÞ is drawn. The CDF of τðxÞ accumulates the probability that the related variable takes values less than or equal to the evaluation point of the CDF. The interpretation of CDFðτðxÞ ¼ 2Þ ¼ 0.5 is therefore, that half of the evaluations of τðxÞ have a value smaller than or equal to 2. For the upcoming analysis a critical value τ c will be defined, such that CDFðτ c Þ ¼ p c . A reasonable value is p c ¼ 0.9, hence 90% of the evaluations of τ have values smaller than or equal to τ c . Note that in the first example of Fig. 11 the particular value of x at which τðxÞ is evaluated has no impact on the distribution of τ. Thus τ c is not sensitive to the range of x. Now consider the projection of τ on the τ−y-plane to be distributed as shown in the mid row of Fig. 11. Obviously the distribution of τ depends on the range of y. In this example, however, the corresponding CDF is calculated unconditionally, meaning that the whole parameter range of y is considered. The resulting CDF is identical to the one shown in the upper example. In the third example in the lower row, conditions on the input space are applied. The range of y is binned into two subranges, y 1 and y 2 , as highlighted with different colors in the lower row. The conditional CDF of τðyÞ depends on the range of the y i . Therefore the critical values of the two conditional CDFs, τ c;1 and τ c;2 , differ from each other. The sensitivity of the model to its parameters can be quantified in the variance of the critical values varðτ c;i Þ. It indicates if a model parameter, in this example y, is influential in determining the distribution of the model output τ (cf., the upper and lower row of Fig. 11). It is worth noting that the particular distribution of τ c;i depends on the parameter limits, the size of the subsamples and on the definition of p c . The previous concept will be applied on the analysis of the intrabunch-train transverse dynamics in the injector module. In this case the different parameters are: the maximum values of the intrabunch-train variation of the amplitude ΔV 0 and phase Δϕ of the accelerating field, respectively, and the variation of the detuning Δf. In addition, cavity offsets Δu cav and -tilts Δu 0 cav in respect to the module design axis as well as the offset of the module Δu mod and its tilt Δu 0 mod in respect to the design axis of the gun section are studied. In the following, ΔV 0 and Δϕ will be referred to as amplitude slope and phase slope, respectively.
The total number of model evaluations is 10 6 . The parameters are varied randomly within the range shown in Table I in the Appendix, giving a seven-dimensional distribution of τ-values. When conditions are applied subsequently to each parameter, the variation of τ c provides information about the distribution of τ in respect to that parameter. Figure 12 points out this one dimensional relationship for p c ¼ 0.9. For each parameter the sample size is divided into 5 subsamples. For example, the conditional limits on the offset of the injector module varies between Δu mod;1 ¼ ð0.5 AE 0.5Þ mm, Δu mod;2 ¼ ð1.5 AE 0.5Þ mm, and so on. The determination of the multi-bunch emittance blowup clearly is dominated by the amplitude slope, the detuning (for the horizontal plane) and the offset of the injector module with respect to the gun section, whereas other parameters are not significant in determining the distribution of τ. Results in Fig. 12 have to be taken cautiously. The range of each parameter is binned into relatively large subranges and single-parameter analysis does not reveal correlations between individual model parameter. Note also that the parameter limit for the module offset is ten times larger than the limit of individual cavity offsets. However, since the individual parameter range was set according to reasonable tolerance limits, this approach clearly indicates the dissimilar influence of parameters and assists in finding key parameters for the upcoming analysis.
For a more precise prediction the size of the subsamples has to be decreased, while conditions can be applied to multiple parameters simultaneously. Thus, analysis of τ c can be used to point out n-dimensional correlations. At first a parameter prioritization follows. The phase slopes of the accelerating field of individual cavities Δϕ as well as their offsets Δu cav and tilts Δu 0 cav are secondary in limiting the multibunch emittance for the current parameter range. Higher order mode based cavity misalignment measurements indicate, that the intended misalignment tolerances for individual cavities have been met throughout the modules, thus providing fixed evaluation limits of Δu cav ¼ 0.5 mm and Δu 0 cav ¼ 0.5 mrad. The maximum module tilt angle will be set to Δu 0 mod ¼ 0.5 mrad. The maximum slope of the phase of the accelerating field will be set to Δϕ ¼ 1°. Experimental observations show that this value is hardly exceeded during a typical user run with long bunchtrains. Fixing the above mentioned values remains a threedimensional parametrization of τ c ðΔV 0 ; Δf; Δu mod Þ, representing most of the variance of the multibunch emittance for realistic scenarios. This representation allows to study the three-parameter-CDF on τ, of which a twodimensional projection is plotted in Fig. 13 as a contour plot of τ c ðΔV 0 ; ΔfÞ for different offset limits of the injector module. The number of subsamples for each parameter is ten, splitting the total amount of 10 6 evaluations for each module offset into 10 × 10 ¼ 100 subsamples, each containing about 10 4 data points.
The superior influence of the amplitude slope compared to the detuning for the vertical plane is obvious. However, a critical amount of multibunch emittance is only exceeded when both amplitude slope and the module offset reach large values simultaneously. It is reasonable to assume that the evaluation limits are chosen large enough to include all possible machine realizations. The results can therefore be used in order to quantify the achievable performance improvement when reducing one of the parameter limits. If the goal, for instance, is to assure by a likeliness of 90%, that the horizontal multibunch emittance blowup is below 30%, one has to limit either the amplitude slope of the accelerating field to 0.7 MV m −1 or fix the module offset in respect to the gun section to better than 2 mm. Limiting the detuning, on the other hand, for example by means of a piezo-tuner to 10 Hz [28], would decrease τ c by about 15%, depending on the remaining module offset and amplitude slope.
Analysis of experimental data at FLASH [21] shows intrabunch-train trajectory variations up to 200 μm being induced in the injector module in both planes. This reflects a multibunch emittance blowup of about 20%. The rmsvalue of the amplitude slope of the accelerating field is about 250 kV m −1 with a detuning of about 100 Hz. The presented simulations are therefore interpreted as an indicator for a significant offset of the injector module with respect to the gun section, especially in the vertical FIG. 12. Single-parameter variance of τ c , calculated for 5 subsamples each, for the horizontal (blue) and vertical plane (red). Amplitude slope of the accelerating field and the offset of the module clearly dominate the variation of τ c . The detuning only impacts the horizontal plane. Note that the plotted results strongly depend on the defined parameter limits.
plane. This result agrees with the cautious estimate from fitting as presented in Fig. 10.

VI. SUMMARY AND CONCLUSION
An efficient model for calculating the beam transport in a TESLA cavity at low beam energy was found, including an analytical expression for describing discrete coupler kicks. Cross-checking between tracking and experimental data was performed. It can be concluded that the presented model is both qualitatively and quantitatively able to reproduce the transverse beam dynamics at low beam energy reasonably well.
The deduction of misalignments of accelerating structures from multibunch rf-and BPM data was considered and a fit algorithm was implemented. Thus far experimental constraints prevent accurate fitting of the misalignments of cavities and the injector module. On that account the developed model was used to analyze the interactions of intrabunch-train rf variations and structure misalignments statistically. Results indicate that the observed intra-bunch-train trajectory variation at the injector section at FLASH is caused by significant misalignments of the injector module with respect to the gun section, which is supported by the preliminary estimates of the fitting results. Investigations on procedures to reduce the misalignment are advised. It was furthermore possible to give accurate bounds for a performance study of the injector module regarding intra-bunch-train trajectory variations. It can be inferred that limiting the amplitude slope of the accelerating field is crucial for an improvement of the multibunch performance. Additionally Lorentz force detuning compensation would decrease the horizontal multibunch emittance blowup by 15% solely by limiting the variation of coupler kicks.

ACKNOWLEDGMENTS
This research was supported by the German Academic Scholarship Foundation. We would like to thank Vladimir Balandin, Siegfried Schreiber, Nicoleta Baboi, Jörg Rossbach and Christian Schmidt for valuable comments and fruitful discussions. We wish to extend our particular thanks and appreciation to Sven Pfeiffer, who with great skill and patience supported the measurements.