Transverse phase space characterization in an accelerator test facility PHYSICAL REVIEW ACCELERATORS AND BEAMS

The transverse phase space of a beam in an accelerator can be characterized using well-established methods based on observation of changes in the beam profile between screens at different locations along the beamline, or on observation of changes on a single screen for different strengths of upstream quadrupoles. Studies on CLARA FE (the Compact Linear Accelerator for Research and Applications Front End, at Daresbury Laboratory, UK) show that where the beam has a complicated (nonelliptical) distribution in transverse phase space, conventional analysis techniques aimed at characterizing the beam in terms of the emittances and Courant – Snyder parameters fail to provide a good description of the beam behavior. Phase space tomography, however, allows the construction of a detailed representation of the phase space distribution that provides a better model for the beam. We compare the results from three measurement and analysis techniques applied on CLARA FE: emittance and optics measurements from observations on three screens; emittance and optics measurements from quadrupole scans on a single screen; and phase space tomography. The results show the advantages of phase space tomography in providing a detailed model of the beam distribution in phase space. We present the first experimental results from four-dimensional phase space tomography, which gives an insight into beam properties (such as the transverse coupling) that are of importance for optimizing machine performance. of charge density from four-dimensional phase normalized phase space of four-dimensional emittance obtained from a Gaussian fitted to the four-dimensional (normalized) phase space distribution. tilt charge distribution in that refer to different degrees of freedom. reconstructed from analysis applied to simulated data based on the measured phase space distribution, technique. differences between analysis results the experimental overall very agreement in the phase space distribution found in each


I. INTRODUCTION
Knowledge of transverse beam emittance and optical properties are essential for the commissioning and performance optimization of many accelerator facilities. There are well-established techniques for emittance and optics measurements, often based on observation of changes in beam size in response to changes in strength of focusing (quadrupole) magnets, or observation of the beam size at different locations along a beam line [1,2]. Beam phase space tomography is also an established method for providing detailed information about the phase space distribution [3][4][5][6][7][8][9][10]. In this paper, we report the results of studies on CLARA FE (Compact Linear Accelerator for Research and Applications, Front End) at Daresbury Laboratory [11,12], aimed at characterizing the transverse phase space of the electron beam. The results of three different measurement techniques are compared, namely: beam-size measurements at three different locations along the beamline ("three-screen analysis"); measurement in the change of the beam size in response to the change in quadrupole strengths ("quadrupole scan"); and beam phase space tomography. At the time that the studies were carried out, the beam in CLARA FE had significant detailed structure in the phase space distribution (i.e., the beam could not be described by a simple elliptical phase space distribution). We find that in these circumstances, phase space tomography provides important insights into the transverse beam properties that cannot be obtained from the other techniques. Quadrupole scans can provide some useful information, but the results from three-screen analysis can vary widely, depending on the precise measurement conditions. Our studies of phase space tomography include the first experimental demonstration of beam tomography in four-dimensional phase space [13,14]. We find that this technique can provide information on coupling in the beam, which can be of value for optimizing machine performance [15]. This paper is organized as follows. In Sec. II we briefly review the definitions that we use for the emittances and optics functions in coupled beams, and the methods that we use for calculating these quantities. In Sec. III we describe the measurement procedures in CLARA FE.
The three-screen analysis method is discussed in Sec. III A, where experimental and simulation results are presented. The results show some limitations of the technique, and we discuss in particular why it does not produce reliable results when the beam has a complicated structure in phase space (i.e., the beam cannot be described by a simple elliptical distribution). In Sec. III B we describe the quadrupole scan analysis method, including application to measurement of the full covariance matrix in two (transverse) degrees of freedom. The quadrupole scan technique has some advantages over the three-screen analysis, but neither method can determine the detailed structure of the beam distribution in phase space. Such information can be provided by the final analysis technique, phase space tomography, which is considered in Sec. III C. We describe the implementation of phase space tomography on CLARA FE, including the use of normalized phase space [16], and show how tomography can be applied to determine the beam distribution in four-dimensional phase space [14]. Simulations are used to validate the technique, and experimental results are again presented. In Sec. IV we show the application of phase space tomography to provide a detailed characterization of the beam in CLARA FE under a range of machine conditions, looking at the dependence of the phase space distribution (including the coupling characteristics) on parameters of the electron source, including the strengths of the focusing solenoid and bucking coil, and the bunch charge. Given the complicated structure generally present in the phase space distribution of the beam in CLARA FE, phase space tomography provides important insights into the beam properties and behavior that would not be obtained from the three-screen or quadrupole scan analysis techniques. Tomography in four-dimensional phase space provides, in particular, information on transverse coupling that is of value for optimizing machine performance. Finally, in Sec. V, we summarize the key results, discuss the main conclusions, and consider appropriate directions for further work.

