Method of Singular Integral Equations for Analysis of Strip Structures and Experimental Confirmation

Abstract

This paper presents a rigorous solution of the Helmholtz equation for regular waveguide structures with the finite sizes of all cross-section elements that may have an arbitrary shape. The solution is based on the theory of Singular Integral Equations (SIE). The SIE method proposed here is used to find a solution to differential equations with a point source. This fundamental solution of the equations is then applied in an integral representation of the general solution for our boundary problem. The integral representation always satisfies the differential equations derived from the Maxwell’s ones and has unknown functions ${\mu }_e$ and ${\mu }_h$ that are determined by the implementation of appropriate boundary conditions. The waveguide structures under consideration may contain homogeneous isotropic materials such as dielectrics, semiconductors, metals, and so forth. The proposed algorithm based on the SIE method also allows us to compute waveguide structures containing materials with high losses. The proposed solution allows us to satisfy all boundary conditions on the contour separating materials with different constitutive parameters and the condition at infinity for open structures as well as the wave equation. In our solution, the longitudinal components of the electric and magnetic fields are expressed in the integral form with the kernel consisting of an unknown function ${\mu }_e$ or ${\mu }_h$ and the Hankel function of the second kind. It is important to note that the above-mentioned integral representation is transformed into the Cauchy type integrals with the density function ${\mu }_e$ or ${\mu }_h$ at certain singular points of the contour of integration. The properties and values of these integrals are known under certain conditions. Contours that limit different materials of waveguide elements are divided into small segments. The number of segments can determine the accuracy of the solution of a problem. We assume for simplicity that the unknown functions ${\mu }_e$ and ${\mu }_h$, which we are looking for, are located in the middle of each segment. After writing down the boundary conditions for the central point of every segment of all contours, we receive a well-conditioned algebraic system of linear equations, by solving which we will define functions ${\mu }_e$ and ${\mu }_h$ that correspond to these central points. Knowing the densities ${\mu }_e$, ${\mu }_h$, it is easy to calculate the dispersion characteristics of the structure as well as the electromagnetic (EM) field distributions inside and outside the structure. The comparison of our calculations by the SIE method with experimental data is also presented in this paper.

1. Introduction

An analysis of practical problems in electrodynamics (implying, the propagation of electromagnetic (EM) waves), aerodynamics (interaction of gases with a moving body), elastodynamics (properties of elastic waves in time), fluid mechanics (for liquids, gases, and plasmas), solid mechanics (especially solids motion and deformation under the action of forces, temperature etc.), fracture mechanics (spreading of cracks in materials) may be to find a solution to a certain differential equation that describes physical processes under some appropriate boundary (and/or initial) conditions.

A boundary value problem (BVP) is a differential equation together with a set of additional limitations, which are called the boundary conditions. In many problems of electrodynamics, as well as in our algorithm, the problem is reduced to a particular case of BVP: an eigenvalue-eigenfunction problem. Under the rigorous solution, we mean the solution that satisfies the differential wave equation and all boundary conditions.

The difficulties in searching for the rigorous solution may occur when a structure contains materials with the specific constitutive relations, for example, characterized by large loss, and/or a structure has an intricate cross–section shape [1,2,3].

Many problems in elasticity (e.g., elastic materials with inclusions) [4,5,6] and crack problems in non-homogeneous materials with stress singularity [7,8,9] of solid mechanics, the electrostatics and magnetostatics (e.g., calculation characteristic parameters of microstrip lines in TEM approximation) [10,11], aerodynamics [12,13], electrodynamical (e.g., dispersion characteristics of 2D (two-dimensional) and 3D (three-dimensional) structures) [14,15], give rise to singular integrals with the Cauchy kernel. Since most Singular Integral Equations (SIE) arising in applications do not have analytical solutions, there is considerable interest in the numerical solution to these problems.

Here we present a simple technique to find a rigorous solution to the wave equation for electrodynamical problems, for example, as the propagation of hybrid EM waves in a waveguide structure consisting of different isotropic materials and elements with complicated configuration, and having limited sizes in their cross-sections. The wave equation, in this case, is the second-order linear partial differential equation. We are using the theory of SIE [16,17,18,19]. It is important to note that the behavior of singular Cauchy type integrals on a contour of integration, contour ends (if the one is open), for corner points of a contour, as well as in the region on the right and left of the integration contour has been carefully researched by mathematicians with a rigorous mathematical proof of theorems for all situations encountered in the applications.

2. Basis of the SIE Method

For finding a solution, we are using the theory of boundary value problems for analytic functions as the Cauchy type integral and its applications to the solution by the SIE method. We will be describing the main properties of the Cauchy type integral in a short form. Those who want to learn a rigorous approach to problem solving should refer to References [16,17,18,19,20,21].

2.1. Sokhotski–Plemelj Formulae

The Sokhotski–Plemelj theorem relates the principal value of the integral over the contour L with the mean-value of the integrals with the contour displaced slightly above and below so that the residue theorem can be applied to those integrals [22]. The integral:

$$\Phi \left(z\right)=\frac{1}{2\pi j}\underset{L}{\int }\frac{\mu \left(t\right)dt}{t-z},$$

(1)

