Dispersion properties of rectangularly-corrugated waveguide structures by the in-house 3D FDTD code COCHLEA in cylindrical coordinates

: An in-house full-wave numerical code COCHLEA based on the Finite Difference in Time Domain (FDTD) method has been developed and used to study the dispersion properties of corrugated waveguide structures. The mathematical formulation of the corresponding electromagnetic problem is presented in detail. Numerical results are derived for axially corrugated waveguides similar to those appeared in gyrotron stacked beam tunnels and are compared with those obtained by the in-house semi-analytical numerical code FISHBONE.


Introduction
High-power high-frequency gyrotrons are probably the most promising microwave source for heating and current-drive in modern fusion devices [1]. In the International Thermonuclear Experimental Reactor 1-MW class gyrotrons operating at 170 GHz will be used for Electron Cyclotron Resonance Heating and Electron Cyclotron Current Drive [2].
In gyrotrons, an annular electron beam is generated by a magnetron injection gun and accelerated by the cathode-anode potential difference. Following the converging magnetic field lines of an axially increasing magnetic field, the electrons are guided through the beam tunnel (compression region) to the interaction cavity, where part of their transverse kinetic energy is converted to microwaves. Parasitic oscillations that usually appear in the beam tunnel degrade the quality of the electron beam and consequently the efficient operation of the gyrotron [3]. For this reason, modern gyrotrons use quite complex 'stacked' beam-tunnels consisting of alternating metallic and dielectric rings [4]. Lately, azimuthal corrugations have also been introduced on the metallic rings to further suppress parasitic oscillations [5]. Alternatively, simpler in construction beam tunnels have been proposed and used, like the conical ceramic cylinders made of SiC [6,7] or fully metallic structures with random corrugations of different width and depth [8].
In these variants of beam tunnels, due to their complexity and large electrical size, it is difficult to study the corresponding electromagnetic phenomena. The semi-analytical methods developed to treat them, incorporate assumptions (mainly on the geometrical properties) that limit their applicability to the real complex structure [9][10][11]. On the other hand, numerical electromagnetic methods, e.g. the finite elements [12], the finite integral technique [13], the transmission line matrix [14], the method of moments [15], can in principle handle such complex geometries, requiring though sufficiently large computational resources. In parallel, they do not always provide adequate insight into the underlying physical phenomena.
In this paper, an alternate numerical method, the Finite Difference in Time Domain method (FDTD) [16], is used to study the electromagnetic problem of the corrugated waveguide structure, similar to that appearing in gyrotrons. The Maxwell curl equations are written in the cylindrical coordinate system and properly discretised in time and space, and the fields are calculated on specific spatial and temporal positions [17]. The FDTD method is flexible, but it requires the computational domain to be finite. To avoid undesired reflections at both ends of the structure, absorbing boundary conditions are imposed on them. For this purpose, artificial anisotropic materials are introduced to cause a gradual decay of the incident waves and effectively absorb them. These materials are implemented with several different layers located perpendicular to the propagation axis, the so-called Perfectly Matched Layers (PML) [18]. Since conventional PML presents low absorption for low-frequency and evanescent waves in waveguides [19], we use a Complex Frequency Shifted PML, which is ideal for studying scattering problems in waveguide structures [19,20].
Most of the numerical codes based on the FDTD method use the Cartesian coordinate system because it is simple and the discretisation of the Maxwell equations is straightforward. However, in structures with cylindrical symmetry (like the beam tunnels), this system is not suitable, because the curved boundaries have to be approximated by straight lines (staircase approximation), leading to errors in the calculations reducing the overall accuracy. To avoid such kind of errors, the FDTD in a cylindrical coordinate system has been implemented in the previous years, either in the form of BOR (Bodies of Revolution) [21], generalised coordinate formulation [22] or by following uniform [23,24] and non-uniform [25][26][27] 3D spatial discretisation. In addition, such algorithms can handle anisotropic materials too [24,27].
Calculation of wavenumbers is realised either by solving the corresponding analytical problems [9,10,28,29] or numerically by calculating the resonance frequencies by monitoring multiple spatial points in time domain [30,31]. Usually, these techniques employ periodic boundary conditions to study periodic systems [30].
In this paper, we present the calculation of the dispersion characteristics of rectangularly-corrugated cylindrical waveguides by using a full-wave 3D numerical code based on the FDTD method in the cylindrical coordinate system. Results are postprocessed by applying the spatial fast Fourier transform (sFFT). As a test case (but not limited to) we calculate the dispersion characteristics of finite-length periodic waveguides terminated with absorbing boundary conditions (PML).
The paper is organised as follows: In Section 2, the mathematical formulation is presented. In Section 3, the numerical code implementation is summarised. In Section 4, numerical results are presented for the simple case of an axially corrugated waveguide as well as for the case of a dielectric-loaded corrugated waveguide, both terminated with CFS-PML. These cases model the infinite length corrugated waveguide structure considered in the inhouse code FISHBONE [9,32]. The dispersion properties of the structure are post-calculated by applying a spatial Fourier Transform to the electromagnetic field distributions and compared with those derived by FISHBONE. In Section 5, the conclusions are given.