II. NORMAL MODE EMITTANCES AND OPTICAL FUNCTIONS
Since various definitions of beam emittance are used in different contexts, we briefly review the definition we use for the studies presented here, considering in particular the case where there is coupling in the beam. For clarity, however, we begin with the case of a single degree of freedom. Considering, for example, the transverse horizontal direction, the covariance matrix at a specified point in a beamline can be written: where x represents the transverse horizontal coordinate of a single particle at the specified location, p x is the horizontal momentum P x (at the same location) divided by a chosen reference momentum P 0 , and the brackets hi indicate an average over all particles in the beam. Note that we assume there is no dispersion in the beamline, so that the trajectory of the beam (nominally passing through the center of each quadrupole) is independent of its energy: for the present studies in CLARA FE, since the layout from the electron source to the end of the section where the emittance measurements are performed is a straight line, this will be a good approximation. The horizontal (geometric) emittance ϵ x and optical functions (Courant-Snyder parameters β x and α x and γ x ) can be calculated from: These relations imply that: and: Equation (6) defines an "emittance ellipse" in phase space. It is straightforward to extend these results to the vertical direction, to find the vertical emittance and Courant-Snyder parameters. If the beam distribution in phase space has elliptical symmetry, then the Courant-Snyder parameters and emittance are sufficient to describe the overall size and shape of the distribution. In general, by an "elliptical" distribution we refer to one for which the phase space density ρ is a function only of the betatron action, defined (in one degree of freedom, e.g., the transverse horizontal) by: For example, a Gaussian elliptical distribution can be described by: where N 0 is the total number of particles in the beam and, in this case, using Eqs. (6) and (8), the emittance ϵ x is equal to the mean action, i.e., ϵ x ¼ hJ x i.
Note that even for nonelliptical distributions, Eqs. (2)-(5) can be used to calculate values for the emittance and Courant-Snyder parameters. However, whether these values are meaningful or useful depends on the extent to which the phase space distribution can be characterized purely in terms of its second-order moments: the emittance and Courant-Snyder parameters are essentially a convenient way to parametrize the covariance matrix.
In principle, for an elliptical phase space distribution in one degree of freedom, the Courant-Snyder parameters and emittance at a particular location can be determined from just three appropriate measurements of the beam properties, for example, the rms beam size at three different locations along the beamline (corresponding to three different phase angles). This is the principle behind the three-screen measurement technique, discussed in Sec. III A. However, if the beam distribution in phase space has a more complicated structure, then the distribution cannot be described by just three parameters, and a larger number of measurements will be needed to determine the distribution. The results presented in Sec. III show that this is the case in CLARA FE. In such situations, a technique such as phase space tomography, discussed in Sec. III C, is needed to provide a good understanding of the beam properties and behavior.
In considering only a single degree of freedom, we assume that there is no transverse coupling in the beam or in the beamline, so that the transverse horizontal and vertical motions may be treated independently. The 4 × 4 covariance matrix for the transverse motion in this case takes block-diagonal form: As a beam passes along a beamline, coupling in the beam can be introduced and modified by skew components in the quadrupoles (for example, from some alignment error in the form of a tilt of the magnet around the beam axis) or from a solenoid field either at the particle source or further down the beamline. Then, the general covariance matrix for the transverse phase space distribution takes the form: The coupling can be characterized by the values of the cross-plane elements (hxyi, hxp y i, hp x yi, and hp x p y i) in the 4 × 4 covariance matrix, or by a generalization of the Courant-Snyder parameters (for example, using the methods in [17][18][19]). In the absence of coupling, the elements of the covariance matrix (10) can be determined from six measurements, such as the horizontal and vertical beam sizes at three different locations along the beamline. However, when coupling is present, the covariance matrix (11) has ten independent, nonzero elements, and ten measurements are therefore needed to determine all the elements. Screens at three different locations would provide just nine measurements (the beam sizes hx 2 i, hy 2 i and tilt hxyi at each screen); measurement of the full 4 × 4 covariance matrix can be accomplished using screens at additional locations (for example, [10,20]). Alternatively, the beam sizes and tilt on a single screen can be measured for a number of different strengths of the upstream quadrupoles (see, for example, [21,22]). This is the technique used in quadrupole scan and tomography measurements, discussed in Secs. III B and III C. With coupling present, the emittance calculated using Eq. (2) will not be the most useful quantity, since it will not be constant as the beam travels along the beamline. The conserved quantities where coupling is present are the normal mode emittances (or eigenemittances) ϵ I and ϵ II , where AEiϵ I;II are the eigenvalues of ΣS, with Σ the 4 × 4 covariance matrix (11), and S the antisymmetric matrix: As already mentioned, various formalisms have been developed for generalizing the emittance and Courant-Snyder parameters from one to two (or more) coupled degrees of freedom. Here, we use the method presented in [19], in which the ði; jÞ element of a covariance matrix Σ is related to the normal mode emittances ϵ I , ϵ II and corresponding optical functions β I ij , β II ij by: The optical functions can be obtained from the eigenvectors of ΣS. If U is a matrix constructed from the eigenvectors (arranged in columns) of ΣS, then: where Λ is a diagonal matrix with diagonal elements corresponding to the eigenvalues of ΣS. If the eigenvectors and eigenvalues are arranged so that, in two degrees of freedom: then the optical functions are given by: where k ¼ I; II, and: The covariance matrix Σ can then be expressed in terms of the normal mode emittances and optical functions using (13). When there is no coupling, the normal mode emittances and optical functions correspond to the usual quantities defined for independent degrees of freedom. For example, where the transverse horizontal motion is independent of the vertical and longitudinal motion, then: and: Since the normal mode emittances and optical functions depend on the eigenvalues and eigenvectors of ΣS, whereas the "uncoupled" emittances and Courant-Snyder parameters are calculated from the 2 × 2 submatrices along the diagonal of the covariance matrix Σ, there is no simple relationship between the coupled and uncoupled quantities in the general case. Finally, we note that if the optical functions β k ij at a given point s 1 in a beamline are known, the optical functions at any other point s 2 are readily computed using: where M 21 is the transfer matrix from s 1 to s 2 (calculated, for example, from a computational model of the beamline). The normal mode emittances and optical functions defined as described here, therefore provide convenient quantities for describing the variation of the beam sizes hx 2 i, hy 2 i (and other elements of the covariance matrix) along a given beamline. However, as discussed earlier in this section, the elements of the covariance matrix only provide a useful description of the phase space distribution in simple cases (e.g., for elliptical distributions). For complicated distributions with significant structure, alternative ways of describing the distribution may be needed.

III. MEASUREMENTS IN CLARA FE
Ultimately, CLARA is planned as a facility that will provide a high-quality electron beam with energy up to 250 MeV for scientific and medical research, and for the development of new accelerator technologies including (with the addition of an undulator section) the testing of advanced techniques and novel modes of FEL operation. So far, only the front end (CLARA FE) has been constructed: this consists of a low-emittance rf photocathode electron source and a linac reaching 35.5 MeV=c beam momentum. The layout of CLARA FE is shown in Fig. 1. The electron source [23] consists of a 2.5 cell S-band rf cavity with copper photocathode, and can deliver short (of order a few ps) bunches at 10 Hz repetition rate with charge in excess of 250 pC and with beam momentum up to 5.0 MeV=c. The source is driven with the third harmonic of a short (2 ps full-width at half maximum) pulsed Ti: Sapphire laser with a pulse energy of up to 100 μJ. The typical size of the laser spot on the photocathode is of order 600 μm. The source is immersed in the field of a solenoid magnet which provides emittance compensation and focusing of the beam in the initial section of the beamline. A bucking coil located beside the source cancels the field from the solenoid on the photocathode in the region of the laser spot.
The studies reported in this paper are based on measurements made in the section of CLARA FE following the linac, at a nominal beam momentum of 30 MeV=c. Measurements were made under a range of conditions including various bunch charges, and different strengths of the solenoid and bucking coil at the electron source. Three techniques were used, to allow a comparison of the results and evaluation of the benefits and limitations of the different methods. The first technique, the three-screen measurement and analysis method (described in more detail in Sec. III A, below), is based on observations of the transverse beam profile at three scintillating (YAG) screens, shown as SCR-01, SCR-02 and SCR-03 in Fig. 1. The quadrupole scan (Sec. III B) and tomography (Sec. III C) methods use only observations of the beam on SCR-03, though observations on SCR-02 were also made, in order to validate the results. For each of the three methods, two quadrupoles (QUAD-01 and QUAD-02) between the end of the linac and SCR-01 were used for setting the optical functions of the beam on SCR-01, and were kept at fixed strengths during data collection. As discussed later in this section, the strengths of these quadrupoles were chosen to produce a beam waist at SCR-02, based on a model using the design parameters for the upstream components (including the electron source and linac). A collimator is located between SCR-01 and SCR-02, but this was not used during the measurements. For all three measurement techniques, the strengths of three quadrupoles (QUAD-03, QUAD-04, and QUAD-05 in Fig. 1) located between SCR-02 and SCR-03 were varied. For the three-screen analysis, only the beam sizes for one set of magnet strengths are strictly needed to calculate the emittances and optical functions; however, as described in Sec. III A, measurements with different sets of quadrupole strengths can be used to validate the results by showing the consistency of emittance and optics values obtained for different strengths.
In the case of the quadrupole scan and tomography methods, SCR-03 provides the necessary data, and is referred to as the "observation point." For ease of comparison, for all three techniques we construct the covariance matrix at SCR-02, which is referred to as the "reconstruction point." Initial values for the quadrupole strengths to be used during a scan were chosen to provide a wide coverage (from 0 to 2π) of the horizontal and vertical phase advances from SCR-02 to SCR-03: in principle, this provides good conditions for reconstruction of the phase space using the quadrupole scan and tomography analysis techniques, since the projection of the phase space onto the coordinate axes is then observed over a wide range of phase space rotation angles. The phase advances are calculated for given quadrupole strengths using the nominal Courant-Snyder functions at SCR-02, based on the design parameters of the machine. Simulations were then performed to optimize the number of quadrupole strengths used in a scan, and the quadrupole strengths themselves, to allow accurate reconstruction of the phase space using the minimum number of different settings for the quadrupoles. Reducing the number of points in a scan reduces the time needed to collect data, but has an adverse impact on the accuracy of phase space measurements.
For the experiments reported here, each quadrupole scan consisted of 38 sets of quadrupole strengths. The quadrupole strengths and corresponding phase advances (from SCR-02 to SCR-03) are shown in Fig. 2.  At each point in a quadrupole scan on CLARA FE, ten screen images were recorded on successive machine pulses (at a rate of 10 Hz, with a single bunch per pulse): this allows an estimate to be made of random errors arising from pulse-to-pulse variations in beam properties. A background image was recorded without beam (i.e., with the photocathode laser blocked), so that any constant artefacts in the beam images, for example from dark current, could be subtracted. The rms beam sizes were calculated by projecting the image onto either the x or y axis, with coordinates measured with respect to a centroid such that: Average quantities are calculated from a beam image by integration of the image intensity with an appropriate weighting, for example: where Iðx; yÞ is the image intensity at a given point on the screen. The screens and cameras are specified to provide images with a linear relation between beam intensity and image intensity (i.e., without saturation of the images) up to bunch charges of 250 pC, with beam sizes in the range achievable on CLARA FE. For the measurements reported here, where bunch charges of up to 50 pC were used (with most measurements at 10 pC) the beam intensity was significantly below the point at which saturation occurs. Between each quadrupole scan, the quadrupole magnets were cycled over a set range of strengths to minimize systematic errors from hysteresis. Remaining sources of systematic errors include calibration factors for the magnets (when converting from coil currents to field gradients), magnet fringe fields, calibration factors for the screens, and accelerating gradient in the linac. It was found that better agreement between the analysis results and direct observations (used to validate the results) could be obtained if the beam momentum in the model used in the analysis was reduced slightly from the nominal 30 MeV=c. In the results presented here, a momentum of 29.5 MeV=c is used. It should also be noted that some variation in machine parameters (including rf phase and amplitude in the electron source and the linac) is likely to have occurred during data collection, and because of the time required to re-tune the machine it was not always possible to confirm all the parameter values between quadrupole scans.
In principle, for each of the three measurement and analysis methods, the strengths of the quadrupoles between SCR-02 and SCR-03 can be chosen randomly. However, if the profile of the beam on any of the screens becomes too large, too small, or very asymmetric (with large aspect ratio) then there can be a large error in the calculation of the rms beam size. Before collecting data, therefore, simulations were performed to find sets of magnet strengths, with fixed QUAD-01 and QUAD-02 and variable QUAD-03, QUAD-04 and QUAD-05, for which the transverse beam profiles on each of the three screens would remain approximately circular, and with a convenient size. It is also worth noting that, from (2), a large value of α x at a given location can indicate a large value for hxp x i at that location: calculation of the emittance then involves taking the difference between quantities that may be of similar magnitude, leading to a large uncertainty in the result. A further constraint, therefore, was to find strengths for QUAD-01 and QUAD-02 that would provide a beam waist in x and y (i.e., with α x and α y close to zero) at SCR-02 (the Reconstruction Point). Finally, quadrupole strengths were chosen to provide a wide range of horizontal and vertical phase advances from SCR-02 to SCR-03: this is a consideration for the tomographic analysis, and is discussed further in Sec. III C. Simulations to find sets of suitable strengths for all five quadrupoles were carried out in GPT [24], tracking particles from the photocathode (with nominal laser spot size and pulse length) to SCR-03, using machine conditions matching those planned for the experiments.
Space charge effects were included in the tracking simulations in GPT [25,26], though these effects are only really significant at low momentum, upstream of the linac. Space charge was not included in the analysis (using any of the three methods discussed in this section). To justify the assumption that space charge can be neglected in the analysis, we can compare the perveance with the beam emittance. In the case of a beam with an elliptical distribution in phase space, the evolution of the rms beam size σ x with distance s along a beamline is described by the envelope equation: where k 1 is the local quadrupole focusing, and the perveance K is: Here, β 0 is the particle velocity (scaled by the speed of light), is the relativistic gamma factor, and I A ≈ 17.045 kA is the Alfvèn current. For CLARA FE, we assume representative values of 30 MeV for the beam energy, 2 A peak current (corresponding, for example, to 10 pC bunch charge in a Gaussian longitudinal distribution with standard deviation 2 ps), beam sizes σ x ≈ σ y ≈ 0.3 mm, and normalized emittance 1 μm. With these values we find K ≈ 1.2 × 10 −9 , and 4ϵ 2 x =σ 2 x ≈ 1.2 × 10 −8 . We therefore expect to see space charge effects become significant only at bunch charges larger by some factor above the nominal 10 pC at which most of the measurements reported here were made.
Additional nonlinear and intensity-dependent effects may come from wakefields. Although wakefields in the linac may be significant, in the section of CLARA FE used for the phase space measurements, the vacuum chamber is essentially smooth and wakefields are small. With the bunch charges used in the experiments reported here (up to 50 pC) wakefield effects in the beamline following the linac are not expected to have any observable impact on the beam.