over the contour L in the plane of the complex variable t is called the Cauchy type integral. The kernel for the fixed value t with respect to the variable z is an analytical function everywhere except the point $z=t$, in which there is a pole [14]. For this case, when $z=t$, we change z on $t_0$ in order to underline that in this singular point there is a pole and integral calculates according to Cauchy’s residue theorem. Therefore, the function $\Phi \left(z\right)$ is analytical everywhere except the points of the contour L, there is no single–valued $\Phi (\infty )=0$. The function $\mu \left(t\right)$ is called a density of the integral (1) and the fraction $1/(t-z)$ is the Cauchy kernel of the one. The Cauchy type integral constitutes a function that is the complex analytic one in the entire plane of the complex variable, except for the points of the contour of integration L. The proof of analyticity of Cauchy type integral consists in establishing the possibility of differentiation with respect to the variable z in the integrand $\mu \left(t\right)/(t-z)$.

We denote ${\Phi }^{+}\left(t_0\right)$ having the limit value of integral (1) when the point $t_0$ is on the left side of the contour L, approaching the point z of the contour, that is, from domain $D^{+}$. The limited value takes integral (1) along the line, which at the finite distances from the point $t_0$ coincides with the path of integration L. At a small distance from the point $t_0$, we take integral (1) along the line of a part of the circumference with a small radius r, when the limit $r\to 0$ (Figure 1). The integral along the circumference part is [14]:

$$\underset{ABC}{\int }\frac{dt}{t-t_0}=j(2\pi -\alpha ).$$

(2)

The value $\alpha$ is the angle between the tangents of a contour L in the point $t_0$. When the contour L is a smooth one then $\alpha =\pi$. In case of a corner point (Figure 1a), the angle $\alpha$ can be any of $0\le \alpha \le 2\pi$. According to the definition, the integral that remains is the principal value integral because $\left|A{t}_0\right|=\left|t_{0}C\right|=\xi$.

Therefore, we can write [14]:

$${\Phi }^{+}\left(t_0\right)=\left(1-\frac{\alpha }{2\pi }\right)\mu \left(t_0\right)+\frac{1}{2\pi j}\underset{L}{\int }\frac{\mu \left(t\right)dt}{t-t_0}.$$

(3)

Therefore, the contour L on the $arcABC$ moves around counterclockwise with respect to its center $t_0$ (Figure 1) when point z is approaching point $t_0$ from domain $D^{+}$. When point z is approaching point $t_0$ from the right side of the contour L, that is, from domain $D^{-}$, the corresponding part of the circumference Figure 2, is on the left side of the contour L. That is why the contour L on the $arcABC$ moves around clockwise with respect to its center. Therefore, the integral [14]:

$$\underset{ABC}{\int }\frac{dt}{t-t_0}=-j\alpha ,$$

(4)

where for the smooth point $\alpha =\pi$. We can write for this case:

$${\Phi }^{-}\left(t_0\right)=-\frac{\alpha }{2\pi }\mu \left(t_0\right)+\frac{1}{2\pi j}\underset{L}{\int }\frac{\mu \left(t\right)dt}{t-t_0},$$

(5)

where the second term of expression is understood in the sense of the principal value of singular integral. The Formulaes (3) and (5) are called the Sokhotski–Plemelj (“jump”) ones. The “jump” means that the value of function $\Phi \left(t_0\right)$ undergoes a leap when crossing the line L. We rewrite these formulae for the smooth contour L (Figure 1b and Figure 2b):

$${\Phi }^{\pm }\left(t_0\right)=\pm \frac{1}{2}\mu \left(t_0\right)+\frac{1}{2\pi j}\underset{L}{\int }\frac{\mu \left(t\right)dt}{t-t_0}.$$

(6)

From the Formula (6), we arrive at the expressions:

$${\Phi }^{+}\left(t_0\right)-{\Phi }^{-}\left(t_0\right)=\mu \left(t_0\right),$$

(7)

$${\Phi }^{+}\left(t_0\right)+{\Phi }^{-}\left(t_0\right)=\frac{1}{\pi j}\underset{L}{\int }\frac{\mu \left(t\right)dt}{t-t_0}\equiv 2\Phi \left(t_0\right),$$

(8)

showing that the difference of the limited values is equal to the density $\mu \left(t\right)$. The integral value in the point $z=t_0$ of the contour L coincides with the half–sum of the limited values.

2.2. Solution of Singular Integral Equations

When we solve an electrodynamical problem by the SIE method, the Cauchy type integral will occur as shown in Reference [14]. We will be considering the solution of integrals when $z\ne t_0$. These SIEs can be numerically solved by using an appropriate numerical integration rule [23,24] and the reduction of the equations to a system of linear algebraic equations, either directly or after the reduction of the Cauchy type singular integral equation to an equivalent Fredholm integral equation of the second kind [25]. Here the SIEs are solved numerically according to the following scheme. In the inhomogeneous Fredholm equation of the second kind:

$$f\left(t_0\right)=\mu \left(t_0\right)+\underset{L}{\int }\mu \left(t\right)K(t,t_0)dt,$$

(9)

$$\underset{L}{\int }\mu \left(t\right)K(t,t_0)dt=\sum _{i=1}^n\underset{t_i-\Delta /2}{\overset{t_i+\Delta /2}{\int }}\mu \left(t\right)K(t,t_0)dt$$

(10)

are written as the sum of integrals over the sections of a small length $\Delta$. Here $K(t,t_0)=1/(t-t_0)$ is the Cauchy kernel of the integral. Every such integral, according to the Krylov-Bogolyubov method [19], is written and calculated by this formula:

