A Geometrical Method of Decoupling

The computation of tunes and matched beam distributions are essential steps in the analysis of circular accelerators. If certain symmetries - like midplane symmetrie - are present, then it is possible to treat the betatron motion in the horizontal, the vertical plane and (under certain circumstances) the longitudinal motion separately using the well-known Courant-Snyder theory, or to apply transformations that have been described previously as for instance the method of Teng and Edwards. In a preceeding paper it has been shown that this method requires a modification for the treatment of isochronous cyclotrons with non-negligible space charge forces. Unfortunately the modification was numerically not as stable as desired and it was still unclear, if the extension would work for all thinkable cases. Hence a systematic derivation of a more general treatment seemed advisable. In a second paper the author suggested the use of real Dirac matrices as basic tools to coupled linear optics and gave a straightforward recipe to decouple positive definite Hamiltonians with imaginary eigenvalues. In this article this method is generalized and simplified in order to formulate a straightforward method to decouple Hamiltonian matrices with eigenvalues on the real and the imaginary axis. It is shown that this algebraic decoupling is closely related to a geometric"decoupling"by the orthogonalization of the vectors $\vec E$, $\vec B$ and $\vec P$, that were introduced with the so-called"electromechanical equivalence". We present a structure-preserving block-diagonalization of symplectic or Hamiltonian matrices, respectively. When used iteratively, the decoupling algorithm can also be applied to n-dimensional systems and requires ${\cal O}(n^2)$ iterations to converge to a given precision.


I. INTRODUCTION
The significance of the symplectic groups in Hamiltonian dynamics has been emphasized for instance by A. Dragt 5 , and it has long been known 6 that the Dirac matrices are generators of the symplectic group Sp(4, R). In Ref. 4 the author presented a toolbox for the treatment of two coupled harmonic oscillators that is based on the use of the real Dirac matrices (RDMs) as generators of the symplectic group Sp(4, R) and a systematic survey of symplectic transformations in two dimensions. This toolbox enabled the developement of a straightforward recipe for the decoupling of positive definite two-dimensional harmonic oscillators. Here we present an improvement of the method that is based on geometric arguments, i.e. on the orthogonalization of 3-dimensional vectors associated via the electromechanical equivalence (EMEQ) to certain linear combinations of matrix elements.
There is a long history of publications covering the diagonalization (and related) problems in linear algebra as well as in linear coupled optics, linear Hamiltonian dynamics and control theory. A (non-exhaustive) list is given in the bibliography (see Refs. [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] , but also Ref. 3,4 a) Electronic mail: christian.baumgarten@psi.ch and references therein). However, none of the previous works (known to the author) takes full advantage of the group structure of the generators of Sp (4). The conceptually closest approach uses "quaternions", the representations of which seems to be identical to the RDMs 27 , but seems to be limited to orthogonal symplectic transformations. The decoupling method of Teng and Edwards has been the starting point for this work, as it turned out to fail in some special cases (see Ref. 3 and App. D).
The method that we present here, is based on a survey of all symplectic similarity transformations. We do not make specific assumptions about the Hamiltonian other than that it is a symmetric quadratic form and we present a geometric interpretation via the EMEQ, which provides a physical notation of otherwise complicated and nondescriptive algebraic expressions 1 . Furthermore we believe that the use of the EMEQ is an interesting example of how elements of classical physics, quantum mechanics, special relativity, electrodynamics, group theory, geometric algebra, statistics 28 and last but not least symplectic theory fit together and allow to use a common formalism.
The simplest classical linear dynamical system with interaction (coupling) has two degrees of freedom and hence a 4-dimensional phase space. It can be considered as "fundamental" and a detailed analysis of its properties will likely be instructive also for n > 2. Indeed it turns out, that the decoupling technique of a two-dimensional system can be iteratively applied to systems with more than two degrees of freedom. A Jacobi-like iteration with pivot-search and is sketched in Sec. V C.
From the viewpoint of coupled linear optics the problem is solved if a symplectic transformation is derived that transforms (constant) Hamiltonian matrices to 2×2block-diagonal form (see below). It has been shown in Ref. 4 , that the same transformation method can be applied to symplectic matrices as well. The arguments will be briefly reported below. When applied to symplectic matrices, the method is equivalent to the computation of the matrix logarithm. A solution for the couterpart, i.e. the computation of the matrix exponential with emphasis on the use of Dirac matrices, has been presented by Barut, Zeni and Laufer in 1994 7 .
If 2 × 2-block-diagonal form has been achieved, the remaining task is completely analogous to the application of the Courant-Snyder theory for one degree of freedom. Nevertheless some arguments require awareness of the eigenvalues and their relation to the properties of the Dirac matrices so that a reference to a complete diagonalization seemed appropriate.