A. Three-screen method
The three-screen analysis method is based on the principle that, given the rms beam size (in either the transverse horizontal or vertical direction) at three separate locations, and knowing the transfer matrices between those locations, it is possible to calculate the covariance matrix characterizing the phase space beam distribution. As discussed in Sec. II above, this technique is appropriate when the beam has an elliptical distribution (or a distribution that is close to elliptical) in phase space, and there is no coupling. In such cases, the 2 × 2 covariance matrices of the second-order moments of the beam distribution in the transverse horizontal and vertical phase spaces can be treated independently, and provide a good representation of the beam properties.
The three-screen measurement and analysis technique has the advantage over other methods (such as tomography, discussed in Sec. III C) that data collection and analysis can be performed rapidly. CLARA FE includes a number of screens in appropriate locations following the linac, allowing this technique to be readily applied. In principle, the beam in CLARA FE should have a simple (close to elliptical) distribution in transverse phase space, with little or no coupling. The distribution should be well characterized by the 2 × 2 transverse horizontal and vertical covariance matrices, and the three-screen analysis technique was therefore used during machine commissioning. However, it was found that the results did not provide a good description of the beam behavior, and that values found for the emittances and optical functions could vary significantly, depending on the strengths of the quadrupoles between the screens. Results from phase space tomography show that at the time the measurements were made, the beam distribution had a complicated structure in phase space. This is also evident from the images observed on SCR-02, which consistently show a lack of elliptical symmetry in the coordinate space projection at this screen for any of the machine configurations that we studied (some examples are given, and compared with predictions from the tomography analysis to be described later, in Fig. 9). Here, we report the results of the three-screen measurements mainly to illustrate the limits of the technique in this case.
Let us consider just the transverse horizontal motion: the application to the vertical motion is straightforward. If the transfer matrix from one location in the beamline s 1 to another location s 2 is M 21 (with transpose M T 21 ), then: where Σ 1 is the covariance matrix at s 1 and Σ 2 is the covariance matrix at s 2 . For the present case where we consider motion in just one degree of freedom, the covariance matrices and transfer matrices are all 2 × 2 matrices. At a third location s 3 : where Σ 3 is the covariance matrix at s 3 , and M 32 is the transfer matrix from s 2 to s 3 . The elements of the transfer matrices can be calculated using a linear model of the beamline, with known quadrupole strengths.
Given measured values of hx 2 1 i, hx 2 2 i and hx 2 3 i (from observation of the beam images on the three screens), using (25) and (26) we can find hx 2 p x2 i and hp 2 x2 i from: where a 21 and b 21 are (respectively) the (1, 1) and (1, 2) elements of M −1 21 , a 32 and b 32 are (respectively) the (1, 1) and (1, 2) elements of M 32 , and: The same formulas can be applied to calculate the elements of the covariance matrix for the vertical motion, and (as discussed in Sec. II) the technique can be extended using additional screens or by changing the strengths of selected quadrupoles, to find the elements of the 4 × 4 covariance matrix for motion in two degrees of freedom. In applying the formulas (25)- (28) to CLARA FE, s 1 , s 2 , and s 3 correspond to the locations of screens SCR-01, SCR-02, and SCR-03, respectively (see Fig. 1). If there is no coupling in the beam, and if the distribution in phase space has elliptical symmetry, then the results may be validated by repeating the measurements for different strengths of the three quadrupoles between screens SCR-02 and SCR-03: since the magnets upstream of SCR-02 remain at constant strength, measurements made with different strengths of downstream magnets should all yield the same values for the emittances and Courant-Snyder parameters at this screen. Figure 3 shows a typical set of results from the threescreen analysis applied to CLARA FE, for the transverse horizontal and vertical directions, respectively. Elements of the covariance matrix [each scaled by an appropriate Courant-Snyder parameter, calculated from the covariance matrix elements using (3)- (5)] are plotted as functions of the phase advance from SCR-02 to SCR-03 in the respective plane (corresponding to different strengths of the quadrupoles QUAD-03, QUAD-04, and QUAD-05). In the case of a simple elliptical distribution with no coupling, from Eqs. (3)-(5) we see that scaling the covariance matrix elements by the Courant-Snyder parameters should give values that are independent of the phase advance, and equal to the geometric emittance. However, the results in Fig. 3 show significant variation in each of the scaled elements of the covariance matrix over the range of the quadrupole scan: this is particularly evident in the horizontal direction, and is reflected in the values calculated for the emittance and optics functions. In both the horizontal and the vertical planes, a number of points in the quadrupole scan lead to imaginary values for the emittance (the covariance matrix has negative determinant), or nonphysical negative values for the covariance matrix element hp 2 x i. These points (21 points in the horizontal plane, and 7 in the vertical plane) are omitted from the plots in Fig. 3 and from the calculation of mean values of emittances and optics functions. The omitted points all have phase advance greater than 1.4 rad in the respective plane.
As discussed in Sec. II, if the beam has a complicated structure in phase space, then its properties and behavior cannot be well characterized in terms of just the elements of the 2 × 2 (transverse horizontal and vertical) covariance matrices. To support the argument that the failure to obtain consistent results using the three-screen analysis technique in CLARA FE was due to the structure of the beam distribution in phase space, we performed tracking simulations using an initial distribution based on that obtained from a tomography analysis (presented in Sec. III C). In the simulations, particles were tracked in a computer model of the beamline from SCR-02 (the reconstruction point) backward to SCR-01, and forward to SCR-03. At each screen, the horizontal and vertical rms beam sizes are calculated, and the same procedure that was applied to the experimental data was used to calculate the covariance matrix at SCR-02. The elements of the covariance matrix were then used to calculate the emittance and optical parameters at this point. The tracking and optical calculations were repeated for different strengths of the quadrupoles, corresponding to those used in the experiment. Results for the transverse vertical plane are shown in Fig. 4. For a Gaussian elliptical distribution in phase space, there are only very small variations in the calculated covariance matrix at SCR-02 and in the optical functions, for different quadrupole strengths (the small variations arise from statistical variation in the distribution, resulting from tracking a finite number of particles). The simulation can be repeated, but using a phase space distribution without elliptical symmetry, instead of a Gaussian elliptical distribution. Since an appropriate distribution is provided by the tomography analysis that we present later, in Sec. III C, we use this distribution (illustrated in Fig. 7) in the simulation for the nonelliptical case. However, we emphasize that the purpose at this point is only to illustrate the impact of the lack of elliptical symmetry on the threescreen analysis, rather than to demonstrate any specific aspects of the tomography analysis. Using the nonelliptical distribution in the simulations, we see much larger variations in the covariance matrix at SCR-02 and in the emittance and optical functions at this point, depending on the strengths of the quadrupoles between SCR-02 and The horizontal axis on each plot shows the phase advance between SCR-02 and SCR-03; each point represents the results for a single set of quadrupole strengths between these two screens. The plots show (from top to bottom): the mean-square beam size observed at SCR-03 scaled by the beta function; the normalized emittance; the Courant-Snyder beta function at SCR-02; the Courant-Snyder alpha function at SCR-02. The emittance, beta function, and alpha function are calculated from the covariance matrix at SCR-02 using (27) and (28). The beta function at SCR-03 (used for scaling the beam size) is calculated by propagating the Courant-Snyder functions from SCR-02 using the appropriate transfer matrix. Points leading to imaginary values for the emittance are omitted. Error ranges show the standard deviation across the measurements, omitting points with large deviation from the mean. SCR-03. For some quadrupole strengths, the calculated covariance matrix is unphysical, and it is not possible to find real values for the emittance or optical functions. The overall behavior is qualitatively similar in some respects with that seen in the experiment, Fig. 3(b). Results of simulations for the horizontal plane show the same behavior as the vertical plane, with almost no variation in the emittance or optical functions as a function of quadrupole strength for a Gaussian elliptical phase space distribution, but large variations in the case of a more realistic phase space distribution based on the results of the tomography analysis.

