Te/tm Field Solver for Particle Beam Simulations without Numerical Cherenkov Radiation

The Yee finite-difference time domain method (FDTD) is commonly used in wake field and particle-in-cell simulations. However, in accelerator modeling the high energy particles can travel in vacuum faster than their own radiation. This effect is commonly referred to as numerical Cherenkov radiation and is a consequence of numerical grid dispersion. Several numerical approaches are proposed to reduce the dispersion for all angles and for a given frequency range, that justifies itself for domains big in all three directions. On the contrary, in accelerator modeling the transverse dimensions and transverse beam velocity are small, but it is extremely important to eliminate the dispersion error in the well-defined direction of the beam motion for all frequencies. In this paper we propose a new two-level economical conservative scheme for electromagnetic field calculations in three dimensions. The scheme does not have dispersion in the longitudinal direction and is staircase-free (second order convergent). Unlike the FDTD method, it is based on a ''transversal-electric/transversal-magnetic'' (TE/TM)-like splitting of the field components in time. The scheme assures energy and charge conservation. Additionally, the usage of damping terms allows suppressing high frequency noise generated due to the transverse dispersion and the current fluctuations. The dispersion relation of the damping scheme is analyzed. As numerical examples show, the new scheme is much more accurate on the long-time scale than the conventional FDTD approach.


I. INTRODUCTION
The particle-in-cell (PIC) method [1] is an effective approach for simulation of beam dynamics in accelerators.The electromagnetic fields in many PIC codes are computed using the finite-difference time domain (FDTD) method [2,3].As any numerical mesh approach the FDTD method suffers from an anisotropic numerical dispersion.The numerical wave phase speed is slower than the physical one.Hence, the high energy particles can travel in vacuum faster than their own radiation.This effect is commonly referred to as numerical Cherenkov radiation [4], which (due to its accumulative character) corrupts the simulation.Hence, the electromagnetic field computation for short relativistic bunches in long structures remains a challenging problem even with the fastest computers available.
Several approaches [4 -8] have been proposed to reduce the accumulated dispersion error of large-scale threedimensional simulations for all angles and for a given frequency range.These methods require the usage of larger spatial stencils and a special treatment of the material interfaces.The increased computational burden justifies itself for computational domains large in all three dimensions.However, in the accelerator applications the domain of interest is very long in the longitudinal direction and relatively narrow in the transverse plane.Additionally, the electromagnetic field changes very fast in the direction of bunch motion but is relatively smooth in the transverse plane.Hence, it is extremely important to eliminate the dispersion error in the longitudinal direction for all frequencies.
As it is well known, the FDTD method at the Courant limit is dispersion free along the grid diagonals and this property can be used effectively in numerical simulations [9].However, the only reasonable choice in this case is to take equal mesh steps in all three directions.Alternatively, a semi-implicit numerical scheme without dispersion in the longitudinal direction with a simpler conformal treatment of material interfaces and the usage of nonequidistant grids has been developed in [10 -13].
The scheme described in [11] allows one to solve the scalar problem and to calculate the wake potential for fully axially symmetric problems with staircase approximation of the boundary.In [12,13], a three-level scheme R y n1 ÿ 2y n y nÿ1 Ay n f n for the vectorial problem was suggested.The scheme is based on a second order hyperbolic wave equation for vector potential.A modification of the uniformly stable conformal method [14] is used to avoid the ''staircase'' problem and to obtain a second order convergent algorithm.However, the operator R in the scheme is not selfconjugate; and therefore an ''energy'' conservation cannot be proven theoretically by the standard techniques [15].
Additionally, the scheme is not economical for general three-dimensional geometries.The last drawback can be overcome by splitting methods [16].However, the absence of a theoretical proof for an energy conservation has stimulated us to look for an alternative approach in the threedimensional case.
In this paper, a new two-level economical conservative scheme for electromagnetic field calculations in three dimensions is presented.The scheme does not have dispersion in the longitudinal direction and is staircase-free (second order convergent).Unlike the FDTD method [2] and the scheme developed in [12,13], the new method is based on a ''transversal-electric/transversal-magnetic'' (TE/TM)-like splitting of the field components in time.Additionally, it uses an enhanced alternating direction splitting of the transverse space operator that renders the scheme computationally as effective as the conventional FDTD method.Unlike the FDTD alternative direction implicit (ADI) [17] and low-order Strang [18] methods, the splitting error in our scheme is only of fourth order.
The new scheme assures energy and charge conservation.Additionally, the usage of damping terms allows one to suppress a high frequency noise generated due to the transverse dispersion and current fluctuations.The dispersion relation of the damping scheme is derived and analyzed.
Numerical examples show that the new scheme is much more accurate in long-time simulations than the conventional FDTD approach.For axially symmetric geometries, the new scheme performs at least 2 times faster than the scheme suggested in [12,13] while achieving the same level of accuracy.