II. COUPLED LINEAR OPTICS
The Hamiltonian of a n-dimensional harmonic oscillator with arbitrary coupling terms can be written in the form where A is a symmetric matrix and ψ is a state-vector or "spinor" of the form ψ = (q 1 , p 1 , q 2 , p 2 , . . . , q n , p n ) T . Even though the matrix A is time-dependent in the general case, it is well-known practice to use the Floquettransformation to reduce it to a constant matrix for the treatment of periodic systems (see app. B and for instance Ref. 30,33 ). The symplectic unit matrix (usually labeled J or S) is a skew-symmetric matrix that squares to the negative unit matrix. For n = 2 it is identified with the real Dirac matrix γ 0 . As described in Ref. 4 , it is possible to freely choose the order of the variables in the state vector. However, the order of the variables fixes the form of the symplectic unit matrix γ 0 2 . We prefer the use of a ordering system in which the phase space coordinates (q i , p i ) are grouped as pairs of canoni-cal conjugate variables, so that γ 0 has the form Using the over-dot to indicate the derivative with respect to time (or path length), the equations of motion (EQOM) have the familiar forṁ or in vector notation: where the force matrix F is given as From the definition of F one quickly finds that 30,33 where the superscript "T" denotes the transposed matrix. Matrices that obey Eqn. (6) are usually called "infinitesimally symplectic" or "Hamiltonian" 30 . Both terms are -in the opinion of the author -misleading: The former because F is neither symplectic nor is it infinitesimal, the latter since F does not appear in the Hamiltonian while the symmetric matrix A does. In addition A and not F is in the view of the author the classical counterpart of the Hamiltonian operator (see App. A in Ref. 4 ). Furthermore, Eqn. (6) is a purely formal property and not necessarily connected to a Hamiltonian. Therefore the author decided to use the term "symplex" (plural "symplices") when referring to its formal definition (i.e. Eqn. (6)) and its relation to the symplectic transfer matrix and to call it "force matrix" when referring to its physical content -especially with respect to the EMEQ (see Ref. 4 and below). Accordingly we speak of an antisymplex or cosymplex (i.e. "skew-Hamiltonian" matrix), if a matrix C fulfills the equation If we write S (C) for (co-)symplices, respectively, optionally with a subscript, then it is easy to prove that and A. Dirac Matrices In the following we focus on two degrees of freedom (n = 2), i.e. to a four-dimensional phase space and the use of the real Dirac matrices to describe its dynamics and transformation properties. Often the term "Dirac matrices" is used more restrictively and designates only four matrices, namely γ k , k ∈ [0 . . . 3]. Here we consider the four basic Dirac matrices as the four basic elements of a Clifford algebra Cl(3, 1) with 16 elements derived from the basic matrices (see app. A). For further details see for instance Ref. 8,31,32 .
Any real 4 × 4-matrix M can be written as a linear combination of the RDMs The RDM-coefficients m k are given by 3 where Tr(X) is the trace of the matrix X. Only the first ten RDMs are symplices and since symplices obey the superposition principle 4,5,33 , any force matrix (symplex) can be written as F = It has been shown in Ref. 4 , that the decoupling of the symplex-part M s of a symplectic matrix M automatically decouples the corresponding cosymplex M c . Hence it is sufficient to derive a method to decouple symplices of the above mentioned type. In cases where only the one-turn-transfer matrix is available, Eqn. (16) is used beforehand to extract the symplex-part of the transfer matrix. The decoupling algorithm can then be applied to this matrix (see also the detailed discussion in Ref. 4 ).