B. Quadrupole scan method
One of the limitations of the three-screen analysis method described in Sec. III A is the inability to provide information on beam coupling. This can be overcome, however, by combining observations of the transverse beam size at different screens for various strengths of the quadrupoles between the screens. If a sufficient number of quadrupole strengths are used, then beam size measurements at a single screen provide sufficient data to calculate the 4 × 4 transverse beam covariance matrix at a point upstream of the quadrupoles. The 4 × 4 covariance matrix has ten independent elements: in principle, just four sets of quadrupole strengths provide twelve beam size measurements (values for hx 2 i, hy 2 i and hxyi for each set of quadrupole strengths), and are more than sufficient to determine the covariance matrix. In practice, it is desirable to use a greater number of quadrupole strength settings, to overconstrain the covariance matrix. Although it is again assumed (implicitly) that the beam can be described by a simple elliptical phase space distribution, the quadrupole scan technique differs from the three-screen measurement technique in using larger number of measurements, which allows account to be taken of coupling. Also, depending on the details of the beam distribution in phase space, by overconstraining the covariance matrix, it may be possible to construct a phase space ellipse that better represents the beam behavior.
The quadrupole scan technique that we use is similar to that presented by Prat and Aiba [22]. The theory can be developed as follows. The covariance matrix Σ 3 at a location s 3 in the beamline (SCR-03 in the case of CLARA FE) is related to the covariance matrix Σ 2 at a location s 2 (SCR-02 in CLARA FE) through Eq. (26), where all matrices are now 4 × 4. The relationship between the observable quantities at s 3 (assuming a YAG screen at that location) and the independent elements of Σ 2 can be written: 0 where hx 2 3 i ðnÞ represents the mean square horizontal transverse beam size measured at s 3 for a particular set of quadrupole strengths, and similarly for hy 2 3 i ðnÞ and hx 3 y 3 i ðnÞ . With measurements of the beam distribution in coordinate space at s 3 for N different sets of quadrupole strengths, D is a 3N × 10 matrix. The elements of D can be expressed, using Eq. (26), in terms of elements of the transfer matrices M 23 from s 2 to s 3 (with each set of three rows in D corresponding to a single set of quadrupole strengths). Explicit expressions for the elements of D (for a given transfer matrix) are as follows: where m ij is the ði; jÞ element of the transfer matrix M 23 (for a given set of quadrupole strengths). Given observations of the beam profile at s 3 for a number of different sets of quadrupole strengths, and the corresponding values for the elements of D, the elements of the covariance matrix at s 2 may be found by inverting Eq. (30). Since D is not a square matrix, the pseudoinverse of D (found, for example, using singular value decomposition) must be used instead of the strict inverse. It is worth noting that whereas in one degree of freedom it is possible to obtain the elements of the covariance matrix at the reconstruction point by varying the strength of a single quadrupole between the Reconstruction Point and the observation point, this is not the case in two degrees of freedom. To understand the reason for this, consider the case of a single thin quadrupole with the reconstruction point s 2 at the upstream (entrance) face of the quadrupole, and the observation point s 3 some distance downstream from the quadrupole. The elements of the covariance matrix hx 2 3 i, hx 3 y 3 i, and hy 2 3 i each have a quadratic dependence on the quadrupole strength, with coefficients determined by the elements of the covariance matrix at the reconstruction point. By fitting the quadratic curves obtained from a quadrupole scan we therefore obtain nine constraints (three for each of the observed elements of the covariance matrix at s 2 ); however, the covariance matrix at s 2 has ten independent elements (in two degrees of freedom). The problem is therefore underconstrained: in the context of Eq. (30) this is manifest as the matrix D having fewer nonzero singular values than are required to determine uniquely the elements of the covariance matrix at the observation point. Although it is always possible to "invert" D using singular value decomposition, the procedure in this case would yield a solution for the covariance matrix that minimizes the sum of the squares of the matrix elements: there is no reason to suppose that this least-squares matrix is near the correct solution. To address this problem, however, it is only necessary to collect data from a scan of two quadrupoles at different locations between the observation point and the reconstruction point. This breaks the degeneracy in the system, and (if the system is properly designed) more than ten singular values of D will be nonzero: in other words, the system becomes overconstrained, rather than underconstrained.
The same data collected for the three-screen method can be used in the analysis using the quadrupole scan method, and the same practical considerations (concerning, for example, the desirability of a beam waist at the reconstruction point, and maintaining an approximately round beam at the observation point) apply. However, it should be noted that for the three-screen method, the observed beam sizes at all three screens are used to reconstruct the covariance matrix: an independent reconstruction is obtained for each point in the quadrupole scan. For the quadrupole scan analysis method, on the other hand, we use only the observed beam size at a single screen (SCR-03 in this case) and combine all the measurements for different quadrupole strengths to calculate the elements of the covariance matrix. In effect, we calculate the size and shape of the distribution in phase space based on the widths of projections at many different phase angles: this leads to a more reliable result than is obtained using the three-screen analysis method, for which only three different phase angles are used. Nevertheless, even for a large number of phase angles, the quadrupole scan method does not provide the same detailed information on the phase space distribution that is provided by the tomography method (discussed in Sec. III C). Rather, it attempts to fit a phase space distribution that may have significant detailed structure with a simple Gaussian elliptical distribution. Figure 5 shows the residuals from a fit based on data from a quadrupole scan in CLARA FE made with nominal machine settings. Each point indicates the observed and fitted beam size (hx 2 i, hy 2 i or hxyi) at the observation point for a different set of strengths of the quadrupoles between SCR-02 (the reconstruction point) and SCR-03 (the observation point). The results may also be validated by comparing the beam size predicted at the reconstruction point with the actual beam size observed at this point. For the case shown, the agreement is within about 15%.
To estimate the uncertainty in the elements of the reconstructed covariance matrix, we use the residuals between beam sizes observed on SCR-03 (for a given set of quadrupole strengths) and the beam sizes predicted from the model, using the reconstructed phase space distribution at SCR-02: we treat the residuals as an error on the measured beam sizes. We then construct an ensemble consisting of sets of beam sizes at SCR-03 produced by simulating the quadrupole scan procedure. Within each set, the beam size for given quadrupole strengths has a value chosen randomly from a normal distribution with mean equal to the actual measured beam size for those quadrupole strengths, and standard deviation equal to the corresponding residual. From each simulated quadrupole scan (i.e., for each member of the ensemble) we find the corresponding emittance and optics functions. Finally, we assume that the uncertainty on these values can be found from the standard deviation of the values across the ensemble. The uncertainties found in this way are shown as the range of variation in the results from the quadrupole scans (in one and two degrees of freedom) for the emittances and optics functions in Table I.