$$\underset{t_i-\Delta /2}{\overset{t_i+\Delta /2}{\int }}\mu \left(t\right)K(t,t_0)dt\approx \mu \left(t_i\right)\underset{t_i-\Delta /2}{\overset{t_i+\Delta /2}{\int }}K(t,t_0)dt.$$

(11)

We determined in the Equation (9) that the magnitudes $t_0=t_1$, $t_0=t_2$, $t_0=t_3$…$t_0=t_n$. As a result we arrive at a system:

$$f\left(t_i\right)=\mu \left(t_i\right)+\sum _{j=1}^n\mu \left(t_j\right)\underset{t_j-\Delta /2}{\overset{t_j+\Delta /2}{\int }}K(t,t_i)dt, i=1,2,\dots ,n,$$

(12)

of linear integral equations for determining the values of $\mu \left(t_i\right)$ at separate points $t_i$ of the contour L. The integrals in (12) are calculated in the sense of the principal value. Sokhotski–Plemelj formulae are applied to the Cauchy type integrals (1), when index $i=j$ the integral in (12) is omitted and only the first member $\mu \left(t_i\right)$ is left in the right part of expression (12).

In solving an electrodynamical problem, we need to compute the system of linear Equation (12). It is more convenient to write this system in the matrix form [14]:

$$Ax=B,$$

(13)

where x is a column of unknown coefficients $x_i\equiv \mu \left(t_i\right)$, the index i is the number of equations, the magnitude B is a column of free members $b_i=f\left(t_i\right)$, the term A is the matrix of the coefficients:

$$a_{ij}={\delta }_{ij}+\underset{t_j-\Delta /2}{\overset{t_j+\Delta /2}{\int }}K(t,t_i)dt,$$

(14)

where ${\delta }_{ij}$ is the Kronecker delta and, the Kronecker leads to the appearance of a diagonal member of the matrix. The formulae given in this section are suitable for the integration contours L of any shape, which corresponds to the arbitrary configuration of the boundary separating two materials of a waveguide structure with different constitutive parameters.

3. Scalar Wave Equations for Solution of 2D Problems

The SIE method can be used to investigate the dispersion characteristics of main and higher modes (waves) in regular waveguide structures of arbitrary cross–section geometry containing piecewise homogeneous isotropic materials.

The proposed SIE method is based on the SIE theory [14,16,17,18]. Here we solve the problem using Maxwell’s Equations [26]:

$$\nabla \times \overrightarrow{\mathrm{E}}=-{\mu }_0{\mu }_{\mathrm{r}}\frac{\partial \overrightarrow{\mathrm{H}}}{\partial \mathrm{t}},$$

(15)

$$\nabla \times \overrightarrow{\mathrm{H}}={\epsilon }_0{\epsilon }_{\mathrm{r}}\frac{\partial \overrightarrow{\mathrm{E}}}{\partial \mathrm{t}},$$

(16)

where $\overrightarrow{\mathrm{E}}$ is the electric field strength vector and $\overrightarrow{\mathrm{H}}$ is the magnetic field strength vector. $\epsilon ={\epsilon }_{\mathrm{r}}{\epsilon }_0$ is the absolute permittivity and $\mu ={\mu }_{\mathrm{r}}{\mu }_0$ is the absolute permeability. Value ${\epsilon }_{\mathrm{r}}$ is the relative permittivity and value ${\mu }_{\mathrm{r}}$ is the relative permeability of an isotropic medium. The electric and magnetic constants ${\epsilon }_0\approx 8.85\times {10}^{-12}\left(F/m\right)$ and ${\mu }_0\approx 1.26\times {10}^{-6}\left(H/m\right)$ are called the permittivity and permeability of the vacuum.

The dependencies on time t and that on the longitudinal coordinate z are assumed in the form $e^{i(\omega t-hz)}$, where complex value $h=h^{\prime }-i{h}^{″}$ is the longitudinal propagation constant, $h^{\prime }$ is the phase constant and $h^{″}$ is the attenuation constant. In a lossless medium, the propagation constant h is equal to the phase constant $h^{\prime }$. It is important to emphasize that two different kinds of the imaginary unit j and i are involved in our solution. The imaginary unit on coordinates j of the complex plane when the coordinate of a point is $z=x+jy$ and $j^2=-1$. The imaginary unit i related to time in the multiplier that described the field component phase shift in time and $i^2=-1$.

The magnitude $\omega =2\pi f$ is the cyclic signal frequency and i is the imaginary unit related to time in the multiplier $e^{i(\omega t-hz)}$. The transversal components $E_x$, $E_y$, $H_x$, $H_y$ of the EM field are being expressed through the longitudinal components $E_z$, $H_z$ from Maxwell’s equations as follows:

$$E_x=\left({\mu }_0{\mu }_{r}i\omega \frac{\partial H_z}{\partial y}+ih\frac{\partial E_z}{\partial x}\right)/\left(h^2-k^2{\epsilon }_r{\mu }_r\right),$$

(17)

$$E_y=\left(-{\mu }_0{\mu }_{r}i\omega \frac{\partial H_z}{\partial x}+ih\frac{\partial E_z}{\partial y}\right)/\left(h^2-k^2{\epsilon }_r{\mu }_r\right),$$

(18)

$$H_x=\left(-{\epsilon }_0{\epsilon }_{r}i\omega \frac{\partial E_z}{\partial y}+ih\frac{\partial H_z}{\partial x}\right)/\left(h^2-k^2{\epsilon }_r{\mu }_r\right),$$