III. BLOCK-DIAGONALIZATION AND EIGENVALUES
The force matrix F is by definition a product of a symmetric matrix A and of a skew-symmetric matrix γ 0 . Hence it has zero trace and the sum of all eigenvalues is zero. We restrict ourselves to systems with realvalued force matrices and therefore real-valued transfer matrices. The eigenvalues of real-valued 2 × 2-symplices are either both real or both purely imaginary (since they are the square root of a real expression). Blockdiagonalization (in the case of the variable ordering as described above) means to find a symplectic similarity transformation R such that the matrixF = R F R −1 has the formF whereF k are real 2 × 2-matrices. Since similarity transformations preserve the eigenvalues, a symplex is blockdiagonalizable in the form that we are going to describe, if the (pairs of) eigenvalues are either real or imaginary. In case of imaginary eigenvalues, the corresponding degree of freedom (i.e. pair (q i , p i )) is stable (or focused), while a pair of real eigenvalues belongs to an unstable (non-focused) degree of freedom. The corresponding betatron motion is unstable in the sense, that no sufficient focusing is present. However -in the general coupled case without further assumptions -F is a general 4 × 4-symplex (or larger). Using the RDMs it is relatively easy to construct matrices with complex eigenvalues. An example is which has the complex eigenvalues ±i (B x ± i E x ). Since the eigenvalues are complex, also the 2 × 2-blocks are complex. They can be block-diagonalized, but the generalization to the 2 n × 2 n-case requires a general treatment of the complex case, which goes beyond the scope of this paper. As in Ref. 4 the author speaks of regular or massive systems, if the Hamiltonian is positive definite and of irregular or magnetic systems in case of indefinite Hamiltonian, respectively. Both types may be stable or unstable and this distinction should not be confused with the question of stability. A detailed discussion of stability would go beyond the scope of this paper and we refer the reader for instance to Ref. 33 or Ref. 27 and references therein.

A. The S-matrix
The matrix of second moments σ of a charged particle distribution has the time derivativė Multiplication from the left with γ 0 and the use of Eqn. (6) leads toṠ where the matrix S is defined by If Eqn. (23) is compared to Eqn. (5), then it is obvious that S is also a symplex as it is also the product of a symmetric and a skew-symmetric matrix and obeys Eq. 6. From Eqns. (13), (14) and (20) it follows that The second moments of a matched distribution are unchanged after one turn (or sector) of period L so that σ(L) = σ(0) so that one obtains in a few steps 4 :

