Radiation of a Charge Exiting Open-Ended Waveguide with Dielectric Filling

We consider a semi-infinite open-ended cylindrical waveguide with uniform dielectric filling placed into collinear infinite vacuum waveguide with larger radius. Electromagnetic field produced by a point charge or Gaussian bunch moving along structure's axis from the dielectric waveguide into the vacuum one is investigated. We utilize the modified residue-calculus technique and obtain rigorous analytical solution of the problem by determining coefficients of mode excitation in each subarea of the structure. Numerical simulations in CST Particle Studio are also performed and an excellent agreement between analytical and simulated results is shown. The main attention is paid to analysis of Cherenkov radiation generated in the inner dielectric waveguide and penetrated into vacuum regions of the outer waveguide. The discussed structure can be used for generation of Terahertz radiation by modulated bunches (bunch trains) by means of high-order Cherenkov modes. In this case, numerical simulations becomes difficult while the developed analytical technique allows for efficient calculation of the radiation characteristics.


INTRODUCTION
In recent years, an essential interest is observed in the area of contemporary sources of Terahertz (THz) radiation based on beam-driven waveguide structures loaded with dielectric. Despite of the fact that both ordinary vacuum THz devices (such as classical backward wave oscillator) are widely available and other mechanisms for THz sources are discussed (see, e.g., Refs. [1][2][3]), beam-driven sources are still extremely attractive due to extraordinary THz radiation peak power [4]. According to this idea, Cherenkov radiation should be generated by well-controlled electron bunch passed through a waveguide structure with dielectric filling and open aperture [5,6]. The electron bunch should be modulated so that a high-order Cherenkov frequency is excited, therefore allowing the use of, for example, mm-sized waveguides for THz generation. Another challenge here is efficient extraction of the radiation from inside the structure into free space. The possibilities for using the nonorthogonal end cut for this purpose were theoretically estimated [7] and experimentally confirmed [8]. Nevertheless, rigorous solution for the electromagnetic (EM) field produced by a charged particle bunch passing from the open-ended circular waveguide with dielectric filling is still missing even in the case of orthogonal end cut. In particular, such a solution is required for determination of the area of applicability of the approximate technique used in [7] and it's possible improvement.
General theory for analysis of radiation from openended waveguide structures was actively developed dur- * s.galyamin@spbu.ru † a.tyuhtin@spbu.ru ing several preceding decades [9][10][11]. Typically, the theory of EM processes for the discussed waveguide discontinuity (open end) was constructed for vacuum case and excitation by single waveguide mode, however, vacuum structures excited by a moving charge were also partially investigated [12][13][14][15][16]. It should be also noted here that analytical approaches becomes essentially more complicated while they deal with the structures containing dielectric inclusions [10,17,18].
In a series of recent papers, we started rigorous investigation of the aforementioned problem on EM field in a circular open-ended waveguide with orthogonal cut and dielectric filling excited by the field of a moving charged particle bunch [19][20][21][22][23]. In these publications, the semiinfinite waveguide was placed (embedded) into collinear infinite vacuum waveguide with larger radius. Therefore the considered structure is further referred to as "embedded" structure. Based on a very good agreement between analytical and simulated (using both COMSOL and CST) results observed in the aforementioned papers one can conclude that convenient analytical technique for solution of the described problem was fully approved.
It should be noted that the closed geometry has several advantages compared to the opened one for theory, simulations and possible experiments. First, closed structure possesses discrete mode spectrum, thus simplifying analytical consideration. Second, finite area of EM field existence allows efficient simulations without the need of extremely large amount of computational resources. Third, real experiments on THz generation from the open end can be conducted in circular vacuum chamber, the latter is described by the outer waveguide in the theoretical model.
However, Cherenkov radiation which is responsible for the aforementioned high-power THz emission and there- fore is of most interest in the structure under consideration was not described in details. In the present paper, we give the detailed analytical solution for EM field generated by a charged particle bunch in the "embedded" structure loaded with dielectric. We apply the modified residue-calculus technique residue-calculus technique [17] and describe penetration of Cherenkov radiation into vacuum regions of the structure. The paper is organized as follows. After the Introduction (Sec. I), we present rigorous solution of the problem (Sec. II). Note that this section contains only final analytical results while details of calculations and intermediate derivations are placed into three appendices (App. A, B and C) succeeding the main text. Section III presents numerical results visualizing the obtained rigorous formulas, simulated results (via CST PS package) and comparison between them. The Conclusion (Sec. IV finishes the paper.

II. ANALYTICAL RESULTS
Geometry of the problem under consideration is shown in Fig. 1. A semi-infinite perfectly conducting circular waveguide with radius b filled with a homogeneous dielectric (ε > 1) is put into a concentric infinite waveguide with radius a > b. The structure is excited by a point charge q moving along z-axis with constant velocity V = V e z = βc e z (c is the light speed in vacuum). Corresponding charge density ρ and current density j = j e z have the form Unless otherwise specified, analytical results presented below correspond to the source (1). These results can be easily generalized for the case of a bunch being infinitesimally thin in xy-plane, similar to (1), but having arbitrary charge distribution η(z−V t) along z (longitudinal) direction. In this case, charge and current densities, ρ b and j b = j b e z are: As can be easily shown, to obtain formulas related to the case of the bunch (2) one should substitute whereη(ω/V ) is the Fourier transform calculated for ξ = ω/V . For example, in the case of Gaussian bunch with the rms half-length σ, and one should substitute In this case, the largest essential frequency in the spectrum ω max is determined so that Gaussian exponential term in (6) results in certain predetermined attenuation for ω = ω max . Typical attenuation (for example, used in CST PS code by default) is −20dB, that is This result in ω max = ω σ √ ln 10 ≈ 1.5ω σ . Further the cylindrical frame r, φ, z (associated with the Cartesian frame shown in Fig. 1) is used. The problem will be solved in the frequency domain, so that Fourier harmonic H ωφ will be determined. Other nonzero field components are calculated as follows: Time-domain field dependencies are calculated using the inverse Fourier transform formulas which can be transformed to the following form [24,25]: On the basis of Eq. (10), it is sufficient to consider only positive frequencies in the spectrum.

A. Incident field
Fourier harmonic of the magnetic component of the incident field has the following form [26]: Here 0,1 are Hankel functions of the first order. Equation (12) represents the total field of a point charge (1) uniformly moving in regular waveguide of radius b filled with dielectric ε. Equation (13) represents the total field of the same charge moving in regular vacuum waveguide with radius a.
For εβ 2 > 1, incident field in the inner dielectric waveguide H (i1) ωφ contains field of Cherenkov radiation (so called "wakefield"). Wakefield is the part of (12) with the discrete frequency spectrum, namely the finite set of real "Cherenkov frequencies" ω Ch l which correspond to real positive poles of the expression (12). These poles are determined by the following equation: where j 0l is the zero of the zero-order Bessel function, J 0 (j 0l )=0, l = 1, 2, . . .. It can be shown (for example, by the limiting process from the case with dissipation taken into account in dielectric) that the integration path in (10) passes real "Cherenkov poles" from above. Therefore, these poles contribute to the incident field only behind the charge, i.e. for ζ = z − V t < 0. Contributions of these poles (residues) can be calculated: where In the case of Gaussian bunch (5) due to vanishing exponential term (6) in the spectrum, high-order Cherenkov frequencies are strongly suppressed, therefore one can obtain monochromatic Cherenkov radiation for long enough bunch (this is also true for arbitrary finite length bunch).
In the geometry under consideration (see Fig. 1), the waveguide with dielectric filling has an open end, therefore Cherenkov radiation generated inside the inner waveguide will penetrate both coaxial part of the structure (area 2 in Fig. 1) and wide vacuum part (area 3 in in Fig. 1). Note that penetration of Cherenkov radiation through simple plane infinite interface between two media accompanying generation of transition radiation was investigated previously [25,[27][28][29]. In the case under consideration, the process of penetration occures due to the diffraction mechanism. This process is of main interest in this paper and it can be described by the theory presented below.

B. Scattered field
The unknown additional (scattered) field propagating from the boundary in the domains 1, 2 and 3 can be presented as standard series over corresponding waveguide modes [10]: Note that the first term in the right-hand side of (19) represents the TEM wave with E ωz = 0, in accordance with (9). Here is the transversal eigenfunction of the coaxial region (area 2 in Fig. 1), χ m > 0 is the solution of the dispersion relation for the area 2, N 0 is the Neumann function. Propagation constants are: zm > 0, m = 1, 2, . . .. For readers' convenience, below we discuss the way for solving the problem under consideration just briefly and present the main resulting formulas only. Rather cumbersome details of calculations needed for deep understanding of the used technique are moved into the Appendices. Corresponding references are given in the text.
Performing matching of the components H ωφ and E ωr for z = 0, and eliminating the r-dependence from the resulting relations, after certain analytical transformations we obtain infinite systems for unknown coefficients {A m }, {B m } and {C n } (n = 0, 1, 2, . . .) of mode decompositions (18), (17) and (19), correspondingly. This procedure is explained in detail in Appendix A, where obtained systems (A10), (A11), (A18) and (A19) are presented. Using the modified residue-calculus technique [10,21], these systems can be solved by constructing specific complexvalued function f (w). This procedure is described in detail in the Appendix B. Finally, the coefficients can be expressed through f (w) as follows: Here is the propagation constant of the area 1 in the case of ε = 1, Function f (w) is determined as follows: The correct construction of the function f (w) (33) is the key point of the residue-calculus technique. The details of this procedure is described in detail in the Appendix B. In particular, one should determine zeros {Γ m } shifted with respect to zeros of vacuum problem {γ (1) zm }: where the set {∆ m } determines the unknown shifts.
Shifted zeros of f (w) are the distinguishing feature of the problem with dielectric and this fact complicates significantly the solution (compared to the vacuum case) because the set {Γ m } is determined for each distinct frequency ω. Since this is connected with iterative solution of a certain complicated nonlinear system (see Eq. (B14) in Appendix B), it is difficult to obtain the field spectrum H ωφ for significant range of frequencies [0, ω max ]. Therefore, it is difficult to calculate full time dependencies for the field components using inverse Fourier transform formulas (10). However, the detailed analysis of the function f (w) shows that it contains the same Cherenkov poles ω Ch l as the incident field (12) does (see the Appendix C). Therefore, Cherenkov radiation penetrated all vacuum areas of the structure which is described by contribution of these poles (residues) can be easily calculated, similar to Eq. (15). For example, Cherenkov radiation penetrated areas 2 and 3 (which is of most interest) can be expressed as follows: ωφ which is presented as infinite series over waveguide modes (Eqs. (18) and (19)). Since we suppose that (40) describes radiation, these series should be truncated to contain only propagating modes for given Cherenkov frequency.
Let us discuss the procedure to calculate the contribution of ω Ch l for given l. Since {A m } and {C n } have the pole for ω = ω Ch l , then Here, the term in the numerator is the residue to be found. If we suppose that dielectric possesses some small dissipation, ε = ε ′ + iε ′′ , then Cherenkov pole also becomes complex: Using numerical procedure described in Appendix B, we calculate H (α) In this way one can calculate contributions of all the essential Cherenkov poles relatively simple and fast. Corresponding examples are represented below in Sec. III.

III. NUMERICAL RESULTS AND DISCUSSION
Here we present numerical results obtained via rigorous formulas of the previous section and results of direct numerical simulation in CST PS R package (with the use of wakefield solver). For simulation, we constructed  the model with finite L v and L d (see Fig. 1) and open boundary conditions for z = −L d , z = L v . Also small finite thickness of the inner waveguide wall d w ≪ a, b was taken into account in simulations, i.e. we supposed that coaxial area 2 is determined by inequality b+d w < r < a, z < 0. The adaptive meshing procedure was utilized to obtain optimal simulation parameters and stable results, this point will be explained below using representative example (see Fig. 4).
First, we clarify the statement concerning the frequency spectrum of the fields in vacuum areas of the structure. Figure 2 shows position of the first seven shifted zeros Γ m and "unshifted" zeros γ (1) zm on the complex plane calculated for three Cherenkov frequencies (14). These results are supplemented by Table I where corresponding numerical values (γ zm excluding Γ l (with the number of Cherenkov frequency under consideration). This Γ l is shifted dramatically so that it becomes purely imaginary while initial γ (1) zl was purely real. Moreover, one can learn from Table I that the equality is fulfilled with high accuracy, because Note that in Eq. (45) we present numerical values with 0.01 relative accuracy, similar to Table I.
As it is shown in Appendix B, Eq. (45) leads to conclusion that all sets of unknown coefficients possess poles for Cherenkov frequencies. Contribution of these poles in vacuum areas of the structure can be calculated using Eq. (43).

A.
Single bunch field Here we present numerical results illustrating the field behavior in different subareas of the structure. For all figures, radius of the inner waveguide is the same, b = 0.25 cm. The structure is excited by single relativistic Gaussian bunch (5).
For Figs. 3 -8, the bunch length σ is chosen so that only the first Cherenkov frequency ω Ch 1 lies within essential part of frequency spectrum [0, ω max ] (7) while higher Cherenkov frequencies are suppressed by attenuating Gaussian term (6). In this case we expect monochromatic both Cherenkov radiation in the inner waveguide and Cherenkov radiation penetrating areas 2 and 3. Figure 3 shows transverse electric field E r from the probe located in the inner waveguide for the case of relatively large permittivity, ε = 10, and a = 0.5 cm. The part of the signal enclosed in the dashed line rectangle (0.4 ns < t < 1.1 ns) should be interpreted as the field of Cherenkov radiation. The Fourier spectrum of this part of the signal shown in the inset of Fig. 3 has a strong peak for frequency 15.45 GHz, this is Cherenkov radiation frequency obtained in the numerical experiment. For shortness, this result for simulated frequency will be referred to as "experimental" result throughout this section.
It should be noted that mentioned peak is used for adaptive meshing procedure in CST simulation: mesh is refined (number of lines per wavelength is increased) until the position of the peak becomes stable, i.e. relative difference in position is less than 0.002 for two consequent passes. Figure 4 illustrates this procedure. It shows typical dependence of the experimental Cherenkov frequency on the number of meshlines per wavelength (this is stan-   dard parameter defining mesh density in CST; wavelength is understood as the minimal wavelength which corresponds to ω max ). As one can see, for rare mesh the experimental frequency is considerably larger compared to the theoretical value (15.31 GHz). For 60 lines per wavelength the relative difference is sufficiently small therefore procedure is stopped, the obtained frequency differs from the theoretical by less then one percent. Comparison between first Cherenkov frequencies, theoretical and experimental, for all structures discussed below is shown in Table II. Theoretical value of Cherenkov frequency does not depend on radius of the outer waveguide a, but this is not the case for simulations due to the change in mesh with change in a. In all considered cases, relative difference between theoretical and experimental values is around 1 percent. As one can see below,    Structure and bunch parameters are the same as in Fig. 3 this small difference matters in comparison of the field behavior. Figure 5 shows CST simulated signal from the probe located in area 2 of the structure with ε = 10 and a = 0.5 cm. Solid (green) line corresponds to the field obtained via simulation in CST PS R code. According to the CST curve, with an increase in time t, a strong peak corresponding to the "image" of the bunch can be seen first (this effect was discussed in details in the case of similar vacuum structure in Ref. [23]). After that, some transition process connected with diffraction radiation occurs. For large enough time (t 0.5ns) we see the stationary harmonic process. Top inset in Fig. 5 shows magnified part of the CST curve compared with theoretical curve corresponding to contribution of the first Cherenkov pole, i.e. summand with l = 1 and α = 2 in (39) (red line). Magnitudes correlate well but due to the difference in frequency the curves diverge for large enough time. If we manually adjust the frequency in analytical formulas, i.e. substitute the analytical Cherenkov frequency with the simulated one (see Table II), we will obtain an excellent coincidence between the curves shown in the bottom inset in Fig. 5. The described frequency substitution is further called "frequency adjustment". Based on presented comparison and the tendency for experimental Cherenkov frequency (see Fig. 4) one can conclude on both correctness of the used analytical approach and stable operation of the simulation code for fine enough mesh. Figure 6 shows similar comparison (for the same struc-ture) but for probe located in the area 3 (wide vacuum waveguide). Again, after the frequency adjustment applied the curves correlate very well. Further for all figures the frequency adjustment will be used by default. Note that for given a even the first mode in area 3 is evanescent therefore magnitude of Cherenkov radiation is extremely small. Figure 7, illustrates the case of ε = 10 and larger radius of the outer waveguide, a = 0.9 cm. In this case, Cherenkov radiation penetrated area 3 is more expressed and again it is described very well by analytical formulas. Figure 8, shows signals from symmetrical probes in coaxial and vacuum waveguide areas for the structure with lower permittivity (ε = 2) and correspondingly higher Cherenkov frequency. Again, one can see that pole contribution calculated theoretically describes Cherenkov radiation penetrated vacuum parts of the structure very well.

B. Bunch train field
Here we illustrate the possibilities of the described approach for calculation of Cherenkov radiation at highorder modes. According to the idea of beam-driven THz source described in Sec. I, THz frequencies can be generated in mm-sized waveguides by charged particle bunches with proper charge modulation, i.e. by bunch trains [6]. If we denote byη G the Fourier spectrum of a single Gaussian bunch η G (5), we obtain from (4): The sequence of 2M + 1 identical Gaussian bunches spaced by L and carrying the same total charge has the  Permittivity ε = 2, other parameters are the same as in Fig. 3 following longitudinal charge distribution: and the following spectrum, in accordance with (4): (48) Figure 9 shows comparison of a single Gaussian bunch Fourier spectrum (46) with the spectrum of a bunch train (48) of 15 (M = 7) identical bunches with spacing L > 2σ. Both functions are calculated for ξ = ω/V . Cherenkov frequencies (14) for a mm-sized waveguide are also shown. Parameters σ and L are chosen so that the bunch train spectrum has the expressed maximum exactly at Cherenkov frequency ω Ch 5 . Therefore, this bunch train excites effectively the 5-th Cherenkov mode with the frequency around 0.1 THz falling in the lower part of THz range. In the same manner, other Cherenkov frequencies can be generated separately.
It should be noted that simulation of EM field produced by such bunch trains is rather complicated in CST PS package. In particular, according to Fig. 4, number of meshlines per wavelength required for adequate convergence should be increased considerably. Another issue here is manual determination of bunch profile corresponding to the discussed bunch train. On the contrary, the presented analytical technique allows computation of the EM field properties relatively simple and fast which is illustrated below by Fig. 10. Figure 10 shows behaviour of E r component of Cherenkov radiation at 5-th Cherenkov frequency, ω Ch 5 ≈ 2π·95GHz, in vacuum regions of the structure. Recall that this radiation is generated in the inner waveguide and penetrated vacuum sections of the structure by means of diffraction mechanism. Figure 10(a) shows E r field in three cross-sections of the coaxial region (area 2 in Fig. 1) while Fig. 10(b) shows E r field in three crosssections of the wide vacuum waveguide (area 3 in Fig. 1). Each thin (green) curve shows the E Ch r as a function of r at a given time moment t and given z. In total, each plot contains 151 curves covering the 1.5ns time range with 0.01ns interval. The highlighted solid (red) curve is the curve which provides the maximum field magnitude over the cross-section. Since at the given Cherenkov frequency both coaxial waveguide and wide vacuum waveguide supports 5 propagating modes, field behaviour is rather complicated. As one can see, maximum field in coaxial region is always on the inner waveguide wall. In the wide waveguide, global field maximum is typically at the first or second local maximum.

IV. CONCLUSION
We have considered radiation produced by single Gaussian bunch exiting the open end of a cylindrical waveguide with uniform dielectric filling in the case where this waveguide is put into concentric vacuum infinite waveguide of larger radius. Based on residue-calculus technique, we have constructed the rigorous theory of the electromagnetic process in this structure. Based on this theory, Cherenkov radiation exiting from dielectric waveguide into vacuum parts of the structure, which is of essential interest in the context of beam driven radiation sources, can be calculated easily and fast. We also have performed numerical simulation of the process in CST PS code. It has been shown that simulated Cherenkov radiation spectral peak has correct frequency for only dense enough mesh. In our simulations, we have reached mesh density around 60 lines per minimal wavelength in the spectrum, resulting in around 1% difference between the-  Figure 10. Cherenkov radiation field (Er component) at 5-th Cherenkov frequency, ω Ch 5 ≈ 2π·95GHz, in vacuum regions of the structure: in the coaxial area (a) and in wide vacuum waveguide (b). Each thin (green) curve shows the E Ch r as a function of r at given time moment t and given z. In total, each plot contains 151 curves covering the 1.5 ns time range with 0.01 ns interval. Solid (red) curve corrresponds to the global field maximum over the cross-section. Parameters of the structure and the bunch train are the same as in Fig. 9 so that the 5-th Cherenkov frequency is effectively generated.
oretical and numerical frequencies. In this case, numerical and analytical results for Cherenkov radiation coincided very well therefore proving both the correctness of rigorous approach and good convergence of numerical procedure.
Moreover, we have considered generation of high-order Cherenkov modes by modulated bunches (bunch trains) in vacuum regions of the structure. Since trains of short bunches generate relatevely high frequencies, correct numerical simulations will require large amount of calculating resources. In this case, the presented rigorous approach allowing convenient analysis of the EM field across the structure will be the preferred method of investigation. As an example, we have calculated spatiotemporal distribution of Cherenkov field at the 5-th Cherenkov frequency (around 0.1THz) generated in vacuum regions of mm-sized embedded structure with dielectric filling of the inner waveguide.
In a similar way, we substitute (18) and (19) into (A2) and (A4), integrate these relations over the interval b < r < a with the weight function rZ p (rχ p ) and utilize the property and formulas analogous to (A6) and (A7). Taking into account that after a series of algebraic manipulations one obtains where n = 0, 1, . . ., Note that the case n = 0 is obtained by integration of (A2) and (A4) over b < r < a without any weight function.
In the issue, we obtain four infinite systems (A10), (A11), (A18) and (A19). These systems can be solved simultaneously using the residue-calculus technique [10,[21][22][23]. This procedure is described in Appendix B. In accordance with the residue-calculus technique, to solve systems (A10), (A11), (A18) and (A19), let us consider the following Cauchy-type integrals over the infinite radius circle C ∞ : where f (w) is a complex-valued function that should be found. These integrals equal zero because we suppose that f (w) vanishes for |w| → ∞. Next step is constructing f (w) so that it has certain specific zeros, poles and behavior for |w| → ∞. To solve this problem, it is useful to have in mind the infinite systems and their solution for the corresponding vacuum problem (with the same geometry and permittivity ε = 1 [23]) and point out the differences. First, in the case under consideration, systems (A10) and (A11) are more complicated while systems (A18) and (A19) are the same. Second, the singularity of the longitudinal electric field near the sharp edge r = b, z → +0, E (3) ωz ∼ 1/z 1/2−τ , sin πτ = (ε − 1)/(2ε + 2), becomes weaker in the presence of dielectric [10,21] (see Fig. 11) compared to the vacuum case where we have E ωz ∼ z −1/2 near this edge. Taking into account these points, one should construct f (w) so that: argument {∆ m } of ∆ p , u p and υ p± . This system can be solved numerically using iteration procedure. Possibility to control the convergence of this procedure is connected with Meixner edge condition (B2).
As it was shown in [19,21], condition (B2) dictates the following asymptotic behavior of coefficient A p for p → ∞: This in turn results in the asymptotic behavior of f (w) determined by condition (iv). Since asymptotic of f (w) is determined by asymptotic of Γ p , γ zn and γ zm for large numbers, the asymptotic of γ (2) zn and γ (3) zm can be easily learned from their definitions (24) and (23), Γ p should behave as follows for p → ∞: Therefore, the iteration process for solving (B14) is organized as follows. We fix quantity N of ∆ m , m = 1, 2, . . . N to be found. For zero-order approximation, we put ∆ m = τ for all m in the right hand side of (B14) and calculate first-order approximation for ∆ p , p = 1, 2, . . . N . Then we substitute these calculated {∆ m } in the right-hand side of (B14) and calculate second-order approximation, etc. After these iterations have converged (relative difference in ∆ N for two consequent steps is within the accuracy), we compare ∆ N with τ : if ∆ N ≈ τ within accepted accuracy, process is stopped, otherwise N and/or accuracy of calculations is changed and procedure repeats.