(19)

$$H_y=\left({\epsilon }_0{\epsilon }_{r}i\omega \frac{\partial E_z}{\partial x}+ih\frac{\partial H_z}{\partial y}\right)/\left(h^2-k^2{\epsilon }_r{\mu }_r\right).$$

(20)

The longitudinal components $E_z$ and $H_z$ satisfy scalar wave equations, which are Helmholtz’s equations:

$$\left({\Delta }_{\perp }+k_{\perp }^2\right)E_z=0,$$

(21)

$$\left({\Delta }_{\perp }+k_{\perp }^2\right)H_z=0,$$

(22)

where $k_{\perp }^2=k^2{\epsilon }_r{\mu }_r-h^2$, $k_{\perp }=\sqrt{k^2{\epsilon }_r{\mu }_r-h^2}$ is the transversal wave number, $h=2\pi /\lambda$ is the longitudinal propagation constant in the waveguide structure, $\lambda$ is the wavelength in the waveguide structure, $k=\omega /c$ is the wavenumber in the vacuum, c is the speed of light in the vacuum, $k=\omega /v$ is the wavenumber in a medium and $v=1/\sqrt{\epsilon \mu }=c/\sqrt{{\epsilon }_{\mathrm{r}}{\mu }_{\mathrm{r}}}$ is the speed of light in an isotropic medium. The operator ${\Delta }_{\perp }={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ is the transversal Laplacian. The rigorous approach to solving the Helmholtz equations and its justification are given in References [27,28,29,30].

4. The Integral Representation for the Solution of Scalar Wave Equations

The fundamental solution of the differential Equations (21) and (22) is the Hankel function of the zeroth order in the cylindrical coordinate system. The problem is formulated as follows. We have in the complex plane a piecewise smooth contour $L_m$, $m=1,2,\dots$. For example, we use a simple structure containing only two contours (Figure 3).

The contour subdivides the plane into two domains, the inner $D_m^{+}$ and the outer $D_m^{-}$ (Figure 3a). These domains, according to the physical problem, are characterized by different parameters: $D_m^{+}$ at $m=1,2$ has constitutive parameters ${\epsilon }_1$, ${\mu }_1$ or ${\epsilon }_2$, ${\mu }_2$ and $D_m^{-}$ has constitutive parameters ${\epsilon }_0$, ${\mu }_0$. The positive direction of going round the contour is when the region $D^{+}$ is on the left side. One has to determine in domains $D^{\pm }$ the solutions of Helmholtz’s Equations (21) and (22), which satisfies the boundary conditions for the tangent components of the magnetic and electric fields:

$$E^{+}{\left|{}_{Lm}=E^{-}\right|}_{Lm},$$

(23)

$$H^{+}{\left|{}_{Lm}=H^{-}\right|}_{Lm}.$$

(24)

In the case of perfect metal, the tangent components of the electric and magnetic fields are equalized to zero on the metal surface. When constructing a solution to the Helmholtz Equations (21) and (22), we were guided by following considerations. We intended to use the properties of the Cauchy-type integral. The fundamental solution of the Helmholtz equation in the two-dimensional case is the Hankel functions of the first and second kinds. Since we are interested in studying only open structures, we choose the Hankel function of the second kind $H_0^{\left(2\right)}$ for the sign of the time dependence $exp\left(i\omega t\right)$ adopted in our solution. The Hankel function $H_0^{\left(2\right)}$ corresponds to the behavior of an EM wave that with distance from the waveguide structure should decay and at infinity becomes equal to zero. In each cross section (z = const) of the waveguide structure, only waves radiating from the source of EM waves are considered. As we can see in the further solution, such sources of outgoing (radiating) waves are the points at which the unknown ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$, ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ functions are located.

The fundamental solution to a linear differential Helmholtz operator is a mathematical concept that generalizes the idea of the Green’s function (using the delta function) without relation to any boundary conditions. The fundamental solution is an analytical solution of the Helmholtz equation in all space inside our waveguide structure and outside it, with the exception of the integration line, which coincides with contours of the structure. In our solution, we get an analog of a two-dimensional dipole source, the density of which is determined by the functions ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$, ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ from the left and right sides of the contour L, where the integral (1) becomes singular and exists only in sense of principal value of Cauchy type integral on the contour L.

We write the solution to the Equations (21) and (22) in the form:

$$E_z\left(\overrightarrow{\mathrm{r}}\right)=\underset{L_m}{\int }{\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)H_0^{\left(2\right)}\left(k_{\perp }{r}^{\prime }\right)ds,$$

(25)

$$H_z\left(\overrightarrow{\mathrm{r}}\right)=\underset{L_m}{\int }{\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)H_0^{\left(2\right)}\left(k_{\perp }{r}^{\prime }\right)ds,$$

(26)

where $E_z\left(\overrightarrow{\mathrm{r}}\right)$ and $H_z\left(\overrightarrow{\mathrm{r}}\right)$ are the longitudinal components of the electric and magnetic fields. Here $\overrightarrow{\mathrm{r}}=\widehat{\mathbf{x}}x+j\widehat{\mathbf{y}}y$ is the radius vector of the point, where the electric and magnetic fields are determined (see Figure 3b), where $\widehat{\mathbf{x}}$, $\widehat{\mathbf{y}}$ are the unit vectors. Here $ds$ is an element (segment) of the contour $L_m$. The magnitudes ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$ and ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ are the unknown functions satisfying Hoelder condition [16]. Here $H_0^{\left(2\right)}\left(k_{\perp }{r}^{\prime }\right)$ is the Hankel function of the zeroth order and the second kind, where $r^{\prime }=\left|{\overrightarrow{\mathrm{r}}}_s-\overrightarrow{\mathrm{r}}\right|$. The magnitude ${\overrightarrow{\mathrm{r}}}_s=\widehat{\mathbf{x}}{x}_s+j\widehat{\mathbf{y}}{y}_s$ is the radius vector of the contour point at which the density functions ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$ and ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ are determined. We agreed earlier that these points are taken in the middle of segments $\Delta$ of contour (see Equation (10)). We numerically investigated the dependence of the final calculation results on the location of the point on the contour segments for which the functions ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$ and ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ are sought and the boundary conditions are satisfied. We came to the conclusion that a position of the point on segments does not matter if this point does not coincide with the angle of the contour (if there is an angle), therefore, for convenience, we chose its position in the center of the segment. Since the structures investigated have the rectangular shapes and corners, we have divided the integration contour into horizontal and vertical subcontours, thereby we have excluded the possibility of finding ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$ and ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ functions at the corner points of the contour L. For this reason, the angle $\alpha =\pi$ in Formula (3) and (5) are simplified and transformed to expression (6) in our case.

Let us now consider the calculations of simple rectangular waveguide structures, which refer to the slot line structure. Figure 4 shows the points of the contours $L_m$, where we satisfy the boundary conditions on the boundary line $L_m$, dividing media with the constitutive parameters ${\epsilon }_1$, ${\mu }_1$ and ${\epsilon }_0$, ${\mu }_0$ as well as with a metal strip.

We have to satisfy boundary conditions between two media with the constitutive parameters ${\epsilon }_0$, ${\mu }_0$ and ${\epsilon }_1$, ${\mu }_1$ along the open contours: $L_1$ which is the curve ABCPRT and $L_2$ which is the line FG. The boundary conditions between media with parameters ${\epsilon }_0$, ${\mu }_0$ and the metal have been satisfied on the open contours: $L_3$ which is GMNA and $L_4$ which is TKLF. The boundary conditions between the medium with parameters ${\epsilon }_1$, ${\mu }_1$ and the metal can be satisfied on open contours $L_5$ and $L_6$ which are TF and GA, correspondingly. We agreed to bypass contours in a counterclockwise direction.

When solving electrodynamic problems, we satisfy boundary conditions. On the boundary dividing the air and the dielectric, we require the equality of the tangent components of electric and magnetic fields:

$$E_z^{+}=E_z^{-},E_s^{+}=E_s^{-},$$

(27)

$$H_z^{+}=H_z^{-},H_s^{+}=H_s^{-}.$$

(28)

On the boundary of perfect metal:

$$E_z^{+}=0,E_s^{+}=0,$$

(29)

where s is one of the transversal coordinates x or y and the coordinate z is the longitudinal one. Magnitudes $E_z^{+}$, $E_s^{+}$ are tangent components of the electric field, $H_z^{+}$ and $H_s^{+}$ are tangent components of the magnetic field related to the domain $D^{+}$. Similar designations to the components of the fields related to domain $D^{-}$ are: $E_z^{-}$, $E_s^{-}$, $H_z^{-}$, $H_s^{-}$. We place the origin of the coordinate system in the bottom left corner of the structure (Figure 4), then the coordinates of each point. Therefore, we have to look for unknown functions ${\mu }_e\left({\overrightarrow{\mathrm{r}}}_s\right)$ and ${\mu }_h\left({\overrightarrow{\mathrm{r}}}_s\right)$ and the location of the point where we equate tangent components $E_z\left(\overrightarrow{\mathrm{r}}\right)$, $E_s\left(\overrightarrow{\mathrm{r}}\right)$, $H_z\left(\overrightarrow{\mathrm{r}}\right)$, $H_s\left(\overrightarrow{\mathrm{r}}\right)$ are determined by the radius vectors ${\overrightarrow{\mathrm{r}}}_s$ and $\overrightarrow{\mathrm{r}}$. The location of the origin of the coordinate system can be chosen arbitrarily.

The singularities in the expressions (27)–(29) for the transversal components $E_s$ and $H_s$ appear from the term $ln{r}^{\prime }$ of the Hankel function $H_0^{\left(2\right)}\left(k_{\perp }{r}^{\prime }\right)$ in the representations (25) and (26). When we satisfy the boundary conditions (23) and (24) and approach the contour $L_m$ in the expressions (25) and (26), the distances tend to zero and the expressions become proportional to the function $\phi \left(\overrightarrow{\mathrm{r}}\right)=Re\left(\phi \left(\overrightarrow{\mathrm{r}}\right)\right)+jIm\left(\phi \left(\overrightarrow{\mathrm{r}}\right)\right)$, which can be written as:

$$\phi \left(\overrightarrow{\mathrm{r}}\right)=\underset{L_m}{\int }\mu \left({\overrightarrow{\mathrm{r}}}_s\right)ln{r}^{\prime }ds.$$

(30)

We see that the Hankel function in the function $\phi \left(\overrightarrow{\mathrm{r}}\right)$ is replaced by the natural logarithm due to small distances when the boundary conditions are satisfied. We differentiate the function $\phi \left(\overrightarrow{\mathrm{r}}\right)$ since under the boundary conditions in the Equations (17)–(20) there are derivatives $\partial E_z/\partial x$, $\partial E_z/\partial y$, $\partial H_z/\partial x$, $\partial H_z/\partial y$. The partial derivatives of the function (30) are:

$$\frac{\partial \phi }{\partial x}=Re\underset{L_m}{\int }\mu \left({\overrightarrow{\mathrm{r}}}_s\right)\frac{ds}{t_s-z} ,$$

(31)

$$\frac{\partial \phi }{\partial y}=-Im\underset{L_m}{\int }\mu \left({\overrightarrow{\mathrm{r}}}_s\right)\frac{ds}{t_s-z} .$$

(32)

These partial derivatives are the real and imaginary parts of the Cauchy type integral. The coordinate $z=x+jy$ corresponds to the electric and magnetic fields, the coordinate $t_s=x_s+j{y}_s$ corresponds to the location of the unknown function, where j is the imaginary unit in space.

Using the Sokhotski–Plemelj Formula (6) we write the limit values of the derivatives (31) and (32) with respect to the tangent $\tau$ and the normal n to the contour:

$$\frac{\partial {\phi }^{\pm }}{\partial \tau }=Re{e}^{j{\theta }_0}\underset{L_m}{\int }\frac{\mu \left({\overrightarrow{\mathrm{r}}}_s\right)ds}{t_s-z_0} ,$$

(33)

$$\frac{\partial {\phi }^{\pm }}{\partial n}=\mp \pi \mu \left(r_0\right)-Im{e}^{j{\theta }_0}\underset{L_m}{\int }\frac{\mu \left({\overrightarrow{\mathrm{r}}}_s\right)ds}{t_s-z_0} .$$

(34)

In these formulae, the contour point $z_0$ with the radius vector ${\overrightarrow{\mathrm{r}}}_0$, can be reached from the left side ($z\to z_0$ from the domain $D^{+}$) or the right side from the domain $D^{-}$. The upper plus sign in the Formulaes (33) and (34) corresponds to the reaching point $z_0$ from the left side. The lower minus sign in the index of the left part of these formulae corresponds to the reaching point $z_0$ from the right. The integrals in the Formulaes (33) and (34) can be understood as the Cauchy principal value integrals. Here ${\theta }_0$ is the angle between the x axis and the tangent to the contour (in the integration direction) at the contour point $z_0=x_0+j{y}_0$. When applying the Krylov–Bogolyubov method, we can obtain transversal field components at the contour points:

$$\begin{array}{cc}\hfill E_x^{\pm }& =\mp \frac{2{\mu }_0{\mu }_r\omega cos\theta }{{\left(k_{\perp }^{\pm }\right)}^2}{\mu }_h^{\pm }\left(s_j\right)-Q_1\left[k_{\perp }^{\pm }\sum _{j=1}^{n_x}{\mu }_h{\left(s_j\right)}^{\pm }\underset{\Delta L_x}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{y_s-y_0}{r^{\prime }}d{x}_s-\sum _{j=1}^{n_y}{\mu }_h{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_y\right]\hfill \\ & +Q_2\left[\sum _{j=1}^{n_x}{\mu }_e{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_x-k_{\perp }^{\pm }\sum _{j=1}^{n_y}{\mu }_e{\left(s_j\right)}^{\pm }\underset{\Delta L_y}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{x_s-x_0}{r^{\prime }}d{y}_s\right],\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill E_y^{\pm }& =\mp \frac{2{\mu }_0{\mu }_r\omega cos\theta }{{\left(k_{\perp }^{\pm }\right)}^2}{\mu }_h^{\pm }\left(s_j\right)-Q_2\left[k_{\perp }^{\pm }\sum _{j=1}^{n_x}{\mu }_e{\left(s_j\right)}^{\pm }\underset{\Delta L_x}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{y_s-y_0}{r^{\prime }}d{x}_s-\sum _{j=1}^{n_y}{\mu }_e{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_y\right]\hfill \\ & -Q_1\left[\sum _{j=1}^{n_x}{\mu }_h{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_x-k_{\perp }^{\pm }\sum _{j=1}^{n_y}{\mu }_h{\left(s_j\right)}^{\pm }\underset{\Delta L_y}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{x_s-x_0}{r^{\prime }}d{y}_s\right],\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill H_x^{\pm }& =\pm \frac{2{\epsilon }_0{\epsilon }_r\omega cos{\theta }_0}{{\left(k_{\perp }^{\pm }\right)}^2}{\mu }_e^{\pm }\left(s_j\right)+Q_3\left[k_{\perp }^{\pm }\sum _{j=1}^{n_x}{\mu }_e{\left(s_j\right)}^{\pm }\underset{\Delta L_x}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{y_s-y_0}{r^{\prime }}d{x}_s-\sum _{j=1}^{n_y}{\mu }_e{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_y\right]\hfill \\ & +Q_2\left[\sum _{j=1}^{n_x}{\mu }_h{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_x-k_{\perp }^{\pm }\sum _{j=1}^{n_y}{\mu }_h{\left(s_j\right)}^{\pm }\underset{\Delta L_y}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{x_s-x_0}{r^{\prime }}d{y}_s\right],\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill H_y^{\pm }& =\pm \frac{2{\epsilon }_0{\epsilon }_r\omega cos{\theta }_0}{{\left(k_{\perp }^{\pm }\right)}^2}{\mu }_e^{\pm }\left(s_j\right)-Q_2\left[k_{\perp }^{\pm }\sum _{j=1}^{n_x}{\mu }_h{\left(s_j\right)}^{\pm }\underset{\Delta L_x}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{y_s-y_0}{r^{\prime }}d{x}_s-\sum _{j=1}^{n_y}{\mu }_h{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_y\right]\hfill \\ & +Q_3\left[\sum _{j=1}^{n_x}{\mu }_e{\left(s_j\right)}^{\pm }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)|\Delta L_x-k_{\perp }^{\pm }\sum _{j=1}^{n_y}{\mu }_e{\left(s_j\right)}^{\pm }\underset{\Delta L_y}{\int }{H}_1^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)\frac{x_s-x_0}{r^{\prime }}d{y}_s\right],\hfill \end{array}$$

(38)

where $Q_1=i{\mu }_0{\mu }_r\omega /{\left(k_{\perp }^{\pm }\right)}^2$, $Q_2=ih/{\left(k_{\perp }^{\pm }\right)}^2$, $Q_3=i{\epsilon }_0{\epsilon }_r\omega /{\left(k_{\perp }^{\pm }\right)}^2$. The magnitudes $n_x$ and $n_y$ are the numbers of sections into which the horizontal and vertical parts of the contour $L_m$$(m=1,2,3,\cdots )$ are divided. The magnitudes $\Delta L_x$ and $\Delta L_y$ are the horizontal and vertical segments of the contour $L_m$ correspondingly. The total number of the segments is $n_x+n_y$. The value $s_j$ is the coordinate of the center point of the segment with the subscript j; this subscript indicates the ordinal number of the segments $\Delta L_x$ and $\Delta L_y$ that varies in the range from 1 to $(n_x+n_y)$. The radius vector ${\overrightarrow{\mathrm{r}}}_s=\widehat{\mathbf{x}}{x}_s+\widehat{\mathbf{y}}{y}_s$ is directed from the origin of the system coordinate to the point $s_j$ where j is in a range from 1 to $(n_x+n_y)$. The EM field components and the values of the density functions ${\mu }_e\left(s_j\right)$ and ${\mu }_h\left(s_j\right)$ are noted in the upper–right corner with the sign corresponding to different domains. The functions ${\mu }_e^{+}\left(s_j\right)$, ${\mu }_h^{+}\left(s_j\right)$, the transversal wave number $k_{\perp }^{+}$ relates to $D^{+}$ domain and functions ${\mu }_e^{-}\left(s_j\right)$, ${\mu }_h^{-}\left(s_j\right)$, the transversal wave number $k_{\perp }^{-}$ relates to $D^{-}$ domain. These functions at the same contour point are different for the field components in the domains $D^{+}$ and $D^{-}$ which are filled with different materials, for this reason ${\mu }_h^{+}\left(s_j\right)\ne {\mu }_h^{-}\left(s_j\right)$ and $k_{\perp }^{+}\ne k_{\perp }^{-}$. The integrals in (35)–(38), having the kernel:

$$H_1^2\left(k_{\perp }{r}^{\prime }\right)=R{H}_1^{\left(2\right)}\left(k_{\perp }{r}^{\prime }\right)+\frac{2i}{\pi }\frac{1}{k_{\perp }{r}^{\prime }}$$

(39)

have been calculated numerically by singling out the singularity analytically. The first term in the expression (39) is the regular part of the function and the second is the singular part. However, first terms in the expressions (35)–(38) contribute only to the diagonal elements of the system as follows from the Sokhotski–Plemelj Formula (6). The computations of electrodynamical problems show that absolute values of diagonal elements by one–two orders exceed the absolute values of non-diagonal elements. The system being well-conditioned leads to the stability of the solution algorithms. The longitudinal field components after applying the Krylov–Bogolyubov method take the form:

$$E_z^{\pm }=\sum _{j=1}^{n_x}{\mu }_e^{\pm }\left(s_j\right)\underset{\Delta L_x}{\int }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)ds+\sum _{j=1}^{n_y}{\mu }_e^{\pm }\left(s_j\right)\underset{\Delta L_y}{\int }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)ds,$$

(40)

$$H_z^{\pm }=\sum _{j=1}^{n_x}{\mu }_h^{\pm }\left(s_j\right)\underset{\Delta L_x}{\int }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)ds+\sum _{j=1}^{n_y}{\mu }_h^{\pm }\left(s_j\right)\underset{\Delta L_y}{\int }{H}_0^{\left(2\right)}\left(k_{\perp }^{\pm }{r}^{\prime }\right)ds.$$

(41)

The system of the algebraic equations for the unknown functions ${\mu }_e^{+}\left(s_j\right)$, ${\mu }_e^{-}\left(s_j\right)$, ${\mu }_h^{+}\left(s_j\right)$, ${\mu }_h^{-}\left(s_j\right)$ obtained from the boundary conditions is homogeneous. The condition of solvability is obtained by equalizing the determinant of the system to zero. The root of the equation obtained determines the propagation constant. After obtaining the propagation constant the determination of the electric field of the waveguide poses no difficulty. For the correctly formulated problem the solution is one–valued and stable with respect to small changes in the coefficients and the contour form [31].

5. Testing the Algorithm by Comparison with Experiments

Figure 5 presents geometry and designations of the slot line with the final dimensions. In our calculations, all sizes w $w_1,w_2,s,l_1,l_2,d$ should be taken as finite values, due to integration along all contours.

Figure 6 and Figure 7 show the dispersion characteristics for the main and higher modes of slot lines (SLs), that is, the dependence of the real part of the longitudinal propagation constant h in SLs on the signal frequency f. The real part of the propagation constant matches here with the phase constant, which equals $h=2\pi /{\lambda }_{SL}$, ${\lambda }_{SL}$ is the wavelength of the certain mode in SL. The dielectric permittivity of the substrate material equals ${\epsilon }_r=11.8$. This permittivity corresponds to the parameters of the high ohmic semiconductor n-Si in experiments; for this reason, in our calculations, we did not take into account the losses. Also, the material of the strips $w_1$ and $w_2$ corresponded to the perfect metal inside of which the electric $\overrightarrow{E}$ and magnetic $\overrightarrow{H}$ fields are assumed to be zero. The property of $\overrightarrow{E}$ and $\overrightarrow{H}$ fields by which they disappear into the metal we used to satisfy the boundary conditions according to equality (29).

The solid lines in Figure 6 and Figure 7 are our calculations and the circles are the experimental results from References [32,33]. Every circle shows the average mean over several measurements of the wavelength ${\lambda }_{SL}$, and the phase constant h was determined by using the average mean value ${\lambda }_{SL}$. The dimensions of the slot wavequide structures were chosen based on the sizes of the ones from articles [32,33] as well as illustrating the ability of the SIE method to explore the higher modes. This is why we study not only ones–mode SLs but also multi–mode SLs (Figure 6b and Figure 7b). All the sizes in the cross-section of the SLs which were used for the measuring and sizes taken in our calculations were the same.

Figure 6 presents the dependence of value h on the signal frequency f for the SL with the size of slot s between two metal strips: $s=2\times {10}^{-4}$ m (Figure 6a) and $s={10}^{-3}$ m (Figure 6b). The size of the slots differ five times. In our case, an increase in the slot size leads to the appearance of higher modes.

The SL sizes (Figure 5) are as follows: the thickness of n-Si substrate is $d={10}^{-3}$ m, the metal strip width $w_1=w_2=2\times {10}^{-3}$ m, the substrate shoulder $l_1=l_2=2\times {10}^{-3}$ m, and the thickness of the metal strips $t=3\times {10}^{-6}$ m. The dispersion characteristics $h=F\left(f\right)$ for the main mode and three higher modes of the SL with substrates protruding out of the metal strip limits as well as different size s are shown in Figure 6. Only the main mode propagates in the SL with a narrow slot (Figure 6a). When comparing Figure 6a with Figure 6b, we see in graphs that the last SL is multi–mode and has very different dispersion dependencies. The last SL has a denser mode spectrum and the one is characterized by the high-frequency cut-off frequencies when the mode curves coincide with each other at higher frequencies. The dispersion curve of the first higher mode merges with the dispersion curve of the main mode at 18 GHz. The dispersion curves of the second and third higher modes merge with the dispersion curve of the main mode at 23 GHz and 25 GHz, correspondingly.

We can conclude, from a comparison of characteristics in Figure 6a,b, that the size of the slot s strongly affects the number of modes that can propagate in the SL. As is well known, the EM field of modes is usually concentrated in the gap s between the strips; for this reason, the size s strongly affects the dispersion characteristics.

Figure 7a,b shows the dispersion characteristics for the main mode and the first higher mode propagating in the SLs with the metal strips reaching to the edge of the n-Si substrate. The figures also present the comparison with the experimental data (circles) for both cases. The common sizes for both SL (see Figure 5) are: $d={10}^{-3}$ m, $s=2\times {10}^{-4}$, $t=3\times {10}^{-6}$ m and $l_1=l_2=0$. The sizes of metal strips are: $w_1=w_2=2\times {10}^{-3}$ m (Figure 7a) as well as $w_1=3\times {10}^{-3}$ m, and $w_2={10}^{-3}$ m (Figure 7b).

Figure 7b shows the dependence of the phase constant on the signal frequency for the SL with the different sizes of the metal strips. When comparing the theoretical and experimental results in Figure 7b, we see that with the frequency increasing, the experimental data moves from the theoretical dispersion curve for the main mode (curve 1) nearer to the curve 2 of the first higher mode. The measurement sensor, which has been moved over the slot s of SL, reacts to the superposition of the main mode and the first higher mode. When the frequency increases, the influence of the higher mode becomes stronger. As a result, the contribution of the higher mode in the experimental signal increased, and the wavelength value recorded by the sensor became to the first higher mode value.

We would like to note that after solving the system of Equation (13), created on the basis of boundary conditions (27)–(29) and the determination of density functions ${\mu }_e\left({\mathrm{r}}_s\right)$ and ${\mu }_h\left({\mathrm{r}}_s\right)$ in the central points of all contour segments, it is easy to calculate all components of the electric and magnetic fields according to the expressions (35)–(38), (40) and (41).

6. Conclusions

We have proposed a simple and universal method for calculations of dispersion characteristics of different regular waveguide structures with the arbitrary configuration and final sizes of cross-section elements containing isotropic dielectric, semiconductor and metal fillings. We compared here our calculations of dispersion characteristics with experimental results. We see the same character of experimental and theoretical dependencies on the phase constant from frequency for different SLs. The difference in these results is within the margin of experimental error.

The novelty of the work lies in the use of the method of the Singular Integral Equations for solving a specific electrodynamical problem and confirmation of the numerical results by comparison with experiments.