B. The Eigensystems and Matching
Hence one finds that the matrices M, F and S have the same eigenvectors -but in general different eigenvalues 35,36 : ω i are the oscillation frequencies and ε i the emittances. If E is known, the second moments of the matched distribution can be computed by replacing the eigenfrequencies by the emittances. If a sympletic transformation R is known, that brings F (and hence S and M) to blockdiagonal form, then one can simply use the usual Courant Snyder theory for one-dimensional systems 39 . In this case an explicit computation of the eigenvectors is not required.
It was shown in Ref. 4 , that the ten coefficients of the force matrix F or the S-matrix can be identified with energy E and momentum P of a particle and with electric and magnetic field ( E and B, respectively) seen by a charged particle in external fields. The meaning of this identification is, that the corresponding coefficients of F or S transform under symplectic transformations in the exact same way as the fields and the momentum transform under the corresponding boosts and rotations.
It was also shown that the envelope equations of coupled linear optics are isomorphic to the Lorentz force equation. The Lorentz group was found to be a subset of the two-dimensional symplectic group. The so defined "fields" ( E and B) of the EMEQ should not be confused with the real fields of the beamline elements or accelerator components.
This isomorphism has been named electromechanical equivalence (EMEQ). The ten possible symplectic transformations are identified with spatial and phaserotations, Lorentz boosts and so-called "phase boosts". The transformation properties are analogous to those in Minkowski space-time.
This structural analogy is the basic idea behind the electromechanical equivalence (EMEQ). Naturally, γ 0 is associated with the time-like components of 4-vectors (i.e. energy), the spatial matrices γ = (γ 1 , γ 2 , γ 3 ) T are associated with the momentum, the matrices γ 0 γ with the electric field and γ 14 γ 0 γ with the magnetic field. The pseudoscalar has been named γ 14 = γ 0 γ 1 γ 2 γ 3 (instead of γ 5 , as convention in QED). The remaining six matrices are γ 10 , which is the time-component of the pseudovector, (γ 11 , γ 12 , γ 13 ) T = γ 14 γ are the spatial components of the pseudo-vector and γ 15 = 1 is the unit matrix. A complete list is given in App. A, further details in Ref. 4 and in textbooks on quantum electrodynamics.
The EMEQ is given by the following nomenclature: with the f k given by Eqn. (12). Using the EMEQ, the eigenvalues of F (Eqn. 26 and Eqn. 27) can be expressed by: Force matrices of stable systems have purely imaginary eigenvalues 37 , so that for stable systems one has K 2 > 0 and K 1 > 2 √ K 2 . Using the notation of the EMEQ a general symplex F is given explicitely by In the following we describe a geometrical approach of decoupling that is inspired by the observation, that in the decoupled force matrix, the scalar products E · B and P · B vanish 4 . In Hamiltonian form (see Eqn. 38 below), also the product P · E is zero and only the components E, P x , E z and B y remain. It is therefore instructive to analyze the symplectic transformation properties of these scalar products. The product E · B is known to be invariant under rotations and Lorentz boosts. Formally it is a pseudo-scalar in contrast to the scalar component representing the mass. Hence one might loosely speak of "mass components" and use the abbreviations: The "mass components" are invariant under spatial rotations. We may therefore proceed with phase rotations and boosts. We introduce the following auxiliary vectors: so that K 2 from Eqn. 29 can be written as It is easy to see that g, r and b transform under spatial rotations just like usual vectors. It is also quite obvious that the vector g equals the usual Lorentz force and the vector b equals the "Lorentz force" of a particle with magnetic charge, as the role of E and B is exchanged compared to g in the algebraic way that corresponds to a duality rotation through an angle of π 2 4 . One finds the following products: We introduce the following abbreviations for a better readability The phase rotation generated by γ 0 yields: The transformation of the mass components is listed in Tab. I. From the discussion of the normal form of the force matrix in Ref. 4 it follows, that decoupling to blockdiagonal form is done by a transformation that makes Geometrically this means, that B has to be aligned along the y-axis and the vectors P and E should be in the plane perpendicular to B. In a first step, the decoupling of a two-dimensional harmonic oscillator requires the (partial) orthogonalization of the (3-dimensional) "vectors" E, B and P : which can be interpreted as a geometrical "decoupling". The alignment of B along the y-axis in a second step is simple. A transformation to what we call "Hamiltonian" form requires additionally to make E x = P z = 0, which can again by done in two steps, orthogonalization and subsequent alignment of E and P . The general form of symplectic transformations has been described in some detail in Ref. 4 , here we give only a brief summary. A symplectic transformation matrix R b is generated by a basic symplex γ b with b ∈ [0 . . . 9] and controlled by a parameter ε: The effect of a basic symplex γ b depends on its "signature", which is positive for symmetric and negative for skew-symmetric γ b : where the bold printed 1 is the unity matrix. Note that transformations with γ 2 b = −1 (+1) are called rotations (boosts), respectively. Explicitely, γ 0 is the generator of a "phase rotation", γ b b ∈ [7,8,9] are "spatial rotations" with respect to the x, y and z-axis and γ b b ∈ [4,5,6] generate "Lorentz boosts" with respect to the x, y and z-axis. The "phase boosts" generated by γ b b ∈ [1, 2, 3] are combinations of phase rotations and Lorentz boosts. The parameter ε is called "angle" in case of rotations and "rapidity" in case of boosts. As the decoupling requires a sequence of transformations, we emphasize that the RDM-coefficients have to be updated according to Eq. 11 after each transformation.
Inspection of Tab. I shows that a straightforward strategy is the following: • M g → 0: Make a phase rotation generated by γ 0 with angle ε = arctan ( Mg Mr ). This will always work independent on the size of M i .
• b → | b| e y : Align the vector b along the y-axis by the spatial rotations with R 7 and an angle of ε = arctan ( bz by ) and with R 9 through an angle of ε = − arctan ( bx by ). Such rotations can always be done.
The last transformation is only possible, if |M r | < |b y | = | b|: 42) The first transformations result in P · B = 0, so that Eq. 42 is identical to the requirement that K 2 ≥ 0 (see Eq. 29). This means that the eigenvalues are either located on the real or imaginary axis, but not off-axis in the complex plane. If this condition is fulfilled, then the vector-components ( g) y and ( r) y , ( b) x and ( b) z are zero after the decoupling transformations have been applied. It follows from M r = E · B = 0 and M g = P · B = 0 and Eq. 32 that E · b = 0 and P · b = 0, and since we aligned b along the y-axis, we have E y = 0 and P y = 0, so that with b also B is aligned along the y-axis and B x = B z = 0. If we compare this with Eq. 30, then we note that the matrix F is now block-diagonal.
That is: we found a symplectic decoupling algorithm for both -systems with purely imaginary eigenvalues, which are called "strongly" stable 37 , and unfocused systems with purely real eigenvalues. That the algorithm works in both cases equally well, is important for instance in the case of transverse-longitudinal coupling with space charge in cyclotrons 3 .
We continue the discussion of force matrices with eigenvalues off axis in the complex plane in Sec. V B and assume for now, that K 2 > 0. Using the abbreviations the RDM-coefficients of the block-diagonal (decoupled) force matrix are given by: In order to bring the block-diagonal force matrix to Hamiltonian form, one may apply the following transformations: • M b → 0: Use another phase rotation with γ 0 with ε = 1 2 arctan ( 2 M b E 2 − P 2 ) • P z → 0: Use rotation about y-axis with γ 8 with ε = − arctan ( Pz Px ).
After these two rotations, the matrix has Hamiltonian form, if K 2 > 0 holds. In charged particle optics this is usually the case and therefore we consider this method as a generally applicable decoupling algorithm.

