The method of lines extension for the analysis of multilayered graphene-loaded structures in cylindrical coordinates.

In this paper the extended method of lines (E-MoL) is proposed for the analysis of multilayer graphene-loaded three dimensional structures in cylindrical coordinates. Accordingly, the impedance and admittance matrices are defined as the ratios of the electric and magnetic fields at each plane of the stack. The impedance and admittance parameters are transformed from the input to the output of the structure through layers and interfaces, from which, the scattering parameters are extracted. It is assumed that there is an anisotropic graphene layer at the interface of two successive layers. The impedance and admittance transformations at the interfaces are extracted in the cylindrical coordinates. Then the impedance and admittance values at all planes of the stack and consequently, the scattering parameters of the whole structure are derived. To validate the presented method, two validation benchmarks are provided at the microwave frequency band. A circular waveguide and a coaxial cable loaded with graphene plates are analyzed and the results are compared with those of CST simulation software which show good accordance. It is observed that the E-MoL, as a semi-analytical semi-numerical method, is much more time-efficient than the CST software numerical procedure.

Introduction

Graphene is a one-atom thick allotrope of carbon arranged in a hexagonal lattice with extraordinary properties[1]. This 2D material is a promising candidate in electronics and electromagnetics and has attracted the researchers’ attention all around the world in the recent years[2,3]. Nowadays, graphene sheets with dimensions up to 30 cm are fabricated, and this has led to the use of graphene in the new microwave and millimeter wave applications[4]. Graphene sheets have remarkable electric properties like controlling their surface conductivity through applying electromagnetic fields or DC voltages. This feature is used to design tunable electronic or electromagnetic devices[5,6]. New nano-sized electronic and electromagnetic components and devices such as antennas, flexible electronic devices, touch screens and high-speed transistors are examples of the states of the art in the graphene world. The unique properties of graphene are due to its gapless electronic band structure[7-10].

Various methods have been used to analyze multilayer three-dimensional graphene-loaded structures. In[11] per unit length equivalent circuits are proposed for the accurate characterization of propagation in different graphene-loaded waveguides. For that purpose, analytical expressions for the graphene conductivity as well as lumped circuit components are calculated through which the propagation of surface waves along spatial graphene waveguides is illustrated. The finite difference time domain (FDTD) method has the advantage of solving transient broadband problems in one simulation. In[12], the FDTD method is exploited to investigate the electromagnetic interaction with graphene by treating it as a surface boundary condition and sampling the computational space. Recently, in[13] the magnetized graphene is characterized by an anisotropic surface boundary condition using FDTD. Physically, the aim of the proposed method in[13] is to incorporate magnetically-biased graphene formulation into electromagnetic wave scattering problems macroscopically by modeling the graphene with a conducting sheet. Li et al.[14] introduce the discontinuous Galerkin time-domain (DGTD) method to characterize the transient response of the graphene from microwaves to terahertz frequencies. Having considered an infinite graphene sheet with its finite surface conductivity[14], deals with resistive boundary conditions for the graphene. Moreover, the method of moments (MoM) is used for the diffraction analysis of graphene nanoribbons in[15]. The authors also try to model non-local effects of graphene by the MoM analysis. In[16] graphene nanoribbons are analyzed using space-domain MoM with the restriction that the method is valid only for low wavenumbers. Lastly, the theory of characteristic modes is used to characterize the input impedance and efficiency of a plasmonic graphene-based antenna[17]. In[18,19], based on integral equations governing the surface current density of magnetically biased graphene, an analytical approach is proposed to analyze an array of graphene ribbons.

