Calculation of wakefields in 2D rectangular structures

We consider the calculation of electromagnetic fields generated by an electron bunch passing through a vacuum chamber structure that, in general, consists of an entry pipe, followed by some kind of transition or cavity, and ending in an exit pipe. We limit our study to structures having rectangular cross-section, where the height can vary as function of longitudinal coordinate but the width and side walls remain fixed. For such structures, we derive a Fourier representation of the wake potentials through one-dimensional functions. A new numerical approach for calculating the wakes in such structures is proposed and implemented in the computer code ECHO(2D). The computation resource requirements for this approach are moderate and comparable to those for finding the wakes in 2D rotationally symmetric structures. Numerical examples obtained with the new numerical code are presented.


I. INTRODUCTION
The interaction of charged particle beams and the vacuum chamber environment can be quantified using the concept of impedance or wakefield [1]. In order to calculate electromagnetic fields in accelerators several different numerical approaches have been suggested and implemented in computer codes. Such calculations for complicated three dimensional structures, however, remain a challenge even for today's parallel computers. The approaches suggested in papers [2]- [4] allow one to obtain reliable results for rotationally symmetric and general three-dimensional structures even when using a personal computer. However, fully three dimensional calculations can be cumbersome and require long computation times. In contrast, for rotationally symmetric structures, the modelling is much simpler, the execution time is much shorter, and the required computational resources are relatively modest.
Many three-dimensional vacuum chamber components used in particle accelerators can be well-approximated with structures of rectangular cross-section whose height can vary as function of longitudinal coordinate but whose width and side walls remain fixed. Three examples of such geometry are shown in Fig. 1: (a) a corrugated structure that can be used as a dechirper [5], (b) a rectangular tapered collimator, and (c) a vacuum chamber in a bunch compressor. Calculation of the wakefields for such structures can be greatly simplified by expanding the field generated by the beam into Fourier series. As it turns out, each Fourier harmonic can be solved separately, and the wakefield calculated as a sum of their contributions. This reduces the original 3D problem into a number of 2D ones, each of which requires much less computing resources than the original problem.
In this paper we derive the equations for the Fourier harmonics and propose a computational scheme for their solution. The scheme was implemented in a computer code that we call ECHO(2D). The paper is organized as follows. In Section II, we apply the Fourier transform in the horizontal direction (the direction of constant width), and derive a system of equations for the Fourier harmonics. In Section III, we prove that the longitudinal wake function satisfies the Laplace equation with respect to the coordinates of the source particle. In Section IV, we show that, in the general case, the wakefield of each harmonic can be described by four scalar functions of one variable.
The description is further simplified in Section V where we assume that the entire structure has a horizontal symmetry plane, and then analyze the longitudinal and transverse wakefields near the system axis. An important symmetry relation for the wake function is established in Section VI. In Sections VII and VIII we present a numerical algorithm for solving equations derived in Section II. Several numerical examples obtained with the computer code ECHO(2D) are described in Section IX.
The results of the paper are summarized in Section X.

II. PROBLEM FORMULATION AND FOURIER EXPANSION
FIG. 2: Charged particle bunch moving through an accelerating structure.
We consider a bunch of particles moving with the velocity of light c through a metallic structure comprising an incoming pipe, a transition or cavity, and an outgoing pipe, an example of which is shown in Fig. 2. The bunch is characterized by the charge distribution ρ and the electric current density j = cρ. We assume that the bunch is moving along a straight line parallel to the longitudinal axis of the system, and we neglect the influence of the wakefields on the bunch's motion.
Our problem is to find the electric and magnetic fields E, H, in the domain Ω which is bounded transversally by perfectly conducting metallic walls, defined by the boundary ∂Ω of the domain Ω. We need to solve Maxwell's equations with the boundary conditions: where E 0 , H 0 , is an initial electromagnetic field in the domainΩ and n is a unit vector normal to the surface ∂Ω. While we consider the propagation of beams through vacuum, for the generality of the computational algorithm we included in Eqs.
We choose a coordinate system with y in the vertical and x in the horizontal directions; the z coordinate is directed along the longitudinal axis of the system.
The structures considered in this paper have constant width 2w and a fixed location of the side walls (in the x-direction).
The charge density is localized within the interval 0 < x < 2w and vanishes at the ends of the interval, ρ(0, y, z) = ρ(2w, y, z) = 0. Such a function can be expanded in Fourier series It follows from the linearity of Maxwell's equations, Eqs. (1) and Eq. (2), that the components of electromagnetic field can be represented by infinite sums: For each mode number m we can write an independent system of equations Hence we have reduced the 3D problem, Eqs. (1), to an infinite set of independent 2D problems, Eqs. (4). Our approach is similar to one used in [6] for solving Maxwell's equations in the frequency domain.