B. Complex Eigenvalues
Even though the problem of complex eigenvalues has not yet been solved for the general 2 n × 2 n case, it is possible to give a solution for the 4 × 4 case as we are going to describe here. The more general case of arbitrary 2 n × 2 n-symplices with arbitrary (complex) eigenvalues can presumably be solved by a block-diagonalization with 4 × 4-blocks for each set of complex conjugate eigenvalues and 2 × 2-blocks for each pair of real or imaginary eigenvalues.
If K 2 < 0 the eigenvalues are complex and a blockdiagonalization with 2 × 2-blocks is not possible (within the reals). However a simplification of the matrix is possible with the aim, that the RDM-coefficients of the transformed matrix have the following structure: so that one finds g = 0 and b = 0 and the auxiliary vector r has only a single non-vanishing component r x . We distinguish two cases, the first with E 2 < Max( P 2 , E 2 ) and the second with E 2 > Min( P 2 , E 2 ). In both cases the goal is to let "energy" and "momentum" vanish by appropriate Lorentz or phase boosts. Then one may align B along the y-axis and rotate about the y-axis to make E x = 0. Then the conditions of Eqn. 45 are fulfilled.

The Low Energy Case
The decoupling strategy for the first case, i.e. for E 2 < Max( P 2 , E 2 ): • M g → 0: Apply a phase rotation R 0 with angle ε 1 = arctan ( Mg Mr ). Note that this maximizes M r = E · B.
• E → | E| e y : Align the vector E along the y-axis by the spatial rotations with R 7 and an angle of ε 2 = arctan ( Ez Ey ) and (after an update of the RDMcoefficients and a recomputation of the auxiliary vector and mass components) with R 9 about an angle of ε 3 = − arctan ( Ex Ey ).
• E → 0: Boost using R 2 and rapidity ε 4 = arctanh( E Ey ). According to the assumptions, this is possible and does not change E x = 0 or E z = 0.
• P z → 0: Boost using R 1 and rapidity ε 6 = arctanh( Pz By ). Since E = E z = E x = 0, the energy E as well as E are unchanged by the boost.
• B → | B| e y : Align the vector B along the y-axis by the spatial rotations with R 7 and an angle of ε 7 = arctan ( Bz By ) and (after an update of the RDMcoefficients and a recomputation of the auxiliary vector and mass components) with R 9 about an angle of ε 8 = − arctan ( Bx By ).
• E x → 0: Rotate about the y-axis with R 8 with an angle of ε 9 = arctan ( Ex Ez ).

