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.

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][21][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 k 0 according to z = k 0 z and r = k 0 r . 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 H z and E z are obtained as follows: 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 D •.• r and D •.• φ , 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 use the following abbreviations: Figure 1. The discretization scheme in the cross section of a cylindrical waveguide with inhomogeneous layers. in which In Eqs. (15) and (16), M 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 where F = Eor H . 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: www.nature.com/scientificreports/ 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 k and k + 1 . This common interface is denoted by the letter B. Besides, the two sides of this interface are distinguished by + and − indices. As shown in Fig. 2, the graphene is placed in xy or equivalently rϕ plane in which the GTL equations are discretized. Due to the semi-analytical nature of MoL, the equations are solved analytically in z 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. Y  www.nature.com/scientificreports/ 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 E r and E ϕ are discretized in different positions and hence we cannot add or subtract them directly to obtain H r . 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 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. Z 0 = I . At the input of the structure there is an input waveguide with a wave impedance Z in 0 or a wave admittance Y in 0 . Using the admittance transformation formulas, the input admittance Y A 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 S 11 we have If the modes are ordered properly, the first element in the S 11 matrix gives the fandamental mode reflection coefficient. www.nature.com/scientificreports/ 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 f and b 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 Z (k) 0 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 as Assume 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: (1) Using Eqs. (50)  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.
The values of the admittances Y A + .Y A − and impedances Z A + .Z A − shown in Fig. 3, are determined by the aforementioned admittance and impedance transformation formula Eq. (43). Also, the values Y B and Z B are calculated from Eq. (25). Assuming that these values are known, the transformation coefficient S 21 is derived. Using Eq. (53), the total field at this plane is written as follows www.nature.com/scientificreports/ The electric field at the interface of the plane A is continuous and hence E A + = E A − . Due to normalization, the impedance and admittance values at the plane A − are equal to the characteristic impedance of the waveguide Z A − = Y A − = I . From Eq. (48), the values of the forward E fA − and backward E bA − fields at the plane A − are obtained as Therefore, we have Considering that the applied field at the port B is equal to E fB , then the forward propagating field (towards z direction) at the plane A + will become E fA + = e −Ŵd AB E fB , hence Then, the transmission coefficient matrix is calculated as The first element of the matrix S 21 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 E fC , then the forward electric field at the plane B + is calculated as E fB + = e −Ŵd BC E fC . Considering the continuity of the tangential electric field at the plane B, we have Therefore, the forward and backward fields at the plane B − are obtained as The forward field at the plane A + is written as Consequently, using Eq. (53), the total electric field at the plane A + can be written as Considering the continuity of the electric field components at the plane A we have Since the electric field at the plane A − is completely propagating forward, from Eqs. (61)-(63) we can write Hence, in general, for an (N + 1)-layered structure shown in Fig. 5, the transmission coefficients can be written as follows Note that the scattering matrices S 11 and S 21 are square matrices of size 2N ϕ N r × 2N ϕ N r . The entry in the i th row and the j th column, i.e.S ij 11 /S ij 21 , represents the reflection/transmission coefficient of i th mode when j th mode is incident on the input port. For example, the entry (1, 1) denotes that the excitation mode in the input port is the first mode, i.e. the fundamental TE 11 mode in a circular waveguide, and the reflected/ transmitted mode is the same mode (TE 11 ) . The entry (2, 1) denotes that the excitation mode in the input port is the first mode (TE 11 ) and the reflected/transmitted mode is the second mode of the circular waveguide, i.e. TM 01 mode. In this paper the reflection and transmission coefficients for the entry (1, 1) , i.e. both incident and reflected modes being the fundamental mode, are plotted which pertains to the TE 11 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 grapheneloaded 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 a = 10 mm which is filled with a dielectric permittivity of ε r = 60 . 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 . . . www.nature.com/scientificreports/ and thus the surface conductivity reduces to a Drude model 14 . In this model the relaxation time τ = 0.1ps , and magnetic bias intensity B 0 = 0 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 6c depicts the reflection coefficient |S 11 | and transmission coefficient |S 21 | of the fundamental mode of the circular waveguide, i.e. TE 11 , for the graphene chemical potentials µ c = 0.05 , µ c = 0.3 , µ c = 1.2 and µ c = 2 eV. The values |S 11 | and |S 21 | denote the electric field ratios at the ports of the waveguide. It is evident that |S 11 | 2 , |S 21 | 2 , and 1 − |S 11 | 2 − |S 21 | 2 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 d = 1mm inside a circular waveguide while Fig. 7b shows the transmission and reflection coefficients of the fundamental mode for the chemical potentials µ c = 0.05 , µ c = 0.3 and µ c = 2 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 d = 1mm for the chemical potentials µ c = 0.05 , µ c = 0.3 and µ c = 2 eV. 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  www.nature.com/scientificreports/ provided by Table 1 the chemical potential is assumed to be µ c = 0.3 eV. As expected, the semi-analytical E-MoL technique is much more efficient than the full wave CST simulation software.
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 a = 10 and b = 2.5 mm, respectively, and the cable is filled with the permittivity ε r = 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 µ c = 0.05 , µ c = 0.3 and µ c = 2 eV. Figure 10a shows the coaxial cable loaded with two graphene plates with a spacing of d = 1 mm. Figure 10b depicts the transmission and reflection coefficients of the cable for different values of the graphene chemical potential µ c = 0.05 , µ c = 0.3 , and µ c = 2 eV. Finally, Fig. 11 shows the transmission and reflection coefficients of the cable loaded with four graphene plates with spacings d = 1 mm for different values of the graphene chemical potential µ c = 0.05 , µ c = 0.3 , and µ c = 2 eV.
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  www.nature.com/scientificreports/ 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 quasiplanar 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.  www.nature.com/scientificreports/

Data availability
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.