II. FORMULATION OF THE PROBLEM
At high energies the particle beam is rigid.To obtain the electromagnetic wake field, the Maxwell equations can be solved with a rigid particle distribution.The influence of the wake field on the particle distribution is neglected here; thus, the beam-surroundings system is not solved selfconsistently and a mixed Cauchy problem for the situation shown in Fig. 1 should be considered.
The problem reads: for a bunch moving with the velocity of light c and characterized by a charge distribution , find the electromagnetic field Ẽ; H in a domain which is bounded transversally by a perfect conductor @.The bunch introduces an electric current j c and thus we have to solve for r H @ @t D j; r Ẽ ÿ @ @t B; r D ; r B 0; where Ẽ0 ; H0 is an initial electromagnetic field in domain and ñ means the normal to the surface @.In accelerator applications, the studied structure is usually supplied by an ingoing pipe, and the analytical solution in a perfectly conducting cylindrical pipe [19] can be used as initial condition.If the ingoing pipe is not cylindrical the initial field can be found numerically.

A. Finite integration technique
Following the matrix notation of the finite integration technique [20], the Cauchy problem (1) can be approximated by the time-continuous matrix equations on a grid doublet (G; G) completed by the discrete form of the material relations (constitutive equations) with the discrete inverse permittivity matrix M " ÿ1 and the inverse permeability matrix M ÿ1 .In the following the material matrices are assumed to be real and symmetric.On Cartesian fx; y; zg coordinate grids (like the Cartesian grid shown in Fig. 2) with an appropriate indexing scheme the curl and divergence matrices have the block structures: The two-banded topological P fu;v;wg matrices play the role of discrete partial differential operators [21].
With changing of variables e M ÿ1=2 " ÿ1 ê, h M ÿ1=2 ÿ1 ĥ, j c ÿ1 M 1=2 " ÿ1 ĵ, q q, ct, system (2) and ( 3) reduces to the skew-symmetric one (4) with the new matrices  4) is a time-continuous and space-discrete approximation of problem (1).For representation (4) a discrete analogue to the analytical equation div curl 0 holds: It allows one to show the charge and energy conservation in the scheme: W 0:5e T e h T h: Note that this energy has a direct correspondence to the total physical energy of the continuous electromagnetic fields given by 0:5 R V "jEj 2 ÿ1 jBj 2 dv [22].The next step is a discretization of the system in time.The field components can be split in time and the ''leapfrog'' scheme can be applied.Below, two kinds of the splitting are considered: E=M and TE/TM schemes.

