Singular Integral Equations in Scattering of Planar Dielectric Waveguide Eigenwaves by the System of Graphene Strips at THz

We consider the scattering of the H-polarized eigenwaves of a planar dielectric waveguide by a coplanar system of graphene strips in the THz range. The strips are placed along the centreline of the waveguide. Our treatment is based on the singular integral equations with the Nystrom-type discretization algorithm. Dependences of the scattering characteristics, near and far fields, are studied. Frequency scanning radiation patterns are presented. Maximum of the radiated power is observed near the plasmon resonances. The resonant frequency and main lobe level can be controlled by variation of the chemical potential. Applied optimization procedure allows to obtain the radiation pattern with the side-lobe level less than − 20 dB. The presented results can be used in designing of graphene leaky-wave antennas.


Introduction
Radiating structures based on gratings embedded into a dielectric waveguide transforming eigenwaves to free space waves are promising as low-cost, low-profile, and easy-tofabricate elements of millimetre-wave devices such as filters or antennas with frequency scanning ability [1,2].
New materials are available today, such as graphene. Graphene can be considered as 2d plane of carbon atoms, ideally 1 atom thick. It is famous for being electrically conductive, mechanically strong, and optically transparent. Graphene strips can support plasmon polariton surface waves and corresponding plasmon resonances at THz. The resonant frequency can be tuned by applying electrostatic or magnetostatic doping [3,4]. This feature makes graphene perspective in designing of tunable THz devices. For example, one can adjust the equivalent impedance of the structure and reduce the reflection. Together with the ability of graphene to absorb the electromagnetic field, it becomes a very promising element of absorbers [5][6][7]. The use of graphene for sensing application is also possible [8][9][10]. For the THz range, a numerical analysis of the plasmon resonances on a single graphene strip and corresponding sensing properties in the case of variation of the refractive-index of the host medium is given in Shapoval and Nosich [11], and for several strips with substrate is given in Nejat and Nozhat [12].
The mutual interaction of graphene strips in the array is smaller than metallic ones because of the short surface plasmon polariton wavelength [13,14]. This effect can be used in antennas to reduce their size. The configuration of tunable leaky-wave antenna with graphene strips embedded into or placed on dielectric slab is discussed in previous studies [15][16][17].
Graphene strips may be considered as zero thickness resistive surfaces with the restriction that the strip width is > 100 nm. Here one can assume that the edge effect on the graphene conductivity is negligible, and the electrical conductivity model developed for infinite graphene surfaces can be used. For the considered frequency range 1…4 THz, permittivity, and chemical potential values, the spatial dispersion can be neglected [18]. In the absence of spatial dispersion and magnetostatic bias field, the conductivity of graphene is a scalar function = (f , c , , T) of frequency f, chemical potential c electron relaxation time and temperature T. It can be obtained from the Kubo formalism [19,20].
Finite-difference time-domain method is one of possible computation instruments for graphene structures [5,[21][22][23][24]. However, as in the case of commercial packages based on the finite-elements method [25], it leads to the huge system of equations due to the large domain of discretization. In addition, the time step is defined by the thickness of graphene which is only one-atom. Another disadvantage is inability to analytically satisfy the radiation conditions that limits the accuracy of the results to several digits. Approximate techniques, which work much faster, can give reliable results and describe physical effects [17,26].
In previous studies [27,28,32], several types of resonances are discussed which graphene strip grating embedded into dielectric slab can support. Near the plasmon resonances, an increase in the total scattering and absorption cross section is observed. Except the plasmon resonances, which can arise on a single graphene strip, the grating inside the dielectric slab can support the grating-mode resonances. It is shown that in the case of plane wave incidence, near these resonances significant increase in scattering and absorption is also observed.
In the present work, we study how the radiation from the dielectric waveguide with finite graphene strip grating can be controlled in the THz range. For the solution, we use the method of SIEs with discretization based on the Nystromtype method of discrete singularities [33,34]. Preliminary results were presented in the conference papers [35,36]; here, results are significantly extended.

Singular Integral Equation
Let us consider planar dielectric waveguide with coplanar graphene strips placed as shown in Fig.1. We denote the set of strips as L = N ∪ n=1 L n , where L n is the n th strip. The width of the waveguide is 2 h; its permittivity is 0 .
We will assume that the field is independent of x and consider only the H-polarized waves. Hence, the electromagnetic field has the following components: (0, E y , E z ), (H x , 0, 0) , which satisfy the Helmholtz equation, the boundary conditions at the dielectric-vacuum interface the boundary conditions on graphene strips and outside the strips the radiation and the edge conditions. The signs " ± " indicate the limit values of the field components from above (below) the interfaces. Outside the dielectric waveguide, the radiated field is the field of cylindrical waves which propagate to the infinity. The amplitude of this field decreases as 1∕ √ k , if k → ∞ , where is the distance. However, inside the dielectric waveguide, −h < z < h , the scattered field can be represented as a sum of the evanescent and propagating eigenwaves. The amplitude of the propagating eigenwaves is non-decreasing function, if k → ∞ . Here, we should use the radiation condition in a special form [37,38].
The incident eigenwave H p of the planar dielectric waveguide with number p propagates from the negative direction of the y-axis. The total field can be represented as a sum of the incident and scattered field, The scattered field can be expanded in terms of Fourier integral in each domain as follows: H s (y, z) = Let us consider the radiation conditions inside the waveguide, −h < z < h Function 1 − g( ) has zeros at the points m , m = 1, 2, … , which correspond to the propagation constants of the eigenwaves of the dielectric waveguide. For the evanescent waves, points m are complex numbers with Im m ≠ 0 . For the propagating waves, points m are real numbers, Im m = 0 . Points m > 0 ( m < 0 ) correspond to the waves which propagate to the positive (negative) direction of the y-axis. As a result, integrands in (5) have singularities on the integration path. They can be treated as poles. After that the Cauchy integral formula or residue theorem can be used. To calculate these integrals, the integration path should be transformed. First of all, we should note that integrands in (5) are meromorphic functions which satisfy the asymptotic relation (5) should give us waves that propagate to the positive (negative) direction of the y-axis or decay exponentially. The integration over the real axis can be exchanged by the integration over the contour in the form of semicircle in the upper or lower half-space that will enclose our initial path. From (12) it follows that to obtain convergent integral for y > 0(y < 0) we should take semicircle in the upper Im > 0(lower Im < 0) half-space. For y > 0(y < 0), points m > 0 ( m < 0) must be inside the contour and points m < 0( m > 0) must be outside the contour. Thus, the integration path should be transformed so it coincides with the real axis everywhere except points m . Negative points m should be bypassed from above, and positive points m should be bypassed from below (see Fig. 2). In this case, the radiation condition is satisfied. Finally, the scattered field inside the waveguide, if y → ±∞ , is The enforcement of the remaining boundary conditions (2) and (3) gives us dual integral equations relatively unknown function b( ): where Z 0 = 120 Ohm. The application of Hilbert transform [39] allows us to reduce the dual integral equations (14) and (15) to SIE over the set of strips where Q(y, ) is the kernel function, F(y) can be expressed in terms of b( ) with the help of The first integral in (16) is understood as Cauchy principal value integral. The kernel-function Q(y, ) (18) in the second integral in (16) contains poles. The integration contour is shown in Fig. 2. However, the integral over infinite semicircle vanishes. Thus, the integration path in the kernel-function coincides with the real axis everywhere, except the poles m , and the poles are bypassed from below. After application of the regularization procedure and Cauchy integral formula, the kernel-function becomes the regular integral. It can be calculated numerically with the use of the quadrature formula such as the Simpson's or Gaussian rules.
For discretization of SIE, the quadrature rule for the Cauchy principal value singular integrals is used [34,35]. Taking into account the edge condition and representing unknown function as a product of new unknown regular function n u n (t) and function with the inverse square root singularity on every strip, F(Ψ n (t)) = u n (t)∕ √ 1 − t 2 , we can obtain system of SIEs on the standard interval (-1;1). Here Ψ n (t) is linear transformation function of (-1,1) to the segment which corresponds to the n th strip L n . After that, integrands are replaced with interpolating polynomials. The nodes coincide with the zeros of the Chebyshev polynomials of the first kind. The values of y, co-called collocation points, are taken from the set of zeros of the Chebyshev polynomials of the second kind. As a result, the set of algebraic equations can be obtained.

Convergence
The solution of (16)-(18) is unique and convergence is based on the theorems [40,41]. To control the error of numerical results, we used the error-function defined as err(M)=|S(M)-S(2M)|/|S(2M)|, where S is equal to the radiated power and M is the order of the Chebyshev polynomial. Numerical results are obtained using C + + .
Let us consider the system of identical strips of a width 2d placed equidistantly. The period is l. Notice that (12) is valid for y > > 1. To calculate the near field, we regularize integrals in (5) in the same way as in (17) and used composed Gaussian quadrature for integrals without singularities. Firstly, it follows from (12) that values T ampl = 1 − 4 a 1 and R ampl = −4 a −1 are reflection and transmission coefficients (amplitude) of the dominant eigenwave of the waveguide. Secondly, amplitude of the field inside the waveguide also equals to the reflection (or transmission) coefficient. Thus, We checked (19) numerically by increasing |y|. Equation (19) is satisfied at the level of machine accuracy. This allows to validate our home-made program block connected with the elimination of singularities in (5) and (17).
In Kaliberda et al. [34,35], we demonstrate that the results of commercial software HFSS agreed well with our results. At the same time, the radiation patterns calculated with the help of HFSS show slight instability: the width and the angle of the main lobe vary in the interval ±2 0 ...3 0 ; the angle and magnitude of the side-lobes significantly depend on the size of the "vacuum box".  the resonances that the structure under study can support. The maxima of the radiation coefficient correspond to the plasmon resonances (marked as P i ). They depend on the parameters of individual graphene strip. The position of these resonances on the frequency axis can be controlled by variation of the chemical potential. The most effective radiation is observed for c = 1eV . For c = 0.3eV near the first plasmon resonance, the radiation efficiency is poor, < 20%. With an increase of the chemical potential, the radiated power also increases near the first plasmon resonance and reaches its maximum value of about 79% for c = 1eV . For greater number of strips N we do not obtain better radiation, if d = 10 m . Almost all nonradiated power is absorbed by the graphene strips, so the radiation efficiency is also about 79%. For c = 1eV near the first plasmon resonance frequency, the reflected and transmitted power does not exceed 2%.As one can know, the natural waves and corresponding resonances can be excited in periodic structures. Except the plasmon resonances, the structure under study can support resonances, which are caused by the periodicity of the displacement of the strips. Such resonances are marked as N qi . The position of these resonances mostly depends on the parameters of the dielectric slab and the period and is slightly perturbed by the parameters of graphene strips. Thus, it cannot be controlled by variation of the chemical potential.

Scattering Characteristics
To identify these resonances, we considered two values of the waveguide width h = 50 m and h = 70 m . These resonances give maxima in the dependences of R and minima in the dependences of Rad. In the case of the plane wave incidence from domain z > 0, resonance N 12 is identified as grating-mode resonance [27,32]. Near the grating-mode resonances, in extremely narrow frequency band, almost total absorption is observed in the infinite graphene grating. These resonances do not arise in the free space. They are also observed in the case of the perfectly electric conducting gratings inside the dielectric slab [42]. The propagation constant χ of natural waves of periodic structure can be obtained from the following equation where t and r are reflection and transmission matrixes of a single strip, I is the identity matrix, and diagonal matrix e has elements exp(ik 1 l m ) . Fig. 7 shows dependences of Im χ as a function of period l. Extremes of Im χ correspond to the resonances of the natural waves of the periodic part of the structure. Figure 8 shows the total field distribution near the plasmon resonances P i . Field has maxima along the strip. The number of maxima is equal to the number of resonance i. As a rule, the first plasmon resonance is much more pronounced. One can observe much more pronounced maxima near P 1 in Figs. 4-6 as well as a high field concentration near the graphene strips. Figure 9 shows the total field distribution near the resonances of the natural waves of periodic part of the structure. N qi The field distribution near the resonances N qi significantly differs from that near the plasmon resonances. It is clearly seen that the field near resonances N qi has different number of variations along the y and z axes inside the waveguide. The total field has q maxima of the amplitude on a single period (l ⋅ n;l ⋅ (n + 1)) and i maxima on the interval 0 < z < h.

Far Field
Let us study the far radiated field for various parameters. Fig. 10 shows the radiation patterns (power) near the first plasmon resonance P 1 and near the resonance of natural wave of periodic structure N 12 . The patterns are normalized by the global maximum which corresponds to the first plasmon resonance, f = 2.46 THz, c = 1eV . We take the parameters of the structure so that only one dominant eigenwave of dielectric waveguide can propagate. The results are presented for the first plasmon resonance P 1 for c = 0.3eV and c = 1eV as well as for the resonance N 12 . If the  parameters of the structure correspond to the plasmon resonance, amplitude of the main lobe reaches its maximum. The strong dependence of the graphene conductivity on the chemical potential allows to tune the antenna resonant frequency and the main lobe magnitude. Despite the large reflection near the resonances of natural waves of periodic structure N qi for practical applications, resonances with even indexes can be interesting. Near N 12 , in-phase excitation of currents on strips is observed. As a result, the angle of the main lobe of the radiations pattern is 90°. As we mentioned above, the position of resonances N qi does not depend on the parameters of graphene strips. However, the amplitude of the main lobe here can be controlled by variation of the chemical potential.
As it is usual for the radiating periodic structures, the considered structure shows the frequency-scanning ability. The angle of the main lobe mostly depends on the normalized period kl. Fig. 11 demonstrates the diffraction patterns with the frequency-dependent main lobe angle for the same value of the chemical potential c = 1eV.The possibility of independent biasing each strip gives additional degrees of freedom. In antenna design, it is required that the radiation pattern has the lowest side-lobe level. By variation of the chemical potential of every individual strip in the array, we are going to reduce the side-lobe level. The results are presented in Figs. 12 and 13 for two values of frequency, f = 2.7THz and f = 3 THz. The values of the strip width and period are the same. Our purpose is to obtain the side-lobe level less that -20 dB. As seen, the actual side-lobe level is in good agreement with the desired one. The values of the chemical potential are given in Table 1. We used the gradient descent algorithm, 150 iterations. A total of -20 dB angular width of the main lobe of the obtained radiation pattern is about 25° for f = 2.7 THz and about 20° for f = 3 THz. The radiation pattern for constant value of the chemical potential c = 1eV is also presented for comparison.

Conclusion
We have studied the scattering of eigenwaves of the planar dielectric waveguide by graphene strip grating in the range from 1 to 4 THz. The scheme of solution is based on the mathematically grounded and effective method of singular integral equations with the discretization by the Nystromtype algorithm. It provides the controllable accuracy of the solution. The calculation time of one curve in Figs.10-12 is about 10 s.
The behaviour of the scattering characteristics of the structure has clear resonant nature. The radiated power and main lobe magnitude is maximal near the plasmon resonances. Variation of the chemical potential allows to tune the resonant frequency and radiated power. Near the resonance of natural waves of periodic part of the structure, the main lobe is perpendicular, but the reflection is high. Every individual graphene strip in the array can be biased separately. In this way, we reduce the side-lobe level.
We believe that presented results can be potentially used in designing of leaky-wave graphene antennas.
Author Contribution Mstislav E. Kaliberda carried out the derivation of equations, wrote the program code and performed computations, and participated in data analysis; Leonid M. Lytvynenko made the problem statement and participated in the derivation of basic equations; Sergey A. Pogarsky participated in data analysis, critically revised the manuscript, and helped draft the manuscript. All authors gave final approval for publication.
Funding This work was supported by the Ministry of Education and Science of Ukraine, grants 0119U002535 and 0117U004964.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code Availability
The home-made code that supports the findings of this study is available from the corresponding author upon reasonable request.