C. Phase space tomography
Finally, it is possible to use phase space tomography to construct a more detailed representation of the beam properties than is provided by just the emittance and optical functions. In principle, the tomography method is similar to the quadrupole scan, in that by observing the beam image on a screen for different strengths of a set of upstream quadrupoles, it is possible to reconstruct the phase space distribution at a point upstream of the quadrupoles. The difference is that for the quadrupole scan, only the rms beam sizes are used in the analysis: tomography uses all the information from the (observed) beam distribution to produce a more detailed reconstruction of the phase space distribution of the beam. When tomography is carried out in coordinate space, two-dimensional images (projections) on a screen for different orientations of an object are used to reconstruct a three-dimensional representation of the object. In phase space tomography, different "orientations" correspond to rotations in phase space, which are achieved by changing the horizontal or vertical phase advances between the reconstruction point and the observation point. Mathematically, the analysis is essentially the same as in coordinate space, and standard algorithms developed for tomography in coordinate space (such as filtered  back-projection, or maximum entropy [27][28][29]) can be applied to phase space tomography. For analysis of data from CLARA FE, we have used a form of algebraic reconstruction. The procedure may be outlined as follows. For simplicity we consider just a single degree of freedom: the generalization to two (or more) degrees of freedom is straightforward. Let ψ be a vector in which each element ψ j represents the beam density at a particular point ðx j ; p xj Þ in phase space at the reconstruction point. Assuming that the points are evenly distributed on a grid in phase space, then the projected beam density at the observation point, ρ ð1Þ i at a point x ¼ x i in coordinate space, can be written as a matrix multiplication: where the matrix P ð1Þ has elements: 12 are elements of the transfer matrix from the reconstruction point to the observation point, for a given set of quadrupole strengths. If the vector ρ has N elements ρ i , and there are N 2 points ðx j ; p xj Þ in phase space, then P ð1Þ is an N × N 2 matrix. If the transfer matrix from the reconstruction point to the observation point is changed (e.g., by changing the strengths of the quadrupoles between the two points), then we construct a new vector ρ ð2Þ from the new image at the observation point, corresponding to the new transfer matrix. In general, for the nth transfer matrix, we have: Note that the phase space density ψ is constant, because ψ refers to a point upstream of any quadrupoles whose strength is changed during the measurements. We can combine the observations simply by stacking the vectors ρ ðnÞ and the matrices P ðnÞ : . . .
ρ is a vector with nN elements, and P is an nN × N 2 matrix. In terms of the pseudoinverse P † of P, we have the following formula for the phase space density at the reconstruction point: We perform the analysis in normalized phase space [16], in which the phase space variables ðx N ; p xN Þ are defined by: where α x and β x are the Courant-Snyder functions at the given point in the beam line. The transfer matrix in normalized phase space between any two points in the beam line is represented by a pure rotation matrix, with rotation angle given by the phase advance. This simplifies the implementation of the algebraic tomography method described above. A further advantage of working in normalized phase space is that if the Courant-Snyder functions at the reconstruction point are chosen to match the beam distribution, then the beam distribution in phase space at this point will be perfectly circular: this improves the accuracy with which parameters such as the emittance may be calculated. Note that, since we do not know in advance the actual Courant-Snyder parameters describing the beam distribution at the reconstruction point, we need to make some estimate based on (for example) simulations or a quadrupole scan analysis. In practice, it is not essential for the estimated parameters to match exactly the actual beam parameters: any discrepancy will simply lead to an elliptical distortion of the beam distribution in normalized phase space. To transform experimental observations into normalized co-ordinates is straightforward: all that is necessary is to scale the coordinate axis for the observed beam projection by a factor 1= ffiffiffiffiffiffiffi β OP x p , where β OP x is the Courant-Snyder beta function at the observation point calculated from the estimated (fixed) Courant-Snyder functions at the reconstruction point and the transfer matrix from the reconstruction point to the observation point.
Rather than compute the pseudoinverse of P, we solve Eq. (34) iteratively, using a least-squares method. For the computation of the phase space in a single degree of freedom, we apply a constraint that the particle density must be positive at all points in phase space. However, applying this constraint carries a large computational overhead, and for computation of the phase space in two degrees of freedom, which has considerably greater computational cost than the case of a single degree of freedom, we do not constrain the least-squares solver in this way. This can result in negative (unphysical) values for the particle density at some points in phase space; however, when a good fit is achieved, the negative values make a relatively small contribution to the overall phase space distribution.
Although there is no need for the phase advances between the reconstruction point and observation point to be evenly distributed over the set of observations for different quadrupole strengths, it generally improves the accuracy of the tomography analysis to use as wide a range of phase advances as possible, with roughly uniform spacing: this maximizes the overall constraints on the phase space distribution for a given number of observations. The sets of quadrupole strengths identified in the preparatory simulations (described above) were specifically chosen to provide a wide range of phase advances. The same data (screen images at the observation point, for a range of different quadrupole strengths) can be used for the three-screen analysis (described in Sec. III A), the quadrupole scan analysis (described in Sec. III B) and the tomography analysis described here. Figure 6 shows a set of results from tomography analysis for the nominal machine settings, and in which the horizontal and vertical phase spaces are treated independently. As was the case for the quadrupole scan method, the results may be validated by comparing the predicted beam size at the reconstruction point with the beam observed directly at this point (SCR-02): the results of the comparison are shown in the lower plots in Fig. 6. There is good agreement, and it can be clearly seen that the tomography analysis reveals some features of the charge distribution in phase space that are not obtained from the quadrupole scan analysis. For example, in the horizontal phase space, we see three distinct "peaks" in the charge density, which are associated with peaks in the variation of the intensity observed on SCR-02, projected onto the horizontal axis. The quadrupole scan analysis provides only the parameters of a Gaussian elliptical distribution in phase space, which would lead to a simple Gaussian variation in the intensity projected onto the horizontal axis.
Although the phase space distributions are not perfectly elliptical, we can obtain indicative values for the emittances and Courant-Snyder parameters using a number of different methods. For example, given a phase space distribution, it is possible to calculate the second-order moments, and then to find the emittance and Courant-Snyder parameters using (in one degree of freedom) Eqs. (2)- (5). However, this approach may not give useful results (in terms of predicting beam behavior) if the density does not fall to zero rapidly with distance from the center of the beam. Although this may be addressed by imposing a "cutoff", where the density beyond some amplitude is set to zero, it is not always clear where such a cutoff should be imposed. Also, it is not clear how to estimate uncertainties (or errors) on the values obtained.
An alternative approach is to fit an elliptical distribution function to the phase space density. For example, an elliptical Gaussian of the form (9) could be used. Using a standard algorithm, such as nonlinear least-squares regression, it is then also possible to estimate the uncertainty on values obtained for the emittance and Courant-Snyder parameters. This is the approach used to obtain the values presented in Sec. IVA (Table I). A Gaussian elliptical distribution in horizontal phase space can be written: where ρ 0 is the peak density, ⃗ x T ¼ ðx; p x Þ is a phase-space vector, and the 2 × 2 symmetric matrix Σ −1 is the inverse of the covariance matrix. The peak density ρ 0 and the three independent components of Σ −1 are used as variables in fitting the Gaussian elliptical function to the phase space distribution. We use the FITNLM function in MATLAB [30], which performs nonlinear least squares regression using the Levenberg-Marquardt algorithm [31]. This provides values for the variables in the fit that give the best match (in terms of minimizing the squares of the residuals) to the given data, and the standard error on each variable. The fitting procedure also provides an indication of the quality of the fit in terms of the coefficient of determination, or r 2 value (in the range 0 to 1, with a value of 1 indicating a perfect match of the model to the data). Values for the emittance and optics functions can be obtained from the elements of the covariance matrix Σ describing the fitted Gaussian. In practice, we perform the fit in normalized phase space, and apply the appropriate transformation to find the covariance matrix for the Gaussian elliptical distribution that matches the given distribution in ordinary (not normalized) phase space.
To estimate the uncertainty in the values obtained for the emittances and optics functions, we construct an ensemble of matrices fΣ −1 n g, where, for each member Σ −1 n of the ensemble, the value of each element is chosen randomly from a normal distribution with mean equal to the corresponding fitted value of Σ −1 , and standard deviation equal to the corresponding standard error. From the ensemble fΣ −1 n g we construct ensembles of values for the emittances and optics functions: in Table I, the values obtained using the tomography techniques are the mean of each ensemble, with uncertainty given by the standard deviation across the ensemble.
Treating the horizontal and vertical phase spaces separately in the analysis means that no information is provided on coupling in the beam, which may arise (for example) from incorrect setting of the bucking coil at the electron source. It is possible to extend the tomography analysis from a single degree of freedom, to treating two degrees of Projections of beam density in normalized phase space, found from phase space tomography in two degrees of freedom in CLARA FE. Each plot shows a different projection of the charge density from four-dimensional phase space, using normalized phase space variables. The black ellipses show projections of the four-dimensional emittance ellipse obtained from a Gaussian fitted to the four-dimensional (normalized) phase space distribution. Coupling in the beam is evident in the tilt of the charge distribution in the cases that the axes refer to different degrees of freedom. The left-hand set of plots (a) shows the phase space distribution reconstructed from experimental data; the right-hand set of plots (b) shows the results of the tomography analysis applied to simulated data based on the measured phase space distribution, to validate the technique. Although there are some differences between the analysis results from the experimental data and the results from the simulated data, there is overall very good agreement in the phase space distribution found in each case, and in the emittances and optical functions corresponding to fitted emittance ellipses (see Table I).
freedom simultaneously [14]. Applying this technique to the case considered here, the resulting four-dimensional phase space reconstruction includes information about the betatron coupling in the beam. Some results from experimental data (screen images) are shown in Fig. 7(a). Generally, the fit using Eq. (36) of the phase space beam density to the observed images is very good: the residuals from a typical example are shown in Fig. 8. The beam emittances and optics functions in two degrees of freedom can be found from the reconstructed phase space using a generalization of the method described above for a single degree of freedom, based on fitting a Gaussian to the reconstructed phase space distribution. In two degrees of freedom, the Gaussian function is still given by Eq. (38), but the phase space vector is now ⃗ x T ¼ ðx; p x ; y; p y Þ, and the matrix Σ −1 is now a 4 × 4 symmetric matrix with ten independent elements that are used as variables in the fit. Values for the normal mode emittances and optics parameters can be found from the eigenvalues and eigenvectors (respectively) of the covariance matrix Σ, as described in [19]. Uncertainties in these values can again be obtained from an ensemble of matrices fΣ −1 g, constructed using the standard errors on the variables used in the fit. One drawback of applying phase space tomography in two degrees of freedom is that the matrix P in Eq. (32) becomes very large: in one degree of freedom, to reconstruct the phase space distribution with resolution N in each dimension using n observations, P will be an nN × N 2 matrix. In two degrees of freedom (four-dimensional phase space), P will be an nN 2 × N 4 matrix: even for a relatively coarse reconstruction, with N of order 50, computing P and applying its inverse can require significant computational resources. The situation is eased somewhat by the fact that in practice, P is a sparse matrix, and this allows a significant reduction in the computer memory that would otherwise be required; nevertheless, the required computational resources can be a limit on the resolution with which the phase space in two degrees of freedom may be reconstructed. The results shown here use a four-dimensional phase space resolution N ¼ 69.
Projections from the four-dimensional phase space density found from experimental data in CLARA FE are shown in Fig. 7(a). To validate the technique, we take the four-dimensional phase space distribution, and use it in a simulation to construct a set of images on SCR-03 corresponding to different quadrupole strengths. We then take the simulated images, and again apply the tomography analysis: the results are shown in Fig. 7(b). Although there are some differences between the original and reconstructed distributions they are sufficiently close to indicate that the technique potentially has good accuracy. We also find that there is good agreement between the emittances and optics functions obtained by fitting ellipses to the projections of the phase space into the horizontal and vertical planes (see Table I).
We can further validate the results by reconstructing the two-dimensional distribution in coordinate space at the reconstruction point (by projecting the four-dimensional phase space distribution onto the coordinate axes), and comparing this with the image that is observed directly. Some examples for such comparisons are shown in Fig. 9. In general, we find that the images reconstructed from phase space tomography reproduce reasonably well the general shape and some of the more detailed features of the images that are observed directly. However, the tomography does not reveal the same level of detail as can be seen in the observed image. This may be due in part to the limited resolution of the tomography analysis: as already mentioned, for the analysis presented here, we used a phase space resolution of 69 points on each of the four axes (which was at the upper limit set by the available computer memory). However, it is also likely that measurement errors also play a role. We note that the residuals of the fits to the images at the observation point are typically very small (so that there is no discernible difference between the directly observed images at this point and the images reconstructed from phase space tomography: see the example in Fig. 8). However, there are systematic differences between the reconstructed images at SCR-02 (the reconstruction point), and the images observed directly on that screen. In particular, the vertical size of the reconstructed beam (projecting the phase space distribution onto the vertical axis) is generally of order 10% larger than the vertical size of the image observed directly. Work is in progress to understand and correct the systematic errors: possible sources include calibration errors in the quadrupoles and in the diagnostics used for collecting beam images. It is important to have accurate values for the quadrupole Similarly, the analysis depends on accurate knowledge of the calibration factors of the diagnostic screens. Hysteresis in the quadrupole magnets used to change the optics between the reconstruction point and the observation point may also lead to errors in the analysis: to try to minimize hysteresis effects, the quadrupoles were routinely degaussed (cycled) between scans, but the time taken for this procedure made it impractical to degauss the quadrupoles at each point in a single scan.