Mathematical formulation
We solve Maxwell's curl equations, which include sources as well as conducting and dielectric materials, i.e.
with ɛ, μ and σ being the real values of medium's electric permittivity, magnetic permeability and conductivity, respectively. Following the Yee algorithm [16] in the cylindrical coordinate system, the electric and magnetic fields are staggered on the grid, i.e. each field component is located on different grid positions. The latter allows us to apply the centre-differences in order to calculate the spatial derivatives discretely. Time is also discretised to finite time intervals (timesteps). The electromagnetic fields (E and H) are shifted by half timestep. In Fig. 1, the unitary Yee cell in cylindrical coordinates is depicted. The fields are located in appropriate grid positions in order Gauss' laws to be explicitly enforced, as it is explained in detail in [17]. By applying the central difference scheme for the spatial partial derivatives and the leapfrog scheme for the temporal ones, the field update equations are given by (2) and (3) for E r and H r , respectively, while the other field components can be similarly obtained The quantities Δr, Δφ and Δz are the spatial steps on the r, φ and z directions, respectively, the i, j, k are the grid indices referring to the positions (iΔr, jΔφ, kΔz), respectively, and n is the timestep referring to the nΔt temporal position. Constant matrices C, D and F are time-independent and they only depend on the material properties at each point (grid) of the structure. For the grid point with coordinates i, j, k the corresponding matrix elements are given by where ɛ, σ and μ are the electric permittivity, electric conductivity and magnetic permeability, respectively, on the (i, j, k) grid position. To overcome the singularity, which appears at the radial position r = 0 in E φ , Ε z and H r field updates, the discretised Ampere's and Faraday's laws are applied [23]. Note that due to the degeneration of the coordinate system at r = 0, H r and E φ are calculated near r = 0. The respective equations for the field updates on r = 0 + (E φ and H r ) and r = 0 are In order to eliminate the reflections from both ends of the waveguide, two CFS-PML layers, which are located axially at both ends of the structure with an interface perpendicular to the z-axis, are used. The CFS-PML layer formulation is given in [17][18][19]33]. Due to the orientation of the CFS-PML layers in the structure, the stretch of the coordinates will be applied only on the z-axis [19]. The stretching factor is given by where k z is the arbitrary stretching parameter, usually set equal to one. The layer's conductivity σ z and the attenuation factor a z are calculated using the polynomial scaling [17] with the maximum values to be given by the expressions The quantities R(0), n and d denote the desired reflection coefficient on the PML-PEC interface, the wave impedance of the mode under consideration and the depth of the CFS-PML layer, respectively. Note that the polynomial factors m and m a are free parameters. The maximum value of the attenuation factor is given by where f c is the cutoff frequency of the mode under consideration or the geometric mean of cutoff frequencies in case of multiple mode propagation [19]. The values of the spatial steps Δr, Δφ and Δz are also free parameters. Mesh resolution, or cells per wavelength, is given by where λ 0 is the free-space wavelength, D q the size of the q (r, φ, z) dimension of the structure and Δ q the cell size in this dimension. The mesh resolution should be larger than eight cells per wavelength for acceptable accuracy [17]. Note that the spatial step Δr, the axial step Δz and the number of azimuthal cells N φ depend on the frequency. A usual value of the number of cells per wavelength, to sample accurately the spatial dimensions, is about 20 for the r and z dimensions and 10 for the φ dimension.
Furthermore, we calculate the numerical dispersion of the code. Sample runs with variable Δz are performed and the normalised group velocity is calculated for each Δz. A hollow cylindrical waveguide with radius R = 14 mm and length L = 1200 mm is excited with TE 11 mode driven by a Gaussian pulse of central frequency f 0 = 8 GHz and standard deviation σ f = 0.2 GHz for a total simulation time of 10 ns. These characteristics ensure that the pulse will not be reflected from the end of the waveguide and its frequency content is above the cutoff frequency of the excitation mode. The temporal field distribution of E φ is obtained at two different points p 1 = (7 mm, 45°, 12 mm) and p 2 = (7 mm, 45°, 72 mm). The maximum of their cross correlation is calculated and corresponds to the group delay Δt g . Then, by knowing the spatial difference of the probes positions ΔL g the group velocity is calculated as u g = ΔL g /Δt g . In Fig. 2, the normalised phase velocity for a broad range of Δz is depicted. Δz is inverse proportional to the cells per wavelength (12). Radial and azimuthal spatial steps are kept constant at (dr, Nφ) = (1 mm, 64), which correspond to values of cells-per-wavelength 33 and 23, respectively. It is clear that the numerical dispersion has been minimised for cells per wavelength larger than 20.