B. Explicit FDTD method based on ''electric-magnetic'' splitting of the field components in time
Suggested by Yee [2], the ''electric-magnetic'' (E=M) splitting of the field components yields the explicit FDTD scheme (E=M scheme) e n0:5 e nÿ0:5 4C 0 h n ÿ 4j n ; h n1 h n ÿ 4C 0 e n0:5 ; ( where 4 is the time step, and the update of the electric components is shifted by 0:54 relative to the update of the magnetic components.Scheme ( 6) is a two-layer scheme B y n1 ÿ y n Ay n f n ; (7) where y n e nÿ0:5 We study the stability of scheme ( 7) by the energy inequalities method [15].Let us take the inner product of both sides in Eq. ( 7) with y n1 y n : hBy n1 ÿ y n ; y n1 y n i hAy n ; y n1 y n i hf n ; y n1 y n i: Using the formula y n 0:5y n1 y n ÿ y n1 ÿ y n we rewrite relation (8) in the form hB ÿ 0:5Ay n1 ÿ y n ; y n1 y n i 0:5hAy n1 y n ; y n1 y n i hf n ; y n1 y n i: The second term in the left-hand side is equal to zero since the operator A is skew symmetric and, therefore, hQy n1 ; y n1 i ÿ hQy n ; y n i hf n ; y n1 y n i; (9) where the self-adjointness of the operator Q B ÿ 0:54A is used.
The last relation allows one to prove that the condition is necessary and sufficient for the stability of the scheme.Following [22,23], a discrete energy of electromagnetic fields can be defined as W n E=M 0:5hB ÿ 0:54Ay n ; y n i 0:5he nÿ0:5 ; e nÿ0:5 i hh n ; h nÿ1 i; (10) and relation (9)  The charge conservation holds in the form: q n0:5 ÿ q nÿ0:5 4 S e j n 0; S h h n1 ÿ h n 4 0: The condition (9a) can be rewritten as where f i g are eigenvalues of the matrix C 0 C T 0 .It was proven in [23] that relation (9a) holds if 4y ÿ2 4z ÿ2 s ; i; j; k: Scheme ( 6) is widely used in electromagnetic modeling.However, the FDTD algorithm causes nonphysical dispersion of the simulated waves in a free-space computational lattice.The phase velocity of discrete wave modes can differ from the light velocity by an amount varying with the wavelength, direction of propagation in the grid, and grid discretization.With an equidistant mesh, a homogenous material and the time step equal to the right-hand side of inequality (11), the scheme has zero dispersion along the grid diagonals.Hence, the zero dispersion in a desired direction can be achieved by rotation of the mesh.However, this approach awakes limitations on discretization: the only reasonable choice in this case is to take equal mesh steps in the all three directions.The next difficulty arises with the attempt to use a conformal method.
Why is zero dispersion for a special direction important?Unlike plasma problems, the charged particles in accelerators are organized and a direction of motion (the longitudinal direction) can be identified.Hence, the computational domain is very long in the longitudinal direction and relatively short in the transverse plane.Additionally, the electromagnetic field changes very fast in the direction of motion but is relatively smooth in the transverse plane.
Note also that to be able to model smooth transitions in geometry we should use a conformal approach without time step reduction [14].

C. Implicit FDTD method based on transversal-electric/transversal-magnetic splitting of the field components in time
The arguments, stated in the preceding section, force us to look for a numerical scheme, which (i) does not have dispersion in the longitudinal direction, (ii) allows the use of nonhomogeneous meshes in the transverse plane, (iii) allows the use of a moving mesh without interpolation, (iv) allows accurate geometry modeling without a time step reduction.
In [12,13], a three-level implicit conformal scheme R y n1 ÿ 2y n y nÿ1 Ay n f n was suggested.The scheme is based on a vector potential formulation and allows an economical realization for axially symmetric geometries.However, the absence of a theoretical proof for an energy conservation has stimulated us to look for an alternative approach in the threedimensional case.
To find an alternative scheme, let us consider Fig. 2 and subdue an update procedure to the motion of the bunch.We assume that a charged particle is moving in the z direction with the velocity of light.Additionally, let us assume that our numerical scheme allows to take a time step 4 equal to the mesh step 4z in the z direction.If at the time 0 the particle has the position aligned with the left z facet of the primary grid (see Fig. 2), then at time 0 0:54 it will be aligned with the left z facet of the dual grid and in the time 0 4 it will be again aligned with the next z facet of the primary grid.This suggests that we should replace the E=M time splitting of the field components in scheme ( 6) by a more adequate TE/TM splitting.Indeed, at time 0 it is reasonable to update the ''TE'' components e x ; e y ; h z and half a time step later, namely, at time 0 0:54, we have to update the ''TM'' components h x ; h y ; e z .
Following these considerations, let us rewrite scheme (4) in the equivalent form where Applying the suggested TE/TM splitting of the field in time to system (12), the following numerical scheme is obtained: Just like scheme (6), scheme ( 13) is also a two-layer one where Analyzing relations ( 14) we conclude that just as for Yee's scheme the following relations hold: Likewise we can prove that condition (9a) is necessary and sufficient for the stability of scheme (14).
As for the E=M scheme the discrete energy in the TE/ TM scheme can be defined by the relation

W n
TE=TM 0:5hB ÿ 0:54Ay n ; y n i 0:5hu nÿ0:5 ; u nÿ0:5 i hv n ÿ 4D 21 u nÿ0:5 ; v n i W n E=M O4 2 : The relation ( 9) can be interpreted as energy conservation law Note that the energy W n TE=TM , just like the energy W n E=M defined by relation (10), is a second order accurate approximation to the total physical energy of the continuous electromagnetic field.If the right-hand side in scheme (14) vanishes, the scheme is energy conserving: W n TE=TM W 0 TE=TM : The charge conservation holds in the form: The last relation resembles the well-known stability condition of the explicit FDTD scheme for the onedimensional problem.In the following an equal mesh step 4z in the z direction will always be assumed.Then for a vacuum domain with staircase approximation of the boundary the stability condition reads With the time step 4 equal to the longitudinal mesh step 4z, scheme (13) does not have dispersion in the longitudinal direction (see the dispersion relation in Sec.III E).Relation (15) does not contain information about the transverse mesh.Hence, the transverse mesh can be chosen independently from stability considerations.For a relativistic bunch a mesh moving together with the bunch can be used.The field ahead of the bunch is zero and, as the scaled time step is equal to the longitudinal mesh step, the complete information for updating of the last mesh layer is available, too.This means that interpolation procedures are avoided and the dispersion in the longitudinal direction is equal to zero.The results with the moving mesh for staircase approximation of the geometry are fully equivalent to the stationary global mesh approach.
So far we have found a scheme which with staircase geometry approximation fulfills the first three requirements formulated above.However, in a general case the staircase scheme is only first order accurate.In order to overcome this problem and avoid reduction of the stable time step, the uniformly stable conformal approach described in [13,14] will be used.
With the latter approach the scheme possesses the desired features.However, it is implicit and noneconomical.An economical scheme modification, based on operator splitting, will be considered in Sec.IV

D. Noise control in the TE/TM scheme with implicit enforcement of charge conservation
The derived TE/TM scheme is energy conserving.It does not have dispersion in the longitudinal direction.However, the dispersion in other directions and current fluctuations can result in high frequency noise.To overcome the problem we use the idea of transverse current adjustment [24,25].Our modification of this approach allows us to suppress the noise without introducing dispersion in the direction of bunch motion.
In order to obtain a damping the TE/TM scheme is changed to the form u n0:5 ÿ u nÿ0: 5  4 D 11 u n0:5 u nÿ0:5 2 d u n0:5 ÿ u nÿ0:5 D 12 v n j nd u ; The above equations are equivalent to changing the time centering of the transverse part of the curl operator.These modifications have the additional benefit that they are trivial to implement in the implicit scheme (only some coefficients in the scheme are changed).
In scheme (17) the matrix Q B ÿ 0:54A is not selfconjugate and the energy dissipates.However, the charge conservation holds in the form: q n0:5 e ÿ q nÿ0:5 e S e j n 0; q n1 h ÿ q n h 0; q n0:5 It can be proven from consideration of the equation The modified TE/TM scheme allows for controllable by parameter d damping of the high frequency noise.

E. Dispersion in the damping TE/TM scheme
Following the conventional procedure [26], the dispersion relation for the damping scheme can be obtained in the form where 0:5!, K x 0:5k x x, K y 0:5k y y , K z 0:5k z z.
With the magic time step z, the scheme does not have dispersion in the longitudinal direction.Indeed, the scheme is able to propagate the full field pattern of the relativistic bunch very accurately in the vacuum perfectly conducting pipe.It was tested for the shortest rectangular bunch with the length equal to the mesh step z.Yee's E=M scheme shows in this test the high frequency noise and fast degradation of the pattern due to the longitudinal dispersion.
In order to estimate the dispersion in the transversal plane we consider a plane wave with the wave vector k x ; k y ; k z 1; 0; 0. For the mesh size x the dispersion relation reads ! 2 arccos cosec0:5k x ÿ 2id 1 ÿ 4d 2 ÿ 4id cosec0:5k x cosec 2 0:5k x p : Figure 3 shows the dispersion for different values of the damping parameter d.As can be seen from the figure, the case d 0:125 allows an effective damping of high frequency waves without deterioration of the dispersion curve.
The transversal velocity components of the charges in accelerators are nonrelativistic and the transverse dispersion does not cause numerical Cherenkov radiation.

IV. AN ECONOMICAL TE/TM SCHEME BASED ON TRANSVERSE OPERATOR SPLITTING
In order to find an economical scheme, three different schemes based on an operator splitting were considered in [27] and it was shown that the ADI2 approach results in an accurate numerical method with moderate restriction on the time step.
The numerical scheme using the ADI2 splitting in three dimensions has the form As for the TE/TM scheme (13) the relations hold.However, the stable time step does now not only depend on the longitudinal mesh step 4z but also on the minimal mesh step in the transverse plane and the stability condition reads [27] 4 min24x i ; 24y j ; 4z: The last condition does not reduce the applicability of the scheme, since the field is relatively smooth in the transverse plane and a much coarser grid can be used here.Constraint (19) can be relaxed and replaced by (16) if the ATIp method is used [27].So far we have not defined the ADI2 terms in relation (18).Instead of doing this, let us rewrite scheme (18)  If the material matrices M ÿ1 ; M " ÿ1 are diagonal, then systems(20c) and (21c) only have products of tridiagonal matrices on the left-hand side and can be solved easily.For example, Eq. (20c) leads to the set of equations where the vector F n denotes the right-hand side of Eq. (20c) and M " ÿ1 z P y M ÿ1 x P y ; M " ÿ1 z P x M ÿ1 y P x : However, the conformal scheme with the diagonal material matrices reduces the stable time step.To restore stability condition (19) and the possibility to use the time step 4 4z, we will use a modification of the uniformly stable conformal method [14] as described in detail in [13].The last approach results in modified nondiagonal but symmetric matrices M ÿ1 x ; M ÿ1 y .Other material matrices in scheme (20) and (21) remain diagonal.This means that we do not encounter difficulties in the solution of Eq. (21c).However, solving Eq. (20c) requires additional efforts since the matrices (23)are not tridiagonal.
To overcome the problem we modify system (22) as follows: where and M 0 I 0:5A 0 ên0:5;i z I 0:5A 0 ênÿ0:5 z u i 2 ; i 0; 1; 2; . . .; p: (25) Note that the equation for zero iteration [which we refer to as the TE/TM-ADI2(0)], just as schemes ( 22) and ( 24), results in an approximation of the continuous problem (1) with the error Oj4rj 2 4 2 .However, the TE/TM-ADI2(0) scheme may be unstable for general geometries and the required time step 4 4z.The first iteration [which we refer to as the TE/TM-ADI2( 1) scheme] solves the stability problem for all considered cases.
In the next section we will study properties of the scheme ( 18)-( 25) for the case of rotationally symmetric geometries.In the last section results for the fully three-dimensional scheme will be presented.

A. Realization of TE/TM and TE/TM-ADI2 schemes for rotationally symmetric geometries
In this section we describe the realization of the TE/TM scheme for the case of rotationally symmetric geometries.We consider this case separately since the TE/TM scheme ( 13) is already economical and the application of the splitting method considered in the previous section can be avoided.
For a bunch moving at the speed of light c at an offset a from and parallel to the axis of a rotationally symmetric structure, the source current j can be represented as where s is the longitudinal charge distribution and m is the azimuthal mode number.Numerical scheme (13) for an azimuthal mode number m has the form ĥn0:5 where and the fact that P ' mI is used.
If the material matrices M ÿ1 ; M " ÿ1 are diagonal, then operators (27) are tridiagonal matrices and equations involving them can be solved easily.For the case of nondiagonal matrices M ÿ1 ' ; M ÿ1 r , we will proceed in the same way as described at the end of the previous section.
We rewrite the equation with the operator W e CN in the form where and M 0 matrices.System (28) can be solved iteratively i 0; 1; 2; . . .; p: ( Scheme ( 26)-( 29) will be referred to as TE=TMp.Note that just as for the TE=TM-ADI2p scheme, it is sufficient to perform only one iteration [scheme TE/TM (1)] to obtain a stable solution.
As noted at the beginning of this section, for geometries of revolution we do not need to apply the transverse operator splitting.However, in order to check the achieved accuracy, the TE=TM ÿ ADI2p scheme was implemented for rotationally symmetric geometries, too.
The initial field is converted to the cylindrical coordinates and set in the entire calculation domain.After a period of time T 2 p a we compare the numerical solution with the exact one.A series of equidistant meshes with the cell sizes 4r 4z h is used.
Figure 4 shows the results for the time step 4 4z and for the mesh resolution a=h 10.The left-hand figure shows convergence of the noniterative schemes TE/TM(0) ( 26)-( 29) and TE/TM-ADI2(0) ( 18)- (25).Both schemes achieve the same rate of convergence.The right-hand figure shows conservation of the discrete energy E n TE=TM for the schemes.With an increase of the number of iterations p the discrete energy in the TE=TMp scheme converges to the constant value of the noniterative TE/ TM scheme.However, in order to see the same effect for the TE=TM-ADI2p scheme, the energy norm has to be changed to the one with operator Q from the noniterative scheme (18).We do not show this result here, since it already follows from Fig. 4 that the considered schemes are stable for the time step 4 4z, when they do not have dispersion in the longitudinal direction.As a further test example we use the circular collimator structure shown on the left-hand side of Fig. 5 (with inner radius b not indicated in the figure).Figure 6 shows the results for the dipole wake field (m 1) and compares the TE/TM scheme results to the ones obtained with the classical Yee's scheme (E=M scheme).The latter results are calculated with the help of code ABCI [29] (FDTD method with triangular approximation of the boundary).The geometric parameters are a 35 mm, L 20 cm, and b c 2 mm, where b is the inner radius of the collimator.The left-hand drawing in Fig. 6 shows the transversal dipole wake potential [30] W 1 ? s jW 1 ?s; r; 0jr ÿ1 ; W 1 ?s; r; for the collimator with L 20 cm and a relativistic Gaussian bunch with rms length 1 mm.The solid curves show the results for ABCI and the dashed ones present the results for the new scheme.
In the right-hand figure the transversal dipole loss factor for the collimator is shown for different mesh resolutions =h, where h 4z 4r is the mesh step.The error compared to the reference value (obtained with the finest mesh resolution) is also shown in the figure.The dashed line shows the results for the TE/TM code and the solid line for ABCI.
From the above example we see that the reference code ABCI demands a much more dense mesh for the same accuracy.
Finally, we show in Fig. 7 (left) the dipole wake potentials of a Gaussian bunch with 1 mm for the TESLA cryomodule of total length 11 m [31].The cryomodule contains 8 cavities and 9 bellows whose geometries are outlined in Fig. 8.The iris' radius is 35 mm and the beam tube's radius is 39 mm.
The moving mesh in the last example covers the bunch longitudinally in the range from ÿ5 to 100.The length of the moving mesh is only 0.105 m, which results in a drastic reduction of the computational demands (storage and CPU time) compared to the stationary mesh of total length 11 m.
The right-hand side of Fig. 7 shows the difference between the results obtained by the TE/TM(1) and TE/ TM-ADI2(1) schemes and the reference result calculated with the vector potential method (POT-2.5)described in [12,13].The presented results are calculated with the mesh resolution 4z 4r =5.It can be seen that the new methods introduced in this paper produce numerical results of the same level of accuracy as the vector potential method (POT-2.5).However, the TE/TM method is at least 2 times faster due to the smaller number of required operations.

B. Numerical examples calculated with the threedimensional TE/TM-ADI2 scheme
Finally, we discuss the results of numerical computations with the fully three-dimensional realization of the TE=TM-ADI2p scheme ( 18)- (25).Two test problems are considered.Before presentation of the numerical results we want to discuss briefly the realization of the scheme in the code.To be able to calculate long structures, we should spare the computer memory and not keep all geometric information.For this purpose we cut the long structure in short blocks, discretize them, and keep the geometric information in the external memory.The information will be loaded only at the instant when the head of the bunch arrives at a geometric block and it will be deleted after the bunch (together with the moving mesh) has passed through the block.
To be able to check the accuracy of the threedimensional realization of the TE/TM scheme we have chosen rotationally symmetric structures for numerical tests.However, in the three-dimensional calculations the symmetry of the structures was not exploited.
In the first example we consider a structure consisting of the 20 TESLA cells [31] bounded by infinite ingoing and outgoing pipes with a diameter of 35 mm.
Figure 9 shows the longitudinal wake potential  for a Gaussian bunch with a rms length 1 mm moving on the axis.The solid line (POT-2.5D)corresponds to the reference solution obtained with the vector potential method [13].The two other lines show results obtained with different mesh resolutions from the TBCI code [32,33] based on the classical Yee's scheme (E=M-2:5D).The oscillations that appear are due to the dispersion error of the Yee's scheme.The gray points represent the result obtained by the three-dimensional scheme ( 18)-( 25) (marked as TE/TM-3D).
It can be seen that the three-dimensional TE/TM-ADI2 scheme produces very accurate results even for the coarse mesh.Indeed, the three-dimensional code uses only 2.5 mesh points per in the longitudinal direction.In the transverse direction the mesh steps are even 3 times bigger.
As the next example we use the round collimator again.Figure 10 demonstrates the wake potential for the collimator with parameters a 30 mm, b 2 mm, c 50 mm, L 200 mm and a Gaussian bunch with a rms length 1 mm.Again the high accuracy of the suggested threedimensional scheme can be seen.
Finally, in the last example we calculate the longitudinal wake potential for the fully three-dimensional rectangular collimator shown in Fig. 5 on the right.Figure 11 demonstrates the wake potential for the collimator with parameters a 30 mm, b 5 mm, d 20 mm, c 50 mm, L 200 mm and a Gaussian bunch with a rms length 5 mm moving on the axis.Figure 11 (left) compares the wake potential on the axis for the rectangular (solid line) and the round (dashed line) collimator.The round collimator has the same geometric parameters and an aperture with a radius of b 5 mm.As is well known, the wake potential of a round collimator does not change in the transverse plane.Quite contrary, for the rectangular collimator a variation of the wake potential in the transverse plane is expected.Indeed, this can be observed in Fig. 11 (right), where the energy gain for a test particle moving at the position s behind the bunch center is shown.

C. Cherenkov radiation calculated with the damping scheme
In the preceding sections the numerical examples for the scheme without damping were shown.In this section we show the damping of the numerical noise on the example suggested in [25].
Figure 12 shows Cherenkov radiation emitted from a particle traveling faster than the local phase velocity of light.In this simulation an electron is moving in a perfectly conducting pipe filled with a dielectric (" 4" 0 ).The tube has a radius of 1 cm.The electron travels at twice the local phase velocity of light.This is a very strenuous test as the motion excites very short wavelengths.A noisy nonphysical wake is seen behind the particle.Damping schema (17), in contrary, produces an accurate and quiet wake.

VI. CONCLUSION
A new fully three-dimensional implicit scheme for the calculation of electromagnetic fields in the vicinity of FIG. 9. Comparison of the wake potentials obtained by different methods for the structure consisting of 20 TESLA cells excited by a Gaussian bunch with 1 mm.The solid line shows the reference solution obtained with the help of the scheme described in [13].The dashed line describes the solution obtained by classical Yee's scheme with mesh resolution of 5 mesh steps per .The dotted line describes the solution obtained by Yee's scheme with 2 times denser resolution in the longitudinal direction.The picture shows coincidence of the reference result (solid line) with the results on the coarse mesh obtained from the 3D TE/TM code (gray points).FIG. 10.Comparison of the wake potentials obtained by different methods for the round collimator excited by a Gaussian bunch with 1 mm.The solid line shows the reference solution obtained with the help of the scheme described in [13].The dashed line shows the solution obtained by Yee's scheme with a mesh resolution of 5 mesh steps per .The dotted line describes the solution obtained by Yee's scheme with a 2 times denser resolution.The picture shows coincidence of the reference result (solid line) with the results on the coarse mesh obtained from the 3D TE/TM code (gray points).relativistic charged bunches was introduced.As shown by several numerical examples, the scheme is able to model curved boundaries with high accuracy and allows for a nondeteriorating calculation of the field solution for very long simulation times.
To develop the new scheme we proceeded as follows: first, we replaced the E=M splitting of Yee's scheme by the TE/TM splitting.This resulted in an implicit scheme requiring the solution of the Crank-Nicholson equation for the two-dimensional scalar wave equation.We then introduced the TE/TM scheme based on the ADI2 splitting method and studied several test examples.In order to avoid reduction of the maximal time step and to obtain a scheme without dispersion in the longitudinal direction, the conformal approach with nondiagonal material matrices was exploited.It requires the application of iterative procedures.However, already the first iteration produces an accurate and stable solution for all considered examples.
The implicit TE/TM scheme assures energy and charge conservation.Usage of damping terms does not violate the charge conservation and allows one to suppress the high frequency noise generated due to the transverse dispersion.The dispersion relation of the damping scheme was analyzed.
The high overall accuracy of the scheme was demonstrated for realistic collimator problems.The scheme allows one to use a moving mesh and thus to calculate wake fields of very short bunches for a range of problems, for which presently available 3D codes experience severe problems.

FIG. 1 .
FIG.1.Charged particle bunch moving through an accelerating structure supplied with infinite pipes.

FIG. 3 .
FIG. 3. Dispersion and damping in the transverse k x ; k y ; k z 1; 0; 0 direction.The left-hand figure shows the dispersion curve for d 0 (dashed curve), d 0:125 (solid curve), and d 0:25 (dot-dashed curve).The right-hand figure presents the damping for different values of the parameter d.

FIG. 4 .
FIG. 4. The left-hand figure shows second order convergence of the TE/TM(0) (solid line) and TE/TM-ADI2(0) schemes for the sphere.The right-hand figure presents conservation of the discrete energy by different methods for 4 4z.

FIG. 7 .
FIG.7.The left-hand figure shows the transverse wake potential for the TESLA cryomodule excited by a Gaussian bunch with a rms width 1 mm as obtained from the reference code[16].In the right-hand figure the solid line shows the difference between the reference potential and the one obtained with the help of the TE/TM(1) scheme, and the dashed line shows likewise the result for the TE/TM-ADI2(1) scheme.

FIG. 11 .
FIG.11.The left-hand figure shows the longitudinal wake potential on the axis for a rectangular (solid line) and a circular (dashed line) collimator and a Gaussian bunch with 5 mm.In the right-hand picture the energy gain for a test particle moving at the position s behind the bunch center is shown.
can be interpreted as energy conservation FIG. 2. Positions of the relativistic charged particle in the finite integration technique grid in different moments of time.The scaled time step is chosen equal to the longitudinal mesh step.TE/TM FIELD SOLVER FOR PARTICLE BEAM . . .
TE/TM FIELD SOLVER FOR PARTICLE BEAM . . .