IV. EMITTANCE AND OPTICS MEASUREMENTS UNDER VARIOUS MACHINE CONDITIONS
A. Nominal machine settings Table I shows the emittance and optics parameters obtained under nominal machine settings using the three different techniques discussed in the previous sections: three-screen measurements, quadrupole scans, and phase space tomography. For the quadrupole scan and tomography analysis in four-dimensional phase space, the emittance and optics values in the table are those for the normal mode quantities as described in Sec. II. With the nominal machine settings, the electron source and linac operate with the beam on-crest of the rf (i.e., to give maximum beam acceleration for a given rf amplitude), with amplitudes producing beam momentum 5 MeV=c and 30 MeV=c respectively. The current in the bucking coil is set to cancel the solenoid field on the photocathode, and the laser intensity is set to give a bunch charge of 10 pC. The results in Table I are based on the same data set (i.e., the same set of beam images) in each case; the only difference between the different methods is in the way that the data are analysed. In principle, therefore, for a beam with elliptical symmetry in the phase space distribution, we would expect to see close agreement between the values obtained using different techniques. However, while there is reasonable agreement in some of the cases for the data from CLARA FE, there is also wide variation in the values for some parameters (e.g., the horizontal normalized emittance in two-dimensional phase space). As discussed in Sec. III A, this is likely the result of the complicated structure of the beam distribution in phase space.
For the three-screen and quadrupole scan techniques, the uncertainties in the values shown in Table I provide an indication of the extent to which the given values describe the beam behavior: in effect, the uncertainties in these cases indicate the quality of a fit of a phase space distribution with elliptical symmetry to the actual phase space distribution. For the tomography technique, however, the uncertainties are based on the standard error on the parameters describing a Gaussian fitted to the reconstructed distribution. The error on the parameters may be small, even if the quality of the fit is poor. A better indication of the quality of the fit in this case is given by the coefficient of determination, or r 2 . For the tomography analysis in one degree of freedom, we find for the case shown in Table I that r 2 ¼ 0.60 for the horizontal phase space, and r 2 ¼ 0.80 for the vertical phase spaces. In two degrees of freedom, we find r 2 ¼ 0.10. This shows a poor fit to a Gaussian elliptical distribution in the four-dimensional phase space, but rather better fits to the distributions projected onto the horizontal and vertical (two-dimensional) phase spaces.