Numerical code description
Based on the mathematical formulation presented above, the numerical code COCHLEA has been developed in C. In the code, the geometrical characteristics of the structure under consideration are defined by the user in the input file, along with some other important simulation parameters. The geometry consists of a number of stacked rings, each of them described by its radius and initial and final axial positions. The simulation parameters are the simulation time, the spatial steps, the excitation frequency and the modal characteristics of the excitation mode. The positions, where the field distributions are calculated, are also included. The time drive function is a continuous sine with a predefined rise time, different from zero. A sinusoidal function with zero rise time would introduce a signal with a spectrum containing additional undesired frequencies [34]. The code also supports dielectric materials, the characteristics of which are given in the input file. The output of the code is the fields' values on the grid nodes at specific spatial and temporal positions. The main iterator implements the upgrade equations applied in each cell and for every timestep. To reduce the computational time the spatial loops of the field update equations are parallelised using the OpenMP framework. This parallelisation allows us to study electrically large structures.

Numerical results
We study the dispersion characteristics of a periodic waveguide and compare the numerical results to those obtained by the inhouse code FISHBONE [9,32]. Note that in FISHBONE an infinite length corrugated waveguide is considered, whereas in COCHLEA absorbing boundary conditions are imposed at both ends by using CFS-PML. Due to different boundary conditions imposed in the two codes, discrepancies of the numerical results are expected. The characteristics of the geometry under consideration are given in Table 1. Fig. 3 shows the 2D cross section of the structure. The FDTD mesh and simulation characteristics are given in Table 2. For the grid characteristics given in Table 2, the usual runtime is about 8 h on an Intel Xeon E5-2697 v4 (Broadwell) at 2.30 GHz, 128 GB RAM. The simulation time has been appropriately selected in order that a guided mode of the corresponding smooth waveguide (with radius a) with the same frequency with the excitation performs   [17,35]. Each CFS-PML layer has a length of ten cells along the z-axis. The structure is excited by applying the corresponding transverse spatial distribution of a specific mode of the smooth waveguide at z = 0 and on r-φ plane. The distribution is also multiplied by the sinusoidal time driving function F m, p r, φ, t = S m, p r, φ T t (14) with F m,p (r,φ,t) is the excitation distribution function, S(r,φ) is the mode spatial distribution and T(t) is the temporal driving function. By applying the sFFT on the axial field distributions at a specific position (r, φ), we obtain the wavenumber spectrum. This position should be appropriately selected in order to avoid points, where the field is very small (minima of the field distribution) or nullified by default. The peak positions in the sFFT magnitude spectrum correspond to the wavenumbers of the field oscillating in the structure. To determine these peaks more accurately, a spline interpolation fit is applied. We keep only those peaks that are higher than a certain threshold given by where A p is the sFFT magnitude of each peak and A p,max the maximum magnitude of the spectrum. This threshold is empirical and has been proven to be effective for all cases considered in this paper. Peaks with a magnitude below this threshold commonly appear due to the sFFT error noise. We calculate the dispersion characteristics of the geometry given in Table 1 for different values of the total number of periods. The waveguide is excited with the TE 01 mode of the corresponding smooth waveguide with radius a. The wavenumbers calculated by the sFFT are depicted in Fig. 4 along with the corresponding dispersion results obtained by FISHBONE. As expected, for a small number of periods these results are different from those derived by FISHBONE, since the structure is quite short. In addition, the mean relative error of COCHLEA results with respect to those obtained by FISHBONE is depicted in Fig. 5. Obviously, by increasing the number of periods, the difference between the two codes is significantly reduced and results for 100 periods are almost identical. However, the simulation of a structure with such a large length requires higher computational resources. This is due to a large number of mesh cells used and the larger time that the waves need to perform enough round-trips to reach steady state. In Fig. 6, the computational time normalised to the one required for the geometry with five periods is depicted with respect to the number of periods and it is linearly increasing as expected.
Next, a parametric study for the convergence of the results is performed with respect to the axial cell size Δz. For a structure consisting of 100 periods, the mean relative error (with respect to the FISHBONE results) is calculated for Δz = 0.1, 0.25, 0.5, 0.75 and 1.5 mm and presented in Fig. 7. According to (12), Δz is inverse proportional to the cells per wavelength. It is clear that by increasing the value of Δz, the relative error increases as expected, because a larger spatial step leads to a smaller number of cells per wavelength according to (12), thus to less accurate results. Furthermore, in Fig. 8 we plot the mean relative error with the number of cells per wavelength.
Next, we increase the frequency range to identify a larger part of the spectrum. The TE 01 mode is also used as an excitation with different operating frequencies and the wavenumbers (red diamonds) are calculated from the axial distribution of E φ at r = 5.5 mm (Fig. 9). In the same figure, results obtained by FISHBONE are also depicted (grey full circles) with a small discrepancy between them. The field is a superposition of the eigenmodes of the structure. The magnitude of the sFFT corresponds to the contribution of each expansion term (amplitude and wavenumber) in the total field. Fig. 10 shows the magnitude of the sFFT (colour scale) of the same dispersion points as shown in Fig. 9. It is clear that the largest contribution is from expansion terms, for which the corresponding wavenumbers are very close to the dispersion curve of the TE 01 of the smooth waveguide with radius a (dashed line). The same procedure is followed for TM modes. We excite with the TM 01 and Fig. 11 exhibits the dispersion curves, while Fig. 12 depicts the magnitude of the sFFT for each wavenumber. As previously, the largest contribution in the field is from the TM modes of the corresponding smooth waveguide with radius a.
From Figs. 9 and 11, it is obvious that the wavenumbers at regions, where the dispersion relation has almost zero slope, cannot be calculated. In these regions, the group velocity u g is very small,     therefore the waveguide wavelength λ g is very large, which means that the length of the structure and the simulated time should be chosen extremely large to obtain satisfactory convergence. Therefore, the simulations are not performed in such areas. In addition, TE 11 mode of the smooth waveguide is used to excite the structure under consideration and dispersion results are depicted in Fig. 13. As previously, parts of the dispersion relation (with zero slope) cannot be calculated. Since the structure is explicitly excited by TE 11 of the corresponding smooth waveguide, this mode (both the backward and forward) is mainly present in the dispersion curve. The last result is better depicted in Fig. 14, where the wavenumbers with the largest amplitude are closer to the dispersion curve of the TE 11 .
In Fig. 15, the r-φ contour plot of E φ is depicted for z = 15 mm and f = 17.04 GHz and it is clear that the field distribution corresponds to that of TE 11 . Moreover, in Fig. 16, the azimuthal field distribution of E φ is shown near r = 0. The field follows a sinusoidal form, similar to the excitation one, without any singularity.
Next, the same structure is simulated with a lossless dielectric material with relative dielectric constant ε r = 2 added inside the corrugations without any recess. TE 01 mode is excited at the same position as previously, z = 0 mm. Dispersion results are given in Fig. 17 and present good agreement with those obtained by FISHBONE.

Conclusions
Dispersion characteristics of rectangularly-corrugated waveguides are calculated by using the numerical code COCHLEA, which is based on FDTD. The mathematical formulation of the corresponding explicit FDTD scheme is given, along with details on the post-processing steps. The code is developed in C and parallelised with OpenMP. By employing sFFT and spline interpolation on the fields' distribution, we obtain the dispersion characteristics. The numerical results are in a good agreement with those obtained by the in-house code FISHBONE. This code is a first step towards the simulation of more complex waveguides (i.e.    azimuthally corrugated beam tunnels) as well as to the selfconsistent modelling of the beam-wave interaction in such structures.