The approaches mentioned above have advantages and disadvantages. Pure analytical methods are very time-effective but are able to solve only certain simple structures. On the other hand, fully numerical methods are capable of analyzing complex structures, but generally they are very time consuming. The method of lines (MoL) is a simple, fast, and accurate approach that solves an electromagnetic problem analytically in one coordinate direction and numerically in other coordinate directions[20-22]. This makes MoL especially appropriate to deal with general multilayered structures loaded with complicated boundy conditions of arbitrary shapes, where the wave propagation in the direction perpendicular to the layers stacking is treated analytically. In directions transverse to the layers the structure is discretized and is sloved numerically. MoL can be used to analyze 2D and 3D structures that require 1D and 2D discretizations, respectively[23]. The MoL enjoys the advantages of both analytical and numerical methods while avoiding the shortcomings of each method when modeling multilayered structures. In[24,25], we have proposed MoL for the analysis of 2D structures including graphene plates in Cartesian coordinates. To verify the efficiency of the method, several 2D graphene-loaded structures such as microstrips, striplines and parallel-plate waveguides are presented and analyzed. In[26], we extend the MoL for the analysis of graphene-loaded 3D structures in Cartesian coordinates. Besides, different graphene-based periodic structures again in Cartesian coordinates, are studied in[27] using MoL.

In this paper the MoL is extended in cylindrical coordinates for the analysis of multilayered structures with an arbitrary number of layers loaded with graphene plates of arbitrary shapes. Here the graphene plates are treated as layers represented by surface boundary conditions. To do so, the impedance and admittance transformation formulas through a graphene plate are extracted. Having obtained the transformations, one can compute at every plane of the structure the impedance and admittance values, and subsequently the scattering parameters are derived. For the purpose of validation, the MoL results are compared with those of CST full wave simulation software which are in good accordance. In the first example, we assume a circular waveguide loaded with graphene sheets and perform a parametric analysis in terms of the graphene chemical potential for different numbers of plates. The second example deals with a coaxial cable loaded with annular graphene sheets. The proposed method is general and can be used to analyze generic 3D multilayered stackings in cylindrical coordinates with applications in absorbers and radiating structures. Note that the MoL is a pure classical theory which is applied to Maxwell’s equations and thus it is not able to model quantum mechanical phenomena such as Klein tunneling which occurs in graphene sheets[28,29].

The rest of the paper is organized as follows: in “The method of lines in cylindrical coordinates” section, the E-MoL is presented in its general form for the analysis of multilayer graphene-loaded three-dimensional structures in the cylindrical coordinates. In “Validation of the proposed method” section, the numerical results concerning the reflection and transmission coefficients of a circular waveguide and a coaxial cable loaded with graphene plates are given as validation scenarios using the proposed theory and CST simulation software. Finally, in “Conclusion” section, the conclusions are provided.

The method of lines in cylindrical coordinates

In this section, the generalized transmission line (GTL) equations are first given in cylindrical coordinates to analyze structures that require two-dimensional discretization. Then, the impedance and admittance transformation formulas in a homogenous layer are given. The presence of graphene at the interface of two adjacent layers will result in a new set of boundary conditions through which the impedance and admittance transformation formulas at the interface are extracted. As a result, the impedance and admittance values are derived in each layer of the multilayer structure. Having obtained all impedances and admittances, the electromagnetic fields all over the structure and subsequently the scattering parameters of the whole structure are achievable.

The generalized transmission line equations in cylindrical coordinates

In what follows, we assume that the constitutive parameters are defined by:

First, we present the GTL equations suitable for z direction propagation. The first and second equations of Ampere's law and Faraday's induction law result in[30]:in which the coordinates and dimensions are normalized to the free space wave number according to and . Besides, the transverse field components in Eqs. (2) and (3) are arranged in a way that the inner product of the field vectors with the transverse field components is proportional to the Poynting vector in z-direction. Using the third equations of the Ampere’s Law and Faraday induction’s Law, the longitudinal field components and are obtained as follows: with transverse fields vectors being defined as and . Substituting Eq. (4) in Eq. (2) and Eq. (3), @@the relation between the transverse field components are obtained as[30]:

Then, we combine the two first-order differential equation systems in Eqs. (5) and (6) in order to obtain the second-order differential equations for the transverse electric or magnetic fields as:

These equations are called as wave equations. To solve these equations using MoL, the fields and their derivatives in transverse directions should be discretized.

Discretization of the fields and solutions