B. Effect of varying bucking coil strength
The electron source in CLARA FE is constructed so that the field from the solenoid can be cancelled at the cathode by the field from a bucking coil. If the current in the bucking coil is changed from the value needed to achieve cancellation, electrons are emitted from the surface of the cathode in a nonzero solenoid field: the effect is to introduce some coupling into the beam (as a result of noncompensated azimuthal momentum), which can appear as changes in the beam emittances. In particular, the individual normal mode emittances will vary, though their product should remain constant as a function of the solenoid field strength on the cathode [32,33]. The difference between the normal mode emittances is expected to be minimized when there is zero solenoid field on the cathode: with increasing field strength (parallel to the longitudinal axis, in either direction) one emittance will increase while the other will decrease. Tuning the machine for optimum performance generally involves minimizing the coupling, to achieve the smallest possible emittance ratio [15], and characterizing and understanding the coupling as a function of the strength of the bucking coil is thus an important step in machine commissioning. Four-dimensional phase space tomography offers a powerful tool for providing insight into coupling in the machine, and was used to study the dependence of the phase space distribution on the current in the bucking coil.
The normal mode emittances as a function of current in the bucking coil, found from four-dimensional phase space tomography (as described in Sec. III B) are shown in Fig. 10. Although there is some variation in the product of the emittances with changes in the current in the bucking coil, over a wide range the variation is small. There is also some indication of the expected behavior of the individual emittances. The difference between the emittances is minimized for a bucking coil current of approximately −3.5 A: this is different from the nominal value of −2.2 A for cancelling the field on the cathode. The reason for the discrepancy is being investigated. Note that before collecting data over the range of bucking coil currents, the bucking coil was degaussed with the intention of improving the agreement between the cathode field calculated from a computer model of the electron source and the field that was actually produced for a given current. It is also worth noting that the time taken for data collection over the full range of bucking coil currents took several hours, and it is likely that some variation in machine parameters (such as rf phase and amplitude in the electron source and linac) occurred over this time.
Also shown in Fig. 10 are results from a GPT simulation and from a simple theoretical model: these are included in the figure to illustrate the expected behavior of the normal mode emittances as a function of the solenoid field on the cathode, and are not intended to show results from an accurate machine model (although we see that they match the results from the GPT simulations very well). For the simulations, we use parameters for the electron source corresponding to those in CLARA FE, but with the field from the bucking coil scaled to cancel the solenoid field on the cathode for a current of −3.5 A in the bucking coil (rather than the nominal −2.2 A). Also, the initial distribution of particles in phase space is chosen to give emittances (with zero solenoid field on the cathode) corresponding to the experimental measurements. This requires the beam divergence at the cathode to be scaled to exceed significantly the values believed to be appropriate for CLARA FE; however, it should be remembered that in the simulation, the emittances are calculated immediately after the electron source, whereas the measurements are made in a section of beamline downstream of the linac and numerous other components. Effects (that are not yet well characterized) between the electron source and the measurement section are likely to lead to some increase in emittance. The GPT and theoretical results are therefore included in Fig. 10 purely to illustrate the expected behavior of the emittances as functions of the strength of the solenoid field on the cathode, rather than as a direct comparison of an accurate computational model with the experimental results.
Also shown in Fig. 10 are results from a simplified theoretical (analytical) model. This is based on an assumed beam phase space distribution at the cathode, i.e., immediately after photoemission. If there is no magnetic field on the cathode, then the covariance matrix is characterized by an emittance and beta function in each transverse direction: A solenoid field of strength B 0 on the cathode can be represented by a vector potential: so that the canonical conjugate momenta p x and p y become: where m and e are the mass and magnitude of the charge of the electron, v x and v y are the transverse horizontal and vertical components of the velocity, and P 0 ¼ β 0 γ 0 mc is the reference momentum (which can be chosen arbitrarily). The covariance matrix then becomes: where: Finally, from the covariance matrix (43), we find (using the methods described in Sec. II) that the normal mode emittances are given by: where: The normalized emittances (ϵ N;I ¼ β 0 γ 0 ϵ I , and similarly for ϵ N;II ) remain constant during acceleration of particles in the rf field of the electron source (and in the linac). To apply this model to CLARA FE, giving the results shown in Fig. 10, the initial beam size and divergence are chosen to fit the emittances at their closest approach: the values used are close to those used in the GPT simulation. We also assume that η ¼ 0 for a bucking coil current of −3.5 A, and scale the dependence of η on the field in the bucking coil so as to match the experimental curves. However, we again emphasize that the results from the theoretical model and the GPT simulation are included only to give an illustration of the expected behavior, and are not directly comparable with the experimental results. Direct inspection of the phase space distribution provides a further indication of how the coupling changes with the current in the bucking coil. For example, Fig. 11 shows the projection onto the x-p y plane of the four-dimensional phase space (reconstructed from the tomography measurements) for different values of the current in the bucking coil. The "tilt" on the distribution corresponds to a correlation between the horizontal coordinate and vertical momentum, and indicates the coupling: we see that this changes sign as the bucking coil current is varied from −5 A to −1.5 A. The tilt (and hence the coupling) vanishes for a current of approximately −3.5 A, which is consistent with the current required to minimize the difference between the normal mode emittances.
A more complete characterization of the coupling is given in Fig. 12, which shows the elements of the covariance matrix at SCR-02 (the reconstruction point) as functions of current in the bucking coil. Coupling between motion in the horizontal and vertical directions is indicated by nonzero values of the elements in the 2 × 2 top-right block diagonal. All these elements vanish for bucking coil currents close to −3.5 A. We do not expect this to correspond exactly to the bucking coil current that minimizes the separation between the normal mode emittances, since after leaving the cathode (in zero longitudinal magnetic field) the particles then pass through a section of main solenoid field, not cancelled by the bucking coil. The main solenoid field introduces some coupling in the beam, characterized by nonzero elements off the 2 × 2 block diagonals in the covariance matrix. However, tracking simulations in GPT suggest that in the case of CLARA FE, the coupling in the covariance matrix introduced by the part of the main solenoid not cancelled by the bucking coil is small: the coupling in the covariance matrix is minimized at a current within about 0.1 A of the current that gives the closest approach of the normal mode emittances (see Fig. 10).

