Calculation of the normal modes of closed waveguides

Full Text

1. Introduction The simplest way to model the phenomena of classical electrodynamics is to use the Maxwell equations describing the electromagnetic field and their subsequent discretization by the finite difference method. The development of computer technology currently allows using the finite difference method directly to discretize the Maxwell equations and conduct numerical studies of applied electrodynamics problems considered in limited domains of space, for example, in a resonator, prism, diffraction grating etc. One of the most widely used methods of this kind is the finite-difference time-domain (FDTD) method, described in detail in a number of textbooks on modern computational electrodynamics [1]. A specific feature of waveguide problems, e.g., the problem of waveguide diffraction, consists in the fact that the electromagnetic fields is to be calculated at a considerable distance from the studied object that scatters the electromagnetic field, which leads to the need for huge amounts of computation in the framework of the FDTD method and its modifications. It should also be added that this method introduces “numerical dispersion”, which leads to errors in determining the phase velocity, and “numerical anisotropy”, in which the wave numbers of waves propagating in different directions in an isotropic medium differ [2], [3]. The study of waveguide problems in the full electromagnetic formulation was initiated by the works of A. N. Tikhonov, A. A. Samarsky, P. E. Krasnushkin, and A. G. Sveshnikov, carried out in the second half of the last century. A. N. Tikhonov and A. A. Samarsky investigated the propagation of electromagnetic waves along a cylinder with a constant simply connected cross section, having perfectly conducting walls and filled with a homogeneous substance. In this work, several fundamental theorems were proved that characterize an arbitrary electromagnetic field in such a waveguide, e.g., the theorem on the field decomposition into transverse electric and transverse magnetic (TE and TM) fields and the theorem on the decomposition of a field into normal modes. These results allowed P.E. Krasnushkin to introduce the concept of a normal waveguide wave or mode, and A. G. Sveshnikov [4], [5] to introduce partial conditions of radiation and strictly mathematically pose the problem of diffraction and normal waves in a waveguide. It should be noted that setting perfect conduction boundary conditions does not limit the scope of the developed methods to the research and design of microwave transmission channels only, since the closed waveguide model is also used to simulate open waveguide systems [6], [7]. When modeling the propagation of guided modes along an open waveguide in the optical range, it is natural to assume that the field at a distance of several wavelengths from the boundary of such a waveguide is zero. Therefore, by placing an open waveguide in a box with perfectly conducting walls, we obtain an approximate model of an open waveguide, as A.G. Sveshnikov first pointed out. The model “open optical waveguide in a box” is a correct mathematical model describing the propagation of waveguide modes, and at the moment it is a correct model describing waveguide diffraction in open optical systems [7]. The limits of applicability of this model can be described quantitatively by comparing the results obtained with different distances of the box walls from the boundary of the waveguide. An obvious drawback of the “waveguide in a box” model is the overestimation of energy flow channeled in the direction of the waveguide axis. This is not essential for modeling the propagation of guided modes, but it is important, e.g., for problems of energy flowing out of such a waveguide through an open end. The main difficulty in the development of the theory of waveguides was the spectral problem for waveguides filled with optically inhomogeneous matter. As far back as the middle of the last century, waveguides with cores became actively used in practice, i.e., cylinders, the filling of which varies across the section and remains constant along the axis. Below we will call this structure a regular waveguide filled with optically inhomogeneous material. Modern technologies in the field of creating new materials and metamaterials are able to give waveguides with almost any distribution of dielectric constant into the hands of practitioners. Moreover, they are increasingly trying to use fractal inserts [8], which means that the scalar model, with its assumption that the dielectric constant changes are small and slow, is becoming less and less adequate. The same conclusion can be made with respect to the study of multicore waveguides, the study of which is focused on providing 5G networks [9]-[11]. The method, which goes back to the works of A. N. Tikhonov and 1. A. Samarsky, was based on the possibility of introducing two potentials for a hollow waveguide, which are now interpreted as the electric and magnetic Borgnis functions. This circumstance fundamentally distinguishes the computational complexity of the spectral problems for hollow waveguides and for waveguides filled with optically inhomogeneous matter. In the first case, the problems are scalar and well-developed methods are applicable to them, which are equally suitable for problems of acoustics and quantum mechanics. In the case of a waveguide filled with an inhomogeneous optical medium, one has to solve the problem numerically in full vector formulation. A generalization of the classical problems of the mathematical theory of waveguides, particularly, the spectral and diffraction problems, to the case of waveguides with variable permittivity is much more complicated and still not fully explored. Turning to a discussion of the results obtained in studying the spectral characteristics of waveguides filled with inhomogeneous matter, we agree to work in a Cartesian coordinate system whose axis coincides with the axis of the waveguide. The problem of finding the normal modes of a regular waveguide filled with optically inhomogeneous matter is as follows. Given are: § the waveguide cross section , § the distribution of and over the cross section ; here and below we assume these functions to take only positive values, § the frequency , and therefore, the wave number = /. It is required to find all values of the parameter ∈ ℂ, for which the Maxwell equations have a nontrivial solution of the form ⃗(, )-, ⃗ (, )-, (1) satisfying the conditions of ideal conductivity of the waveguide walls and the joining conditions at the discontinuities of the permittivity and permeability. The parameter is called the phase constant. Traditionally, this problem is formulated as an eigenvalue problem with respect to three field components. The choice of these three components from the set of six components of the vectors ⃗ and ⃗ can be different, which leads to different formulations of the problem. A. N. Bogolyubov and T. V. Edakina [6] and Frank Schmidt [12], [13] used the components of the vector ⃗, E. Lezar and D. Davidson [14] working in the framework of the FEniCS Project used the components of the vector ,⃗ A. L. Delitsyn [15]-[17] used the components , , . Normal modes of an axially symmetric waveguide with a dielectric core were considered by N.A. Novoselova, S. B. Raevsky and A. A. Titarenko [18], as well as by A. L. Delitsyn and S. I. Kruglov. First, for any of the above approaches, we obtain spectral problems for nonself-adjoint operators in function spaces, which are constructed by analogy with Sobolev spaces, but are much less studied. This greatly complicates the proof of the theorem on the expansion of a monochromatic wave in a waveguide in normal modes, without which it is impossible to proceed to setting the partial radiation conditions in the waveguide diffraction problem [15], [19], [20]. For a hollow waveguide, this theorem was proved by A. N. Tikhonov and A. A. Samarsky as a consequence of Steklov’s theorem, while in the case of a waveguide filled with inhomogeneous medium, one has to use the general Keldysh theorem on the completeness of the system of eigenvectors and adjoined vectors [16], [17], [21]-[23]. At present, the completeness of the system of principal vectors of a waveguide has been proved, however, the question about the possibility to use this system as a basis has already been answered in the affirmative sense only for the case of a circular waveguide, the filling of which depends only on the radius [24], [25]. Numerous questions about the distribution of eigenvalues, the conditions for the existence of multiple eigenvalues have remained unexplored due to the difficulty of the spectral theory for non-self-adjoint operators. Are there cases in which the phase constant of normal modes has both the real and imaginary parts? Are there cases in which adjoined normal modes arise? What is their physical meaning? Numerical experiments do not give an unambiguous answer to these questions. For example, in our experiments performed using the software package presented at the annual Saratov Fall Meeting International Conference in 2017, higher modes appeared to possess complex phase constants. However, numerical methods always calculate higher modes worse than the lower ones, so in those experiments it was not clear whether we discovered a new physical effect, or encountered a computational artifact. Second, the spectral problem has a zero eigenvalue of infinite multiplicity, because of which, when solving the eigenvalue problem numerically, fictitious (“ghost”) modes arise [6], [26]. At the end of their review of the current advances in solving the spectral problem of the waveguide theory, A. N. Bogolyubov and T. V. Edakina [6] wrote: “... The appearance of false modes is perhaps the most difficult issue in solving waveguide problems using finite elements or finite differences in a variational formulation, and, in our opinion, researchers will turn to it more than once in search of the simplest and most economical ways to identify the waves actually propagating in the waveguide.” Currently, there are two ways to deal with ghosts: the use of penalties [6] or the use of mixed finite elements [14], [27]. The main drawback of the penalty method is that although an increase in the parameter characterizing the size of the penalty leads to a decrease in the number of fictitious solutions, the accuracy of calculating the characteristics of true modes decreases. Therefore, since the mid-1990s, the mixed finite element method has been regarded as the only reliable means of combating ghosts [26]. Test examples showed the reliability of this technique, however, in our opinion, the issues of its substantiation as applied to waveguide problems were not given due attention. It should be noted that this issue closely relates to the use of FEM, and not, e.g., the incomplete Galerkin method. Third, the conditions on the walls of the waveguide ⃗ × ⃗ = 0,⃗ ⃗ ⋅ ⃗ = 0 are not classical, and FEA Softwares have no built-in elements for such conditions. It is not clear how precisely these conditions are approximated in published papers. In papers on optical waveguides [6], [12], [13] the authors assume that the field is zero on the walls. From a physical point of view, this assumption is very reasonable, but unfortunately it leads to a conflict with Müller’s theorem on a field equal to zero on an element of an analytic surface [28]. This conflict is removed further, at the stage of introducing penalties, and therefore, the result is a correct mathematical problem. Fourth, in the examples most interesting for applications, the dielectric constant has discontinuities at the interface of different media filling the waveguide. At these interfaces, the electromagnetic fields suffer discontinuities, because of which a necessity arises to approximate discontinuous functions using the FEM. The results of numerical experiments convincingly support the legitimacy of this operation; however, for this situation, theoretically, the effect of discontinuities on convergence has not been studied. The study of this issue is substantially complicated by the fact that the approximation is carried out in non-standard, poorly studied functional spaces. For example, in the works of A.L. Delitsyn mentioned above, embedding theorems were proved, which for Sobolev spaces were established as early as at the beginning of the last century. Finally, the achieved accuracy of the calculations is not high. E. Lezar and D. Davidson compared their results with the results obtained earlier by Jin [29] in the same way for the same waveguide, namely, the rectangular half-filled waveguide. Only the first branch of the dispersion curve coincided with graphic precision, while the Jin’s next three branches merged into two. In our works [30], [31], a previously unknown representation of electromagnetic fields in a waveguide using four potentials was proposed. These potentials do not reduce the number of functions to be determined, i.e., they do not “integrate” Maxwell equations. But even in the case when the permittivity and permeability are described by discontinuous functions, they turn out to be quite smooth functions. The Maple system has developed a symbolic-numerical method for finding normal modes based on a combination of this representation of the field and the incomplete Galerkin method [30], [32], [33]. Comparison of the calculation results with the results obtained using the mixed finite element method was significantly complicated by the lack of a public version of the finite element implementation. The program, written by Lazar and Davidson [14] as part of the FEniCS Project, was only partially published by the authors and has not been updated since 2012, in particular, important changes in the syntax of the project were not taken into account. Y. Yu. Kuziv resumed this program and performed the necessary calculations for comparison, which will be presented in this article. The ultimate goal of our research is to create numerical methods for solving the major problems of the theory of waveguides and to implement them as software packages focused on a wide range of practical problems, from classical issues of microwave transmission to the design of optical waveguides and sensors. At the same time, we strive for ease of implementation of the developed methods in computer algebra systems (Maple, Sage) or in software oriented at the finite element method (FreeFem ++). We will not describe finite element methods for solving the spectral problem, since they are described in detail in [14]. As to our method, we intend to dwell on it in more detail. 2. Method of four potentials Let be a singly connected waveguide cross section, be a segment of the waveguide axis, be the considered time interval, and , - piecewise continuous functions on , taking only positive values. The electromagnetic field in the closed waveguide × × with the filling , will be understood as vector fields ⃗, ⃗, whose components are defined on ( - Γ) × × , (2) under the condition that the sections of ⃗, ⃗ and their partial derivatives in и on for any and are piecewise smooth functions that satisfy 1. the Maxwell equations ⎧rot ⃗ = - ⃗ , { { {rot ⃗ = +⃗, ⎨ {div ⃗ = 0, { ⎩ {div ⃗ = 0, inside the waveguide × × , 2. the conditions of perfect conductivity of the waveguide walls ⃗ × ⃗ = 0, ⃗ ⋅ ⃗ = 0 at regular points of the boundary × × , 3. the joining conditions (3) (4) ⎧{[⃗ × ] = 0,⃗ [⃗ ⋅ ] = 0, (5) {⎨ ⃗ ⃗ ⃗ ⎩[ × ] = 0, [ ⋅ ] = 0 at regular points of the boundary where the filling has a discontinuity Γ × × . In practice, it is important to find at least approximately such a field that satisfies some additional conditions. For example, the spectral problem of the theory of waveguides formulated above consists in finding fields that have all 3 of these properties, and in addition to them having the form (1). In this paper, we focus on the calculation of such fields, leaving aside the theoretical questions of their existence. Assume for brevity that ⟂ ⃗ = ( , , 0), ∇ = (, , 0), and ∇′ = (-, , 0). According to theorem 2 in [31], the transverse components of such field can be always presented in terms of four scalar potentials in the following form ⟂⃗ = ∇ + 1 ∇′, 1 (6) ⃗⟂ = ∇ℎ + ∇′ℎ. In this case the potentials и ℎ are solutions of the Dirichlet problems Δ = -, and { | = 0, Δ 1 ℎ = +, (7) { (8) ℎ| = 0, while the potentials и ℎ are solutions of the Neumann problems ⎧Δ 1 = -, { ∣ ⎨ { (9) = 0, and ⎩ ⎧{Δℎ = -, ℎ (10) ⎨{⎩ ∣ = 0. Here and below for the Laplace operator div(∇) we use the notation Δ. The existence of at least generalised solutions of the Dirichlet problems is apparent, for the Neumann problems we have checked the fulfilment of the conditions for the existence of a solution [31]. For simplicity, let us assume that these solutions are classical, which can be proved under certain additional assumptions of the smoothness of boundary and the considered fields based on the Weyl lemma [34]. Remark 1. In a hollow waveguide the potentials , ℎ, as well as , ℎ are linearly expressed via the Borgnis functions, so that this representation can be considered as the development of Borgnis’ ideas. It should be emphasized that, due to condition 2, the components of electromagnetic fields suffer discontinuities at the interfaces, and the factors in the formulas (6) are specially selected in [31] so that the potentials are smooth functions. In fact, Eq. (6) allows passing from discontinuous variables to smooth ones. This is very convenient from the point of view of the further application of numerical methods, since not all approximation methods are applicable to discontinuous functions, and if applicable, the convergence in the discontinuous case is noticeably worse. Remark 2. Note that any field of the form (6) satisfies theorem 1 of [31] at the interface between media, therefore, the appropriate conditions are automatically fulfilled below. The method of four potentials consists of a transition from field components to four potentials. Let us apply it to the spectral problem. 3. Numerical example: normal modes of a hollow waveguide The simplest test of the operability of programs for calculating normal modes is to calculate the normal modes of a hollow waveguide, studied analytically. Consider a hollow waveguide with the side of a square cross section = 3.5 ⋅ 10-6 [m] without filling, i.e., = = 1. The normal modes, i.e., the solutions of the Maxwell equations, have the form -. For brevity let us denote = 2/ and = /, we will carry out a numerical calculation for = 1014 [Hz]. In a hollow waveguide 2 = 2 - , where is an eigenvalue of Dirichlet or Neumann problem on the section. For the Dirichlet problem 2 = 2 (2 + 2), , = 1, 2, … ; (11) for the Neumann problem 2 = 2 (2 + 2), , = 0, 1, 2, … ; (12) the zero eigenvalue corresponding to constant eigenfunction and zero field; it should be removed from the list of normal modes. For comparison, the relative error is presented ̂ = ∣ - ∣ , (13) for the calculation of the first 10 modes. Here is the phase retardation coefficient calculated using an analytical formula, and ̂ is the same quantity found numerically. Figure 1 presents the errors of mode calculation using finite elements - FEM - with 8 points along each side of the square. The error, as expected, grows with the mode number . For the first two phase retardation coefficients relative error is less than 10-6, for the next four coefficients the relative error is less than 10-4 and for the last four calculated phase retardation coefficients the relative error is less than 4 × 10-3. The error growth is non-monotonic as follows from the Figure 1. Figure 1. Relative error in the calculation of the first 10 normal modes: a comparison of finite element method (FEM) and analytical solution Figure 2 presents the errors of mode calculation using four potential method and incomplete Galerkin method (IGM). The sum, which represent the approximate solution in IGM, consists of 1022 terms. The error growth is non-monotonic with the mode number , the relative error is less than 1.1 × 10-5 for the first ten calculated phase retardation coefficients. Figure 2. Relative error in the calculation of the first 10 normal modes: a comparison of incomplete Galerkin method (IGM) and analytical solution The performed numerical experiments show that both methods of calculating normal modes with such a choice of parameters yield reasonable results. 4. Filled waveguide Now let us consider a square waveguide (1 = 1, 1 = 1) of width = 3.5 ⋅ 10-6 [m] with four square cores of width = 10-6[m] with 2 = 2, 2 = 1 inside, as shown in Figure 3. The phase retardation coefficients for such a structure cannot be calculated analytically, so we can only compare the phase retardation coefficients calculated numerically using FEM and IGM. For comparison, a relative error Δ is presented in Figure 4, where Δ is defined as presented below Δ = ∣ - ∣ (14) in the calculation of the first 10 modes using these methods. Here is the phase retardation coefficient calculated by FEM, and is the same quantity calculated by IGM. The relative error grows with increasing mode number and does not exceed the value 4 × 10-2 (see Figure 4). Figure 3. Cross section of a waveguide with 4 cores Figure 4. Relative error of calculating the first 10 normal modes: a comparison of the results obtained by FEM and IGM 5. Conclusion The results obtained demonstrate the coincidence of the numerical results obtained by two different approaches. The first approach is based on the use of the FEniCS Project libraries and on the use of FEM for an approximate solution of the problem of finding normal modes [14]. The second approach is based on the application of the recently proposed four potential method [30]. Numerical calculation is performed using the incomplete Galerkin method [33]. Both methods are tested on a waveguide structure corresponding to a hollow waveguide, for which phase retardation coefficients are analytically known. The accuracy obtained by both methods is not high, but sufficient for technical calculations. Using an example of a four-core waveguide, for which the phase retardation coefficients are not analytically known, both methods also give quite similar results, differing by less than 4%.

About the authors

Mikhail D. Malykh

Peoples’ Friendship University of Russia (RUDN University)

Author for correspondence.
Email: malykh-md@rudn.ru

Doctor of Physical and Mathematical Sciences, Assistant Professor of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation

Dmitriy V. Divakov

Peoples’ Friendship University of Russia (RUDN University)

Email: divakov-dv@rudn.ru

Candidate of Physical and Mathematical Sciences, Assistant of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation

Alexandre A. Egorov

A. M. Prokhorov General Physics Institute Russian Academy of Sciences

Email: yegorov@kapella.gpi.ru

Doctor of Physical and Mathematical Sciences, Chief Researcher of Department of Oscillations

38, Vavilov St., Moscow 119991, Russian Federation

Yaroslav Yu. Kuziv

Peoples’ Friendship University of Russia (RUDN University)

Email: yaroslav.kuziw@yandex.ru

PhD student of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation

References