In this section, we present the discretization of the fields components in the cross section as well as the solution of the GTL equations in the discretized form. The discretization scheme is shown in Fig. 1. In what follows, the 2D discretization quantities are marked with a hat . The discretized field components are collected in vectors. To do so, the discretization is started from the waveguide’s center and continue in the radial direction. The discretized constitutive parameters and radii form diagonal matrices. The difference operators are replaced with central differences in radial and azimuthal directions and form the entries of matrices and , respectively. To obtain the difference operators in the 2D discretization scheme, we use the Kronecker product of difference matrix operators in radial and azimuthal directions, i.e. and , and unit matrices, as follows[30]: in which and are unit matrices of sizes (the number of discretization points in radial direction) and (the number of discretization points in azimuthal direction), respectively. Besides, the difference matrix operators in the radial direction for the circle and dot points are constructed using the Dirichlet (Neumann) and Neumann (Dirichlet) boundary conditions in case of magnetic (electric) walls, respectively. For , the periodic boundary condition is employed in case the complete discretization of the cross section in azimuthal direction is required. To continue with the discretization of GTL equations we define[30]:

Figure 1

The discretization scheme in the cross section of a cylindrical waveguide with inhomogeneous layers.

and use the following abbreviations:

As a result, Eqs. (5)–(8) are discretized as follows:in which

In Eqs. (15) and (16), represents the interpolation matrices detailed in the following sections. Since the fields components are discretized in different positions they cannot be added or subtracted directly. Hence, interpolation matrices are used to compute the field components in the correct positions.

Now we are going to solve Eq. (14). For this purpose, we define the electric and magnetic fields as and , in which and correspond to the electric and magnetic field distributions and are defined as the eigenvectors of and , respectively. Besides, and are the electric and magnetic field amplitude vectors.

Introducing and into Eq. (14), one obtains: with

Since , the following relations hold between the eigenvector and eigenvalue matrices of and :

The general solution of Eq. (18) can be given as forward- and backward-propagating waves:where . Also, due to the normalization introduced in Eq. (20), the characteristic impedance and admittance matrices are equal to identity matrices.

Having known the electric and magnetic fields in each plane according to Eq. (21), the relation between the electric and magnetic fields at two different planes with distance d, are obtained as[30]:with

On the other hand, the electric and magnetic fields at each plane are related to each other through impedance/admittance matrices as follows:

Consequently, the impedance and admittance transformation formulas between planes A and B are as follows:

Impedance and admittance transformation formulas at the graphene interface

To analyze the whole structure, it is necessary to be able to pass through the graphene plates at the interface of two adjacent layers by transforming the impedances and admittances through the graphene interface. As shown in Fig. 2, we assume that there is a graphene plate between layers and . This common interface is denoted by the letter B. Besides, the two sides of this interface are distinguished by + and − indices.

Figure 2

graphene-loaded multilayer structure.