The Intermediate Energy Case
The case where E 2 > Min( P 2 , E 2 ) but K 2 < 0 might be called "intermediate", since the energy is large compared to the "low energy" case, but not large enough to make K 2 > 0. The following procedure leads to the state described by Eqn. 45: . Note that this transformation minimizes P 2 .
• P → | P | e y : Align the vector P along the y-axis by the spatial rotations with R 7 and an angle of ε 2 = arctan ( Pz Py ) and (after an update of the RDMcoefficients and a recomputation of the auxiliary vector and mass components) with R 9 about an angle of ε 3 = − arctan ( Px Py ). Since M b = E · P = 0 one also has now E y = 0.
• B → | B| e y : Align the vector B along the y-axis by the spatial rotations with R 7 and an angle of ε 6 = arctan ( Bz By ) and (after an update of the RDMcoefficients and a recomputation of the auxiliary vector and mass components) with R 9 about an angle of ε 7 = − arctan ( Bx By ).
• E x → 0: Rotate about the y-axis with R 8 with an angle of ε 9 = arctan ( Ex Ez ). In both cases the transformed matrix F then has the form In order to bring it to Hamiltonian form, one applies the transformation R 2 with an "angle" of i π/2: so that F → R F R −1 is: Note that the complex eigenvalues of a force matrix with K 2 < 0 all lie on a circle of radius ρ = (K 2 1 + 4 |K 2 |) 1/4 in the complex plane.

C. Decoupling n-dimensional Symplices
The general eigenvalue problem of symplices (Hamiltonian matrices) is an area of intense research. The algorithm presented above is based on a physical and geometrical analysis of 2-dimensional linear symplectic systems. As described before, the algorithm is limited to symplices that have real or imaginary eigenvalues, but a generalization to include complex eigenvalues might be possible -even though not urgently required in charged particle optics 5 .
In order to decouple symplectic systems with more than two degrees of freedom, the described algorithm can be used in an iterative scheme analogous to the Jacobi method for symmetric matrices 6 . If all eigenvalues are real or imaginary, it is possible to avoid computations using complex numbers. The 2 n ×2 n-symplex is then regarded as a n×n matrix of 2×2-blocks. We tested a pivot search that picks the maximum average square amplitude of all non-diagonal blocks B ij . The blocks B ii , B ij , B ji and B jj are then analyzed as 4×4-symplices and the symplectic similarity transformation that block-diagonalizes this submatrix is applied, so thatB ij =B ji = 0 holds. This iterative scheme allows to compute simultaneously the symplectic transformation matrix and the resulting block-diagonal or Hamiltonian form symplex with high precision.
Given {x} is a sequence of random numbers between zero and one, then one may construct random symmetric 2 n × 2 n-matrices A according to the rule: The increase of the diagonal terms helps to avoid complex eigenvalues. The symplex to decouple is then given by 5 The S-matrix for instance will never have complex eigenvalues as it is derived from the matrix of second moments. Complex eigenvalues would only be possible for correlations with a modulus greater than one. 6 Jacobi introduced a method to iteratively diagonalize real symmetric matrices by a sequence of orthogonal transformations each of which diagonalizes a 2 × 2 submatrix 43 . F = γ 0 A. We tested the algorithm with these random matrices up to n = 12 and logged the number of 4 × 4diagonalization steps. Fig. 1 shows the average number of iterations that is required to compute the transformation that brings a 2 n × 2 n symplex to Hamiltonian form, i.e. into the form:

D. Diagonalization
In order to proceed from Eqn. 38 towards diagonalization, the matrix is scaled using the generators γ 3 and γ 4 : so that one obtains (for stable systems), what we call "normal" form: where the signs of the frequencies ω 1 and ω 2 can both be positive and negative, depending on the signs of α, β, γ and δ. At this stage, all components of g and r as well as ( b) x and ( b) z are zero. Only ( b) y is non-zero.
That is -the last transformation matrix that is required for diagonalization is not only symplectic -it is also unitary.

E. Example
A simplified and idealized cyclotron model with space charge was described, which served as an example for an irregular system 3,4 . Without repeating all details, the constant force matrix has the following form: The RDM-coefficients are then given by: From this one finds for the "mass" terms and the vectors g, r and b: According to the geometrical approach, the first transformation can be omitted, since the "mass" M g is zero.
The second transformation using γ 7 aligns b along the y-axis. The second rotation may again be omitted, since the vector r is already aligned along the x-axis. The last transformation is a phase boost using γ 2 and is sufficient to bring F into block-diagonal form. This transformation would usually change the value of M b , but here it does not, since M b = ( r) y = 0 as can be seen from Tab. I. Hence M b remains zero -M g is invariant under both transformations. Hence, all "mass terms" are then zero after the described two transformations so that the system is decoupled.

F. Operators, Expectation Values and Lax Pairs
Coupled linear optics is in its essence (as quantum mechanics) a statistical theory. Since the reference trajectory is fixed, the coordinates are always taken relative to the local reference frame and the geometry is (only) locally euclidean. Even though the starting point is the description of single particle motion, the orbits of single particles are usually both, hard to access experimentally and of low practical value. The description of the beam by average values in contrast is both -measureable and of high value. The use of symplectic transformations leaves the expectation values unchanged. We can therefore evaluate the expectation values of any operator O in an arbitrary reference frame: since for symplectic R we have The time derivative of the expectation value of an arbitrary operator O, that does not explicitely depend on time, is: Equations of the form (here S = σ γ 0 ) appear frequently in the theory of coupled linear optics and it is worth mentioning that Eq. 60 is a so-called Lax representation and the operators S and F are a so-called Lax pair 40,41 . As a consequence, the expressions are first integrals of motion, where T r() is the trace. Using again the EMEQ to express the elements of S, one finds:

VI. SUMMARY AND OUTLOOK
A powerful method for symplectic decoupling of the ndimensional non-dissipative harmonic oscillator has been developed. The method apparently is stable, of the order O(n 2 ) and works with purely real or purely imaginary eigenvalues, for which a Hamiltonian Schur form does not always exists 20 . The resulting block-diagonal symplex can be used to compute the σ-matrix of matched beam ellipsoids of linear coupled systems in charged particle optics 3,4,36 . Another application is the production of multivariate gaussian distributions for a given covariance matrix 28 .
The presented parametrization gives deep insight into the general nature of coupling and might be instructive also in other areas of physics. The algebraic problem of finding the eigenvalues and eigenvectors of a twodimensional symplectic system was solved using geometrical arguments based on the use of the real Dirac matrices and the electromechanical equivalence.

ACKNOWLEDGMENTS
We would like to mention the work of D. Hestenes, who emphasised the geometrical significance of the Dirac algebra that he called space-time algebra 8 . The idea to introduce the EMEQ is inspired by his work.
Mathematica R has been used for some of the symbolic calculations. Additional software has been written in "C" and been compiled with the GNU c -C++ compiler 3.4.6 on Scientific Linux. The CERN library (PAW) was used to generate the figure.

