Numerical Method for Electromagnetic Wave Propagation Problem in a Cylindrical Inhomogeneous Metal Dielectric Waveguiding Structures

The propagation of monochromatic electromagnetic waves in metal circular cylindrical dielectric waveguides filled with inhomogeneous medium is considered. The physical problem is reduced to solving a transmission eigenvalue problem for a system of ordinary differential equations. Spectral parameters of the problem are propagation constants of the waveguide. Numerical results are found with a projection method. The comparison with known exact solutions (for particular values of parameters) is made.


Introduction
A large class of vector electromagnetic problem concerns electromagnetic wave propagation. In radioengineering and electronics the use becomes waveguide structures of complex cross-sections and with inhomogeneous filling require the construction of mathematical models of the propagation of electromagnetic waves in such structures. It becomes necessary to study a new broad class of problems in electrodynamics, characterized by boundaries with delete, give rise to complex geometry, surfaces, uneven dielectric filling and the presence of thin metal ribs (plates) in the structure. The primary goal here is to construct a numerical method to determine the spectrum of normal electromagnetic waves that propagate in such structures.
The study of the wave propagation in a waveguide filled with inhomogeneous medium are arise a boundary eigenvalue problems for systems of elliptic equations with discontinuous coefficients. On the surfaces of discontinuity are set additional conditions, called transmission conditions. In the simplest formulations the spectral parameter is present only in the equations, and is not included in the transmission conditions. However, in rather complex models a spectral parameter is included not only in the equation, but also in the transmission conditions.
This class of problems has been developed for many years [3,7,12]. The main attention was paid to practical outcomes of the calculation of the characteristics of the main mode of waves which has a greatest interest from the physical point of view, as well as several higher modes. Numerical methods for calculating the parameters of various types of waveguide structures are described in the monographs and review papers [1,4,6,15]. However, it should be said that most of the methods applied to homogeneous waveguides, are not common and are difficult to implement and apply for specific inhomogeneous structures.
In this work the wave propagation in inhomogeneous metal-dielectric cylindrical waveguides is studied numerically using the modification of the projection methods. We consider two types of waveguides: a two-layer dielectric waveguide covered by metal where one of the layers is filled with nonlinear medium and a perfectly conducting cylinder covered by a nonlinear dielectric layer. The obtained numerical results are compared with the data available in the linear theory.

Statement of the problem
Consider three-dimensional space R 3 with a cylindrical coordinate system Oρϕz filled with isotropic medium having constant permittivity ε = ε 0 (ε 0 > 0 is the permittivity of free space), and constant permeability µ = µ 0 ( where µ 0 > 0 is the permeability of free space).
A metal dielectric circular cylindrical waveguide Σ filled with inhomogeneous medium is placed parallel to the axis Oz.  In Figure 1 a perfectly conducting cylinder covered by a dielectric layer, known as Goubau line (GL) is shown. In Figure 2 a circular hollow metallic layered waveguide (HW) is shown.
We will consider monochromatic waves where ( · ) T denotes the transpose operation.
where γ is the real propagation constant (spectral parameter of the problem) and m is an angular integer parameter (which assumed to be known). In what follows we often omit the arguments of functions when it does not lead to misunderstanding.

Differential equations
Substituting E and H with components (2.2) into equations (2.1), we obtain where the prime denotes differentiation w.r.a ρ.
Expressing the functions E ρ ,E ϕ , H ρ and H ϕ through E z and H z from the 1st, 2nd, 4th and 5th equation of system (3.1), we find Substituting the expressions for E ρ , E ϕ , H ρ and H ϕ into the 3rd and 6th equations of system (3.1) and introducing the notation Remark 1. For m = 0 system (3.2) splits into two independent equations, which corresponds to two independent wave propagation problems: for TE and TM guided waves. These problems are well studied (analytically and numerically) e.g. in [8,9,10,13,14,16].
For HW solutions (3.3) takes the form where boundedness of the field at zero is taken into account. For GL solutions (3.3) takes a form where the radiation condition at infinity is taken into account.

Transmission conditions and transmission problem
Tangential components of the electromagnetic field are known to be continuous at the interface. In this case the tangential components are E ϕ , E z , H ϕ and H z . Thus we obtain the following transmission conditions for u E and u M is the jump of the limit values of the function at the point s.
The main problem (Problem P m ) is to findγ such that for a fixed value of angular parameter m = 0, there exist non-trivial functions u E (ρ;γ) and u M (ρ;γ) that satisfy system (3.2), transmission conditions (4.1), and inside the inhomogeneous layer they have the form (3.4) for HW or outside the layer the form (3.5) for GL.

Variation formulation
Let us give the variational formulation of the problem P m . Using the first Green's formula, we obtain quvdρ. (5.1) Taking into account the right-hand side of (3.2), we obtain where u * denotes a replacement by the rule u * E = u M . Let us consider the smooth test functions v E and v M .

Remark 2.
We assume that the test functions v E and v M satisfy the following which coincide with conditions for functions u E and u M at the boundary h 0 .
Multiplying the left and right sides of equations (3.2) by the test functions v E and v M , respectively, and applying (5.1) and (5.2), we obtain which hold for any test functions v E and v M . The solution of (5.6) is equivalent to the original problem P m .