C. Effects of varying main solenoid strength and bunch charge
Although space-charge effects in CLARA FE are negligible in the section of beamline where the emittance and optics measurements are made (with beam momentum around 30 MeV=c), space-charge forces can play a significant role in the electron source, depending on the bunch length and the total bunch charge. In the studies reported here, the photocathode laser was operated with pulse length of 2 ps: as discussed in Sec. III (where we considered the beam perveance) space-charge effects are expected to be weak up to bunch charges of around 50 pC. However, screen images suggested some significant variation in beam parameters even at lower bunch charges. It is planned in the future to use phase space tomography for rigorous studies of the impact of bunch charge (and other parameters) on beam properties; but so far, the limited time available for collecting quadrupole scan data, together with some variability in the machine conditions, has made it impractical to make detailed, systematic measurements. Nevertheless, to provide some information on beam behavior, quadrupole scans were performed for bunch charges of 10 pC, 20 pC, and 50 pC, and for a reduced main solenoid current of 125 A, as well as the nominal 150 A. Some of the results from analysis of these quadrupole scans using phase space tomography in two degrees of freedom are shown in Fig. 9, which compares the reconstructed beam image at SCR-02 with the image observed directly on this screen.
The images in Fig. 9 indicate significant detailed structure in the beam distribution, depending on main solenoid current, bucking coil current, and bunch charge. This is apparent both from the image observed directly at the reconstruction point, and from the phase space Note that in the case of the tomography analysis, the elements describing coupling between the transverse degrees of freedom (shown in the top right 2 × 2 block diagonal) vanish for a bucking coil current of approximately −3.5 A: this is consistent with the current at which the difference between the normal mode emittances is minimized (as shown in Fig. 10). distribution constructed from four-dimensional tomography. In such cases, the phase space cannot accurately be characterized simply by the emittances and optical functions that describe the covariance matrix. Nevertheless, to allow some comparison, we calculate the normal mode emittances and optical functions, using the method described in Sec. II: the values of the normal mode emittances and selected optical functions are shown in Table II. Although there are indications of some patterns (for example, an increase in emittance with bunch charge) no firm conclusions can be drawn because machine conditions between different quadrupole scans were not accurately reproducible. Nevertheless, the measurements that have been made demonstrate the potential value of four-dimensional phase space tomography for developing an understanding of the beam physics in a machine such as CLARA FE, and for tuning the machine for optimum performance.

V. SUMMARY AND CONCLUSIONS
We have presented the first experimental results from four-dimensional phase space tomography in an accelerator. The beam emittance and optical properties obtained from phase space tomography have been compared with results obtained using more commonly employed techniques, such as three-screen analysis and quadrupole scans. We considered the suitability of each method for situations where the beam contains detailed structures in phase space and cannot be described by a simple elliptical distribution. In this case, the three-screen method can give inconsistent and often unphysical results. By contrast, the quadrupole scan method provided approximate values for the emittances which appear broadly to agree with those obtained using tomography (see Table I). However, a proper description of a nonelliptical phase space distribution cannot be given just in terms of a small number of parameters (emittance and Courant-Snyder parameters). Phase space tomography overcomes this limitation by providing the beam density at a number of points in phase space.
Our results for the phase space tomography analysis and the comparisons with other methods are supported by simulation studies. The results of the tomography analysis have been validated by comparing (for example) the beam image at the entrance of the measurement section of the beamline (the reconstruction point) obtained from a projection of the measured four-dimensional phase space, with the beam image observed directly on a screen at this point. In general, the agreement suggests that four-dimensional phase space tomography is providing a useful representation of the beam properties, though the image reconstructed from tomography lacks the same resolution as the image observed directly. There is also evidence for systematic errors in the measurement that need to be properly understood.
A benefit of four-dimensional phase space tomography (compared to tomography in two-dimensional phase space) is that the technique provides detailed information on coupling in the beam. This can be important for tuning a machine such as CLARA FE, for example, where solenoids are used to provide focusing for the beam, but it is desirable to minimize the coupling that can be introduced by those solenoids. Information on coupling can be obtained by applying the quadrupole scan method in two (transverse) degrees of freedom; but information obtained in this way may not be accurate or reliable if there is detailed structure in the beam distribution.
The main drawback of the tomography analysis is that collection of the data may be a time-consuming procedure. In cases where the beam distribution in phase space is smooth and without significant detailed structure (so that it can be well characterized by the emittance and optical functions) then the three-screen or quadrupole scan techniques, using a limited set of observations, may provide sufficient information for machine tuning and optimization relatively quickly. Phase space tomography generally requires data from a larger number of observations, but depending on the level of detail or accuracy required, it may be possible to minimize the number of points in the quadrupole scan used to provide the data: the limits of the technique have still to be rigorously explored, and will likely depend on the specific machine to which it is applied. Regarding practical application of phase space tomography, it is worth mentioning that the requirements in terms of beamline design and diagnostics capability are not demanding. In CLARA FE, the diagnostics section consists of a short (1.661 m) section of beamline between two transverse beam profile monitors, and containing three (adjustable strength) quadrupoles. The design of this section was developed before detailed plans were prepared for phase space tomography studies, and there is limited flexibility in optimizing the phase advances and optical functions over the length of the diagnostics section. Nevertheless, it was possible to identify sets of quadrupole strengths to provide the observations necessary for the analysis and results presented here.
So far, we have used an algebraic reconstruction technique for the phase space tomography. This technique has the advantage (compared to other tomography algorithms) of ease of implementation and flexibility in terms of the input data. However, it is possible that different algorithms may provide better (more accurate, or more detailed) results, and we hope to explore the possible benefits and limitations of alternative tomography methods. A particular issue with tomography in four-dimensional phase space is the demand on computer memory for processing the data and storing the results, especially at high resolution in phase space. However, because of the nature of the problem, the memory requirements will almost inevitably scale with the fourth power of the phase space resolution, and it seems unlikely that other tomography methods would provide significant benefit in this respect. It is possible that more sophisticated computational techniques may allow some reduction in the memory requirements for a given resolution, e.g., [34].
While improvements and refinements in the technique are planned, the results so far show that four-dimensional phase space tomography is a useful technique for detailed beam characterization and for machine tuning and optimization. It is hoped that further studies will include investigation of space-charge effects in the electron source and beam dynamics effects (such as wake fields) in the linac.

ACKNOWLEDGMENTS
We would like to thank our colleagues in STFC/ASTeC at Daresbury Laboratory for help and support with various aspects of the simulation and experimental studies of CLARA FE. This work was supported by the Science and Technology Facilities Council, UK, through a grant to the Cockcroft Institute.