Appendix B: Floquet Theorem
If the matrix A in Eqn. (1) and hence the forces are not constant, but periodic (F(t + T ) = F(t)), then Floquet's theorem can be applied and the solution has the general form 30,33 : where K(t) is symplectic and periodic with period T .
The transfer matrix of one period of length T ("one-turntransfer-matrix") is identical to the transfer matrix for a system with the constant force matrixF and length T . In this senseF is the "average" or "effective" force matrix with respect to one turn and can formally be written as 30 : From Eqs. B1 one derives in a few steps 30,38 : If the canonical transformation represented by K has been applied to the state vector, then with K(0) = 1 it follows: Note that the knowledge of K is not required to solve the decoupling problem, as long as the one-turn-transfer matrix M T is known. M T can either be obtained as a product of the transfer matrices of all beamline elements or simply by numerical integration. If the matched beam distribution has been found at an arbitrary (known) position s = 0 along the closed reference orbit, then the matched distribution can be computed for any position s using: Appendix C: Quick Guide to Decoupling To start with it is required to have either the average or constant force matrix F or the symplectic transfer matrix M that represents a complete turn or (cyclotron) sector. In the latter case one computes an auxiliary force matrix by while the usual (effective) force matrix has the form ω i being the betatron frequencies. The auxiliary matrix has the same structure but different eigenvalues s i = sin (ω i τ ), where ω i τ = 2π Q i with the betatron tunes Q i . Now compute the RDM-coefficients according to: Note that the coefficients for γ k with k ∈ [10, . . . , 15] must be zero -otherwise the system is not symplectic. Then compute the eigenvalues and auxiliary vectors r, g, b according to Eq. 29, 31 and 32. Construct the transformation matrices R b according to: Transform with γ 0 and ε = arctan Mg Mr : Recompute RDM-coefficients, then transform using γ 7 with ε = arctan bz by . Recompute RDM-coefficients, then transform using γ 9 with ε = arctan bx by . Recompute RDM-coefficients, then transform using γ 2 with ε = arctanh Mr by . The (auxiliary) force matrix should now be block-diagonal. Recompute RDM-coefficients, then transform with γ 0 and ε = arctan 2 M b E 2 − P 2 . Recompute RDM-coefficients, then transform with γ 8 and ε = − arctan Pz Px . Now the (auxiliary) force matrix should have normal form, so that the frequencies (or their sines) are given by: The complete transformation is given by: If the auxiliary matrix has been used, then compute the matrixM c according tõ The cosines of the tunes are then given by: cos (ω 1 τ ) + cos (ω 2 τ ) = T r(M)/2 cos (ω 1 τ ) − cos (ω 2 τ ) = T r(M γ 12 + γ 12M )/4 (C10)

Appendix D: The Teng and Edwards Ansatz
Assume that we have an even number of DOF, so that a 4 n × 4 n symplectic matrix R can be written in blockform according to 1,2 : where all quadratic submatrices are of size 2 n × 2 n, then the matrix R is symplectic, if where γ 0 has -in dependence of the context -to be taken as 2 n × 2 n or 4 n × 4 n.
If one now assumes that A = B = C 1, then it follows that If one assumes furthermore with Teng and Edwards, that C = cos (φ), then may define a = sin (φ) a s and b = sin (φ) b s with symplectic matrizes a s and b s , respectively: It has been shown in Ref. 3 , that C = cos (φ) is not the general case, since one might also choose C = cosh (φ), a = sinh (φ) a s and b = sinh (φ) b s . In this case one finds i.e. the matrizes a s and b s can also be antisymplectic (symplectic with multiplier −1). Still the matrix R remains symplectic. Hence Teng and Edwards limited their treatment in two ways: First, they assumed that C = cos (ψ) such that R must be a rotation matrix and secondly, they considered only the case that a and b are symplectic.