Projection method
Using the projection method [2] let us reduce the variational equation (5.6) to a system of algebraic equations. Firstly, split an interval [h 0 , h] into n subintervals with the length l = |h 0 − h|/n. Let us define a set of n subintervals and set of n + 1 subintervals These subintervals we call base finite elements. Denote by a, b, and c an initial point, midpoint and endpoint of subintervals, respectively. In accordance with the scheme of the projection method, it is necessary to introduce basis functions φ i and ψ j in order to approximate the solution. The basis functions are defined on each subinterval Φ i and Ψ j (φ i and ψ j vanishes outside the intervals Φ i and Ψ j , respectively).

The basis functions
The basis functions ψ i defined on Φ i are where j = 3, ..., n and ψ n+1 = 1 c−a (ρ − a). Such defined basis functions takes into account the physical nature of the problem under consideration (see Remark 2 and Figure 3).
We assume an approximate solution with real coefficients α i and β j such that Substituting functions u E and u M with representations (6.1) into the variational equation (5.6), we obtain a system of linear equations with respect to α i and β j (for fixed value of γ) where matrices A and x have the form where used v E = φ i and i = 1, ..., n; where used v M = ψ j and i = 1, ..., n, j = 1, ..., n + 1; where used v E = φ i and i = 1, ..., n, j = 1, ..., n + 1; where used v M = ψ j and j = 1, ..., n + 1. Thus A is a (2n + 1) × (2n + 1) matrix. Let us denote by ∆ the determinant of A Remark 3. If there exists γ = γ such that ∆( γ) = 0, then γ is an approximate spectral parameter of Problem P m . In other words, if an interval [γ, γ] is such that ∆(γ) × ∆(γ) < 0, then this means that there exists γ = γ ∈ [γ, γ] which is an spectral parameter of Problem P m . This value can be calculated with any prescribed accuracy.

Numerical results
Numerical results are obtained with the help of the shooting method for GL (see Figure 1). For the inhomogeneity ε(ρ) of the waveguide the following functions are used to specify the permittivities ε = ε + ρ 20 and ε = ε c + ρ−2 2 in the layer h 0 < ρ < h, where ε c is a positive real constant.
The dielectric constant determined by the expression ε c + ρ 20 defines a weakly inhomogeneous medium, which can be compared with homogeneous. The function ε c + ρ−2 2 defines a filling considerably different from the homogeneous waveguide.
In the figures below the profiles of dielectric permittivity, the spectral parameter γ with respect to angular parameter m and to angular frequency ω, and functions E z and H z are shown.
In Figures 4 and 5 are presented the permittivity profiles. The red and blue curves correspond to inhomogeneous and homogeneous waveguide structures.      Figures 6 and 7 are graphs of the spectral parameter depending on the angular parameter. The color of curves for the homogeneous and inhomogeneous waveguides coincides with the colors for the dielectric permittivity profiles in Figures 4 and 5. As one would expect the graphs of the homogeneous and weakly inhomogeneous waveguides differ only slightly.
Only integer values of the angular parameter m are physically meaningful. But we can investigate this problem numerically for any values of the angular parameter m (for non-integer, too). When m is equal to 0, the problem splits into two well-studied problem -propagation of TE and TM polarized waves. In these, the first and third black dots correspond to the spectral parameters of the TM waves, where as the second black dot is the spectral parameter of TE waves. These graphs show the connection between the TE and TM waves (m = 0) and hybrid waves (m = 0).
For the value of the angular parameter m equal to 1 and 3, the dispersion curves were plotted in Figures 8 and 9. The interior of the dashed region determines the area of the parameter space where the homogeneous problem has a solution.
For the value of the frequency ω = 1 the graphs of the tangential components E z and H z of the electromagnetic field in Figures 10-13 are presented. Color of curves in Figures 10-13 corresponds to the color of spectral parameter in Figures 8 and 9 (points of intersections of vertical dashed line with the dispersion curves). Ez of electric field in case of m = 1, ω = 1. The permittivity is specified by εc + ρ 20 . Figure 11. The tangential components Hz of electric field in case of m = 1, ω = 1. The permittivity is specified by εc + ρ 20 . Figure 12. The tangential components Ez of electric field in case of m = 3, ω = 1. The permittivity is specified by εc + ρ−2 2 . Figure 13. The tangential components Hz of electric field in case of m = 3, ω = 1. The permittivity is specified by εc + ρ−2 2 .
Graphs of the tangential components of the electromagnetic field (Figures 10-13) consistent with the physical formulation of the problem, namely the component E z vanishes on the boundary of the metal, both of components E z and H z are continuous at the interface and decays when ρ → ∞.

Conclusions
The propagation of monochromatic electromagnetic waves in metal circular cylindrical dielectric waveguides filled with inhomogeneous medium was considered. By applying methods of the theory of integral operators, we gave a rigorous mathematical description of the problem, substantiated and implemented a numerical method for its solution.The method allows us to determine approximate eigenvalues with any prescribed accuracy. The approach described in this paper can be applied to other problems, e.g., to multilayered opened waveguides.