As shown in Fig. 2, the graphene is placed in or equivalently plane in which the GTL equations are discretized. Due to the semi-analytical nature of MoL, the equations are solved analytically in direction. For this reason, the analytical formulas must be extracted for the impedance and admittance transformation in this direction. To do so, we must relate the admittance (impedance) on both sides of the graphene interface, i.e. and ( and . The electric current components on the graphene plate are defined by

This equation must be transformed to the cylindrical coordinates system, which will be performed as below. The transverse electric field components must be continuous at the interface. Hence,

The transverse magnetic field components in cylindrical coordinates at are defined as

The boundary conditions for the magnetic field components at plane are

Then, transformation of defined by Eq. (26) into cylindrical coordinates results inwhere the graphene conductivity tensor is obtained in the cylindrical coordinates as

Then, the and components in Eq. (30) are separated as follows

Defining the electric and magnetic field components in cylindrical coordinates asone can rewrite Eq. (32) in the vector form:

In this step, we need to discretize the field components and graphene parameters in their appropriate places as shown in Fig. 1. The discretized form of Eq. (32) is:

As mentioned earlier, the field components, for example and are discretized in different positions and hence we cannot add or subtract them directly to obtain Hence, we need to use interpolation matrices to calculate the fields in the appropriate position. Accordingly, the discretized version of Eq. (34) can be expressed as

interpolates at ° points and calculates its corresponding values at • points. Similarly, interpolates at • points and calculates its corresponding values at ° points.

According to Eq. (33), we rewrite Eq. (38) as

Substituting and in Eq. (39) results in

Multiplying Eq. (40) by , we obtainMoreover, according to Eq. (24) the relation between the electric and magnetic fields is calculated by , substituting it in Eq. (41) results in

Finally, the admittance matrix at is derived as

Thus, transformation of admittance at graphene interface will be done using Eq. (43). Note that if the graphene does not cover the whole cross-section, Eq. (43) can still be used for admittance transformation. However, in those parts of the cross section where graphene is absent, the conductivity parameters are set to zero.

Calculation of the scattering parameters

Reflection coefficient derivation

Starting from the output of the multilayer structure and using the formulas provided for impedance and admittance transformation through layers and interfaces, the input impedance and admittance can be obtained. Usually the load impedance is considered as the output wave impedance of the output waveguide which is normalized to be the identity matrix, i.e. . At the input of the structure there is an input waveguide with a wave impedance or a wave admittance .

Using the admittance transformation formulas, the input admittance is calculated. By dividing the electric field into two parts, the forward and backward propagation components at the input of the waveguide are obtained as

Consequently, for the scattering parameter we have

If the modes are ordered properly, the first element in the matrix gives the fandamental mode reflection coefficient.

Transmission coefficient derivation

In order to calculate the transmission coefficient, the field coupled to the output port should be divided into the injected field at the input port. To do this, the fields should be calculated at the whole structure and they must be transformed in a numerically stable manner. The general solution of the wave equation in the Cartesian coordinates can be written as Eq. (21). The indices and denote the forward and backward propagations, respectively. The electric and magnetic fields are written as

The forward and backward parts are connected by the characteristic impedance

Then, the forward and backward parts can be obtained using the total fields

Using the definition of the impedance at any arbitrary plane, the above equations are written asAssume that the field at the plane A is known. To avoid numerical instability, the following steps should be taken in order to calculate the fields at the plane B:

Using Eqs. (50) and (51), the vectors and are computed from the vectors and .

Using the vectors and , the vectors and are calculated as

The total field at the plane B is computed as

Then, supposing that the graphene plate is placed between the layers, the transmission coefficient is obtained. First, we assume that a graphene plate is located between the two layers of the structure, as shown in Fig. 3.

Figure 3

Two-layered structure with a graphene plate at the interface.

The values of the admittances and impedances shown in Fig. 3, are determined by the aforementioned admittance and impedance transformation formula Eq. (43). Also, the values and are calculated from Eq. (25). Assuming that these values are known, the transformation coefficient is derived. Using Eq. (53), the total field at this plane is written as follows

The electric field at the interface of the plane A is continuous and hence . Due to normalization, the impedance and admittance values at the plane are equal to the characteristic impedance of the waveguide . From Eq. (48), the values of the forward and backward fields at the plane are obtained as

Therefore, we haveConsidering that the applied field at the port is equal to , then the forward propagating field (towards direction) at the plane will become , hence

Then, the transmission coefficient matrix is calculated as

The first element of the matrix denotes the transmission coefficient of the fundamental mode, if the modes are ordered properly.

As shown in Fig. 4, there are two layers of graphene inside the rectangular waveguide. The corresponding transmission coefficient is derived similarly. Assume that the applied field to the port C is , then the forward electric field at the plane is calculated as . Considering the continuity of the tangential electric field at the plane B, we have

Figure 4

Three-layered structure with graphene plates at the interfaces.

Therefore, the forward and backward fields at the plane are obtained asThe forward field at the plane is written as

Consequently, using Eq. (53), the total electric field at the plane can be written as

Considering the continuity of the electric field components at the plane we have

Since the electric field at the plane is completely propagating forward, from Eqs. (61)–(63) we can write

By substitution of Eqs. (60) and (52) in (64)

As a result, from Eq. (65), the transmission coefficient is obtained as

Hence, in general, for an -layered structure shown in Fig. 5, the transmission coefficients can be written as follows

Figure 5

()-layered 3D graphene-loaded multilayered structure.

Note that the scattering matrices and are square matrices of size . The entry in the th row and the th column, i.e., represents the reflection/transmission coefficient of th mode when th mode is incident on the input port. For example, the entry denotes that the excitation mode in the input port is the first mode, i.e. the fundamental mode in a circular waveguide, and the reflected/ transmitted mode is the same mode . The entry denotes that the excitation mode in the input port is the first mode and the reflected/transmitted mode is the second mode of the circular waveguide, i.e. mode. In this paper the reflection and transmission coefficients for the entry , i.e. both incident and reflected modes being the fundamental mode, are plotted which pertains to the mode in the circular waveguide and the TEM mode in the coaxial waveguide.

Validation of the proposed method

To validate the proposed method, two multilayer graphene-loaded structures are examined in the cylindrical coordinates system: a circular waveguide and a coaxial cable. Both structures are analyzed for different values of graphene chemical potential and varying numbers of graphene plates. The results of the proposed method and CST simulations are compared to each other.

Analysis of a graphene-loaded circular waveguide

In this section, as an example of a graphene-loaded multilayer structure, graphene plates are placed inside a circular waveguide. The cross-section of the waveguide is shown in Fig. 6a. The radius of the waveguide is mm which is filled with a dielectric permittivity of . Figure 6b depicts the circular waveguide loaded with a graphene plate. For THz and sub-THz regimes in the Kubo formalism of graphene the term pertaining to the intraband electron tramsitions dominates and thus the surface conductivity reduces to a Drude model[14]. In this model the relaxation time , and magnetic bias intensity are assumed as graphene parameters. Note that the conductivity of an unmagnetized graphene plate will become a diagonal tensor with equal diagonal entries[7].

Figure 6

(a) Cross section of the circular waveguide, (b) The waveguide loaded with a graphene plate, (c) Magnitude of transmission and reflection coefficients of the fundamental mode for different values of the chemical potential (red) = 0.05 eV, (blue) = 0.3 eV and (Green) = 2 eV using E-MOL analysis and CST simulations.

Figure 6c depicts the reflection coefficient and transmission coefficient of the fundamental mode of the circular waveguide, i.e. TE11, for the graphene chemical potentials , , and eV. The values and denote the electric field ratios at the ports of the waveguide. It is evident that , , and represent the reflectance, transmittance, and absorptance, respectively, which are power ratios and can be easily calculated from Fig. 6c. The reflectance is due to the discontinuity imposed by the graphene boundary condition and the absorptance is due to the lossy nature of the graphene.

Figure 7a depicts two graphene plates with a spacing inside a circular waveguide while Fig. 7b shows the transmission and reflection coefficients of the fundamental mode for the chemical potentials , and eV. Comparing Figs. 6c and 7b denotes that adding the second graphene plate decreases the transmission coefficient and introduces frequency dispersion in the reflection coefficient. Finally Fig. 8 shows the transmission and reflection coefficients of the circular waveguide loaded with four graphene plates with the spacings for the chemical potentials , and eV.

Figure 7

(a) A circular waveguide loaded with two graphene plates, (b) Magnitude of transmission and reflection coefficients of the fundamental mode for different values of the chemical potential (red) = 0.05 eV, (blue) = 0.3 eV and (green) = 2 eV using E-MOL analysis and CST simulations.

Figure 8

Magnitude of transmission and reflection coefficients of the fundamental mode in a circular waveguide loaded with four graphene plates for different values of the chemical potential (red) = 0.05, (blue) = 0.3 and (green) = 2 eV using E-MOL analysis and CST simulations.

Table 1 presents the time required for the computation of the reflection and transmission coefficients using the E-MoL (implemented by MATLAB) and the CST full wave simulation software. Both methods are implemented under the same conditions exploiting an Intel Core i7 8700K CPU, 32 GB RAM desktop computer. In the data provided by Table 1 the chemical potential is assumed to be eV. As expected, the semi-analytical E-MoL technique is much more efficient than the full wave CST simulation software.

Table 1

The simulation times (in seconds) for E-MOL and CST for the chemical potential 0.3 eV.

Number of graphene plates	E-MOL	CST
One (Fig. 6c)	0.29	562
Two (Fig. 7b)	0.33	589
Four (Fig. 8)	0.34	616

Analysis of a graphene-loaded coaxial cable

In this section the graphene plates are placed inside the coaxial cable and the structure is analyzed in terms of the chemical potential and the number of graphene plates. The cross-section of the coaxial cable is shown in Fig. 9a. The outer and inner radi are and mm, respectively, and the cable is filled with the permittivity = 60. Figure 9b shows the coaxial cable loaded with a single annular graphene plate. Figure 9c shows the transmission and reflection coefficients of the fundamental mode, i.e. TEM mode, for different values of the graphene chemical potential , and eV.

Figure 9

(a) Cross section of the coaxial cable structure, (b) a graphene plate inside the coaxial cable, (c) magnitude of transmission and reflection coefficients of the fundamental mode for different values of the chemical potential (red) = 0.05, (blue) = 0.3 and (Green) = 2 eV using E-MOL analysis and CST simulations.

Figure 10a shows the coaxial cable loaded with two graphene plates with a spacing of = 1 mm. Figure 10b depicts the transmission and reflection coefficients of the cable for different values of the graphene chemical potential , , and eV. Finally, Fig. 11 shows the transmission and reflection coefficients of the cable loaded with four graphene plates with spacings = 1 mm for different values of the graphene chemical potential , , and eV.

Figure 10

(a) The coaxial cable loaded with two graphene plates, (b) magnitude of transmission and reflection coefficients of the fundamental mode for different values of the chemical potential (red) = 0.05, (blue) = 0.3, and (green) = 2 eV using E-MOL analysis and CST simulations.

Figure 11

Magnitude of transmission and reflection coefficients of the coaxial cable loaded with four graphene plates for the fundamental mode for different values of the chemical potential (red) = 0.05, (blue) = 0.3 and (green) = 2 eV using E-MOL analysis and CST simulations.

Considering the reflection and transmission coefficient plots, increasing the chemical potential increases the graphene conductivity and thus reduces its transparency which results in an increase/decrease in the reflection/transmission coefficient, respectively. In other words, the reflected and transmitted powers of the waveguide can be tuned by an electric voltage controlling the graphene plate. This finds applications in tunable microwave waveguide attenuators. Also, this can be utilized in phased array antennas to control the feeding amplitude of each element of the antenna array. Note that no resonant behavior is observed in the reflection and transmission coefficients because the resonance frequency of a graphene plate is far beyond the microwave frequencies. This nearly flat response offers wideband microwave devices.

Conclusion

In this paper, an extended method of lines (E-MoL) is presented for the analysis of graphene-loaded multilayer structures in the cylindrical coordinates. To do so, first, the impedance and admittance transformations imposed by graphene plates are obtained in the cylindrical coordinates system. Having the impedance and admittance transformations, the impedance and admittance matrices can be calculated in each plane of the structure. Then from the input impedance of the structure the reflection coefficient is derived. Finally by transferring the fields from the input to the output of the multilayered stack, the transmission coefficient is computed. To validate the proposed method, a circular waveguide and a coaxial cable loaded with graphene plates are analyzed and the results are compared with the CST full wave simulations. The proposed theory has the potential to model multilayered absorbers and radiating structures accurately and efficiently which is a topic of the future research. The MoL can be generalized to model complicated boundary conditions other than graphene in planar and quasi-planar waveguide structures in integrated microwave and optical circuits and also in circular and conformal antennas. The graphene-loaded optical fibers are attractive optical devices in which the dispersive properties of graphene emerges and can be characterized using the E-MoL preseted in this paper.