III. LONGITUDINAL WAKE AS A HARMONIC FUNCTION
Let us consider a line-charge beam with vanishing transverse dimensions, W (x 0 , y 0 , x, y, where the electric field on the right-hand side is the solution to Maxwell's equation with the sources of Eqs. (5) (this field, of course, is also a function of x 0 and y 0 omitted in the arguments of E z for brevity).
It is well known [7] that the longitudinal wake potential satisfies the Laplace equation for the coordinates x, y, of the witness particle In order to express the wake potential in structures of constant width through one-dimensional functions we need to prove that the longitudinal wake potential satisfies also the Laplace equation with respect to the coordinates x 0 , y 0 , of the source particle. This statement will actually be proven below for arbitrary structures without any restrictions on the structure geometry. The only requirement imposed on the system is that the incoming and outgoing waveguides have perfectly conducting walls.
The proof is based on the directional symmetry relation between the wake potential W for the case when the particle travels through the system in the positive z direction, and the wake potential W (−) corresponding to the case when the particle travels through the same system but in the negative z direction. It was shown in Ref. [8] that the two wake potentials are related by where 0 is the permittivity of vacuum, r 1 = (x 0 , y 0 ) and r 2 = (x, y) are two dimensional vectors, and φ A , φ B , are the Green functions for the Laplacian in the incoming and outgoing pipes with cross-sections, S A , S B : where ∂S is the boundary of S . Note that r 1 and r 2 are offsets of, respectively, the leading and the trailing particles in W , while r 2 is the offset of the leading particle, and r 1 is the offset of the trailing particle in W (−) . Note also that, according to the general property of the wake potentials (7), W (−) satisfies the Laplace equation where we have used the symmetry of Green's functions, φ A (r 1 , r 2 ) = φ A (r 2 , r 1 ), φ B (r 1 , r 2 ) = φ B (r 2 , r 1 ). We now apply the Laplacian operator to the coordinates x 0 , y 0 in Eq. (12) and use relation (11) to obtain This proves that the wake potential, in addition to being a harmonic function with respect to the coordinates of the trailing particle, is also a harmonic function with respect to the coordinates of the leading particle. To our knowledge, this general property of wakefields has not been previously reported in the literature.

IV. WAKE POTENTIAL REPRESENTATION IN STRUCTURES OF CONSTANT WIDTH
It is well known [7] that for rotationally symmetric structures the wake potential can be represented through one dimensional functions, with only one function for each azimuthal mode number m, As it turns out, a similar, comparably simple representation exists for the wake potential in structures of rectangular cross-section and constant width; however, in the latter case four one-dimensional functions are needed for each mode number m.
Substituting Eq. (5) into (2) we find for the charge distribution It then follows from Maxwell's equations (4) that all components of the fields will be proportional to sin(k x,m x 0 ) and using Eqs. (3) and (6) we find If we insert representation (16) into Eq. (7), we obtain the one-dimesional Helmholz which has a general solution of the form Finally from Eq. (13) it follows that functions W c m (y 0 , s) and W s m (y 0 , s) satisfy equations These equations can again be easily integrated with the result: where Thus, we have proven that, in structures of constant width, for each mode number m four functions are needed to completely describe the longitudinal wake potential.
These functions can be calculated as follows where the m th modal component of the wake potential is excited by a charge distribution that does not depend on x, With a knowledge of the longitudinal wake we can calculate the transverse wakes. For example, the vertical wake potential, W y , can be easily found through the Panofsky-Wenzel theorem [7] ∂ ∂s In the next section we will analyze both the longitudinal and the transverse wakes assuming that the structure under consideration has a symmetry axis.

V. STRUCTURES WITH HORIZONTAL SYMMETRY PLANE
Let us consider a structure of constant width 2w that also has a vertical symmetry plane, at y = 0. Structures in Fig. 1 (a) and (b) possess this symmetry; hence, they have a symmetry axis located at x = w, y = 0. Due to the symmetry, the wake potential satisfies the equation and Eq. (21) simplifies: Note that W m (y 0 , y, s) = W m (y, y 0 , s).
Let us consider the transverse wakes in such structures. We first introduce the integrated wake functions (sometimes called the step function response) It then follows from (25) that the transverse wake function can be written as where For small offsets near the symmetry axis, x = w, y = 0, the transverse wake potential is usually expanded in Taylor series, The first term in (32) is usually called the transverse dipole wake in the y-direction.
It can be calculated as follows The second term in (32) is called the transverse quadrupole wake in y-direction; it is obtained by The transverse wakes in the x direction are obtained by equations corresponding to those of Eqs. (33), (34). Note that W y,q (s) = −W x,q (s).
In numerical calculations of structures with symmetry we can use the approach of paper [9] that allows us to reduce the calculation domain in half. Indeed the charge distribution (24) can be written as a sum of symmetric and antisymmetric parts ρ m (y 0 , y, where ρ H m (y 0 , y, s) = In problems with the symmetric driving charges (36), the tangential component of the magnetic field will be zero in the symmetry plane (the so called "magnetic" boundary condition). In problems with the antisymmetric driving charges (37) the tangential component of the electric field will be zero in the symmetry plane (the "electric" boundary condition). Thus, instead of solving the system of equations (4) in the whole domain, one can solve two independent problems in half of the domain: one problem with the "magnetic" boundary condition at y = 0 and one problem with the "electric" boundary condition at y = 0. This is true not only for the line-charge current distribution (5), but for any arbitrary three dimensional charge distribution ρ(x, y, z, t). From solutions W H m (y 0 , y, s) and W E m (y 0 , y, s) of the two problems we can easily find the one dimensional modal functions in Eq. (27):

VI. SYMMETRY RELATIONS FOR THE TRANSVERSE WAKE POTENTIAL
In the general case, from the directional symmetry relation, Eqs.
Hence the directional symmetry relation (12) for the longitudinal wake potential can now be written as and from the Panofsky-Wenzel theorem (25) it follows that In the last equation we have used the vectorial notation for the transverse wake, Let us now calculate the directional relation for the dipole and quadrupole wakes The Green function for Poisson's equation in a rectangular pipe, φ A (x 0 , y 0 , x, y), has been obtained by Gluckstern, et al [10]. For x 0 = w, their result reads: 1 The last term is singular in the limit x → w and y → y 0 . All other terms are finite for all admissible x and y, and the sum in the first term rapidly converges.
The Green's function of the outgoing pipe, φ B (w, y 0 , x, y), is obtained from (42) by replacing g A → g B . It is then easy to see that the singular terms cancel in the

VII. TE/TM SCHEME IN MATRIX NOTATION
In this section we describe a particular realization of the implicit TE/TM scheme introduced in [4]. The scheme will be discussed in the context of the Finite Integration Technique [12].
where ϑ is a mesh multi-index and L ϑ , S ϑ ∈ G. Solving Faraday's law in integral form for the front surface shown in Fig. 3 yields: −e x,i, j,k + e y,i, j,k + e x,i, j+1,k − e y,i+1, j,k = − ∂ ∂t b z,i, j,k .
Note, that this representation is still exact, as e ϑ is (by definition) the exact electric voltage along one edge of the cell, and similarly b ϑ represents the exact value of the magnetic flux density integral over the cell surface. If we compose column vectors from all voltage and flux components, we can write the combination of all equations over all surfaces in an elegant matrix form as The matrix C picks the affected components out of the long vector to make up the corresponding equation. C is thus the discrete curl operator over the mesh G. With an appropriate indexing scheme the curl matrix has a 3x3 block structure: The double-banded, topological P x , P y , P z -matrices take the role of discrete partial differential operators.
The second important differential operator in Maxwell's equations (1) The discretization of the remaining Maxwell equations requires the introduction of a second cell complexG which is dual to the primary cell complex G. For the Cartesian grid the dual complexG is defined by taking the foci of the cells of G as gridpoints for the mesh cells ofG. We again introduce the computational unknowns as integrals where ϑ is a mesh multi-index andL ϑ ,S ϑ ∈G.
Following an equivalent procedure for the remaining Maxwell's equations, but using the dual cell complexG, we obtain a set of four discrete equations represent-ing Maxwell's equations on a dual grid: where the asterisk denotes the Hermetian adjoint operator.
Equations (51) are completed by the discrete form of the material relations (constitutive equations) which appear (in the simplest linear case) as matrix equations with M −1 the discrete inverse permittivity matrix, and M µ −1 the inverse permeability matrix. In the case of Cartesian grids (or, more generally, whenever the primary and dual grids are orthogonal) all material operators can be represented by diagonal matrices.
Note that the material matrices contain both averaged material parameters, and the lengths and areas of the grid edges and faces, respectively. In the standard staircase approximation the material matrices contain elements (without double indices for simplicity of notation) with p = x, y, z, and the face areas and edge lengths of the primary and secondary grid given by S , L,S ,L, respectively.
In the conformal scheme [14], System (51) is a time-continuous and space-discrete approximation of problem (4). The next step is a discretization of the system in time. The field components can be split in time and the leap-frog method can be applied. With "electric/magnetic" splitting, a well known Yee's scheme [13] will be obtained. However Yee's scheme has large dispersion errors along the grid lines. In the following we consider an alternative TE/TM splitting scheme [4], one that does not have dispersion errors in the longitudinal direction.
Let us rewrite Eqs. (51) in an equivalent form where τ = ct and Applying TE/TM splitting [4] of the field components in time to system (55), the following numerical scheme is obtained Two-layer operator-difference scheme (57) acquires the canonical form [15] B y n+1 − y n ∆τ + Ay n = f n , where It was shown in [4] that stability of the method is insured if If the matrix Q is positive definite, then we can define a discrete energy as W n = 0.5 Qy n , y n (60) and the discrete energy conservation law holds W n+1 − W n ∆τ = − ê n+0.5 ,ĵ n+0.5 , e n+0.5 = 0.5 x + e n x e n+1 y + e n y e n+0.5 z + e n−0.5 The stability condition (59) in the staircase approximation (of the boundary in vacuum) can be rewritten in the form This last condition resembles the well-known stability condition of the explicit FDTD scheme for one-dimensional problems. The maximal eigenvalue of the discrete operator P z P * z fulfills the relation λ max < 4, and the stability condition of the scheme in the staircase approximation of the boundary in vacuum reads On an equidistant mesh the implicit scheme (57) has a second order local approximation error in homogeneous regions.
Relation (63) does not contain information about the transverse mesh. Hence the transverse mesh can be chosen independent of stability considerations. Following the conventional procedure [4], the dispersion relation can be obtained in the form sin Ω ∆τ where Ω = 0.5ω∆τ, K y = 0.5k y ∆y, K z = 0.5k z ∆x. With the "magic" time step ∆τ = ∆z, the scheme does not have dispersion in the longitudinal direction.
Following the approach of [4] it is easy to show that charge conservation holds: whereh n x ,h n y are auxiliary vectors and operator matrix W e is a sparse one, The field components e n+1 x , e n+1 y , h n+1 z at time level n + 1 can be found as e n+0.5 x , e n+0.5 whereẽ n+0.5 x ,ẽ n+0.5 y are auxiliary vectors and operator matrix W h is a sparse one, In the staircase approximation of the boundary, the material matrices M µ −1 , M −1 are diagonal, and the operators (65),(67) are tri-diagonal matrices. The equations involving these operators can be solved easily by the "elimination" method [15].
However, the staircase approximation of the boundary results in first order convergence in the wake potential. In order to obtain second order convergence and avoid time step reduction, we use conformal method in the same way as described for the rotationally symmetric case in [16].
In the original problem formulation (1), we assumed that the structure walls are perfectly conducting. However, in order to use the Fourier expansion (3), it is sufficient to require only that the walls at x = 0 and x = 2w be perfectly conducting.
The bottom and top walls of the structure can have different boundary conditions; for example, they can be specified to be metallic with finite resistivity. A similar approach for the case of rotationally symmetric geometry and the use of the staircase approximation for the boundaries was described in [17]. In ECHO(2D) we have realized a conformal version of the method for including boundaries with finite conductivity. Details of the implementation of our algorithm will be published elsewhere.

IX. NUMERICAL EXAMPLES
In this section we consider two example calculations: a corrugated pipe ("dechirper") with perfectly conducting walls, and a tapered collimator with resistive walls. We compare results obtained with ECHO(2D) with those of several 3D codes.
For our dechirper example we used the parameters of the dechirper experiment performed at the Pohang Accelerator Laboratory [18]. The geometry is shown in In Table I   As the second example problem we consider a symmetric, tapered collimator  Fig. 9 we compare the longitudinal wake for this collimator with one that has the same geometry but is perfectly conducting. The Gaussian bunch in the simulations has an rms length σ z = 0.25 cm. Both wakes were obtained with ECHO(2D). In the right plot of Fig. 9 we compare the ECHO(2D) wake potential with the one obtained using a fully three-dimensional, commercially available code CST [20]. The good agreement between the results indicates a good accuracy of the conformal meshing and the resistive wall modelling in ECHO(2D).
Finally, in order to demonstrate the capabilities of the new code to calculate wakes of very short bunches, we give in Fig. 10  plotting symbols). Note that the shortest bunch length used in the calculations is σ z = 10 µm. For the shorter bunches, the loss factor in a 1 m-long structure is larger than the steady-state (periodic) result. We can see that by plotting on the same figure the difference of κ for a 2 m structure minus that for a 1 m structure (the gray curve). Note that the analytical, asymptotic value of the loss factor, for σ z → 0, is κ = Z 0 cL/(2πa 2 ) · (π 2 /16) = 1234. V/pC [we have taken Z 0 = 377 Ω, structure length L = 1 m, and half-gap a = 3 mm], which agrees well with the linear extrapolation of the ECHO(2D) steady-state results (the gray curve).

X. SUMMARY
In this paper we presented a new method for solving electromagnetic problems and calculating wakefields excited by a relativistic bunch in a structure that is characterized by a rectangular cross-section whose height can vary as function of longitudinal coordinate but whose width and side walls remain fixed. Using the Fourier expansion of the fields, currents and charge densities, we derived a Fourier repre-sentation of the wake potential in terms of four one-dimensional functions for each harmonic. We proved that the longitudinal wake is a harmonic function with respect of the coordinates of the leading charge. For a structure that has a horizontal symmetry plane we also proved an important symmetry relation of the wake function.
A numerical method was proposed for solving Fourier harmonics of the fields.
The method does not generate dispersion in the longitudinal direction and it conserves the energy and charge in the calculations. The computer resources required to solve the system with several tens of the lowest harmonics are moderate and comparable to those needed for 2D rotationally symmetric calculations.