Method of lines for analysis of plane wave scattering by periodic arrays of magnetically-biased graphene strips

In this paper, efficient analysis of the plane wave scattering by periodic arrays of magnetically-biased graphene strips (PAMGS) is performed using the semi-numerical, semi-analytical method of lines (MoL). In MoL, all but one independent variable is discretized to reduce a system of partial differential equations to a system of ordinary differential equations. Since the solution in one coordinate direction is obtained analytically, this method is time effective with a fast convergence rate. In the case of a multi-layered PAMGS, the governing equations of the problem are discretized concerning periodic boundary conditions (PBCs) in the transverse direction. The reflection coefficient transformation approach is then used to obtain an analytical solution in the longitudinal direction. Here, magnetically-biased graphene strips are modeled as conductive strips with a tensor surface conductivity which is electromagnetically characterized with tensor graphene boundary condition (TGBC). The reflectance and transmittance of different multi-layered PAMGS are carefully obtained and compared with those of other methods reported in the literature. Very good accordance between the results is observed which confirms the accuracy and efficiency of the proposed method.

reported for the diffraction analysis of GPSs in the absence of magnetic field bias while less effort has been made in case of magnetically-biased GPSs 20,27 .
In this paper, the semi-numerical, semi-analytical MoL has been provided for the efficient analysis of the plane wave scattering by PAMGS. In this method, we use of PBCs to take the periodicity of the PAMGS into account. Hence, the fields and the derivatives are to be discretized with finite differences according to the PBCs. Besides, the constitutive parameters of the cross section as well as the surface conductivity of graphene strips are discretized. In those parts of the cross-section where there is no graphene, the conductivity value is zero. Hence, there is no need to introduce approximate boundary condition at the interface of two adjacent layers where graphene exists as done by Khoozani et. al 27 .
Moreover, unlike fully numerical methods, the solution of the problem in one coordinate direction is obtained analytically in MoL which means there is no need to use absorbing boundary condition (ABC) such as perfectly matched layers (PML) in this direction. Besides, the analytical nature of the method saves a lot of computing time compared to fully numerical methods. The reflection coefficient transformation approach is introduced to obtain an analytical solution in this direction. In this approach, the structure under the study is divided into homogeneous layers and the known reflection coefficient in the last layer is transformed through layers and interfaces to obtain the reflectance. Having known the reflection coefficient in each layer, the transmission coefficient is obtained by transforming the fields in the opposite direction.
The rest of the paper is organized as follows: First, we present a brief description of the graphene properties. Then, we continue with the basic formulation of the problem. In this section, the MoL is introduced concisely. Afterwards, the plane wave incidence on a PAMGS for the two special cases of normal and oblique incidence is investigated. The reflection and transmission coefficients of a multi-layered PAMGS are obtained in the next subsection. Next, the numerical results of different PAMGSs obtained using the proposed method are presented.

Graphene conductivity model
Consider a graphene sheet in the xy plane with a magnetic field bias B = B 0ẑ as shown in Fig. 1a. As the result of its gapless electronic band structure and one-atom thickness, graphene can be treated as semiconductors 28,29 characterizing with an anisotropic surface conductivity to be obtained through a quantum mechanical analysis using Kubo formulation 8,30,31 . However, in the case of magnetic field bias, there exists an excitonic energy gap which is opened due to the interaction of electrons with the magnetic field that can be significant at low temperature. At high temperature, the carrier density increases and consequently the screening effect enhances. Hence, the gap disappears and the graphene sheet's electronic band structure remains gapless 32 .
The anisotropic nature of graphene's conductivity can be interpreted by the electron motion in the xy plane under an electric field as shown in Fig. 1b 8 . For the case E = E yŷ , the electric force F e will be in the −ŷ direction. The acceleration of the electron in −ŷ direction causes a Lorentz force F m in x direction. Consequently, the combination of these two forces results in two current components in −x and ŷ directions. The case of E = E xx can be treated similarly in which the two current components are produced in x and ŷ directions. Thus, the relation between the resultant surface currents and tangential electric fields is described as follows 8 in which σ xx = σ yy = σ L and σ xy = −σ yx = σ T are the longitudinal (parallel to the electric field) and transversal (perpendicular to the electric field) components consisting of two parts modeling the intraband and interband electron transitions. However, in THz regime where ω < 2µ c , the intraband transition is dominant 8,27 . Here, the optical properties of graphene are similar to those of Drude-type materials and as the result the quantum mechanical closed form expressions for σ L and σ T reduce to the classic Drude model as follows 8,27,33 with (2) www.nature.com/scientificreports/ where is the reduced Plancks constant, e is the electron charge, µ c is the chemical potential, τ is relaxation time, ω c = eB 0 v 2 F /µ c is the cyclotron frequency, v F = 10 6 m/s is the Fermi velocity, n s is the carrier density, and T is the Kelvin temperature.
Alternatively, since the graphene sheet is modeled with a tensor surface conductivity, it can be electromagnetically treated as an special impedance boundary condition at the interface of adjacent layers called tensor graphene boundary condition (TGBC). The TGBCs are given by in which n is the unit vector normal to the interface pointing from media 2 to media 1, ( E 1 , H 1 ) and ( E 2 , H 2 ) are the electric and magnetic fields in media 1 and 2, respectively.
These boundary conditions are then used in the next subsections to derive the reflection coefficient transformation formula through interfaces. Note that a time dependence of e jωt is assumed throughout this paper.

Formulation of the problem
Consider a PAMGS with graphene strips of width W g and period length of D as shown in Fig. 2a in which the graphene strips are assumed to be infinitely long in the y-direction.
Despite graphene sheets with gapless electronic band structures, one-dimensional graphene strips with finite width have been proposed to obtain tunable band gaps where the energy gap is controlled by the strip width 28,29 . Besides, unlike infinite graphene sheets, graphene strips have non-uniform carrier density n s which means that the electrons near the edges of the strips experience different forces in compare to those in the middle of the strips. As the result, two different modes can be supported by the graphene strips namely edge modes and bulk modes with high localization near the edges and in the middle of the strips, respectively 33 . For chemically doped graphene, the net charge is zero. Here, the edge effects are confined to a strip width in order of the dopant's Wigner-Seitz radius. Hence, for wide enough strips these edge effects can be neglected and the graphene strips can be modeled with the same conductivity parameters as those of an infinite sheet 33 .
The method of lines. To start the analysis of a PAMGS and to consider its periodicity in x-direction, we normalize the electromagnetic fields with respect to the phase as follows: in which ψ 0i (x, y, z) is the non-normalized electromagnetic field component, k x is the propagation constant of the incident wave vector in x-direction and ψ i (x, y, z) is either E i or H i with i = x, y, z . These functions are periodic in the x-direction with period D.
Next, we eliminate the longitudinal field components E z and H z in Eqs. (6) and (7) to obtain generalized transmission line (GTL) equations as follows 34 In this step, we discretize the fields and the derivatives in Eqs. (8) and (9) using finite differences. To this end, consider the unit cell of the PAMGS as shown in Fig. 2b with period length of D and N x discretization lines within the period. To start, the periodic functions of Eq. (5) are discretized as follows in which the S e and S h are diagonal matrices of size N x × N x whose elements are e jkx ihx and e jkx (i+0.5)hx for i = 1 : N x , respectively. The discretized quantities are displayed in bold. Accordingly, the discretize field components of E y and H y are saved in column vectors of size N x × 1 . Then, the phased normalized derivatives in the discretized domain are obtained as 34 with in which s = e j kx hx 2 . Hence, D ē x = DxE y and D h x = DxH y are constructed according to Eq. (11). Besides, since there is no variation in the y-direction we put Dȳ = [0] N x ×N x in Eqs. (8) and (9). As the result, the discretized GTL equations are obtained as follows: where I is an identity matrix of N x .
Using the definition E = [E y E x ] t and H = [−H xHy ] t where the superscript t refers to the matrix transpose, the solution to Eqs. (13) and (14) are given by 34 d dz ε r I .
and (H f ,H b ) are forward and backward amplitudes of propagating electric and magnetic fields eigenmodes, respectively.

Plane wave incidence. According to Maxwell's curl equations
, the relation between transverse components of the electric and magnetic fields is given by where the subscript t refers to transverse components. In what follows, we investigate two special cases of normal and oblique incidence on a PAMGS using MoL.
Normal incidence. According to Fig. 3a, for the case of normal incidence we put k x = 0 in Eqs. (13) and (14) and compute the corresponding T E , T H and Ŵ in the two half spaces. Next, we define incident, reflected and transmitted electric and magnetic fields as follows:

By definition
To obtain the unknown reflection and transmission coefficients of the PAMGS, the TGBCs at the interface of layers 1 and 2 must be applied. The discretized form of Eq. (4) reads as For W g = D , the PAMGS shown in Fig. 3a reduces to an infinite graphene sheet at the interface of two infinite layers. In this case the reflection and transmission coefficients can be obtained analytically 8 . Figure 4a compares the obtained Co-and Cross-polarized reflection and transmission coefficients of a free standing infinite graphene sheet using MoL and analytical formulation. Besides, in order to examine the accuracy of the proposed method we calculated the relative error between our proposed method and those of analytical formulations for two  Oblique incidence. For the case of oblique incidence, the two particular cases of TE(H) and TM(E) polarized waves are investigated simultaneously. According to Fig. 3b, we put k z = k 0 cos(θ i ) with i = 1, 2 in Eqs. (18) and (19) which leads to where with Taking Eqs. (29) to (31) into account and following the same procedure as done in the previous case, the reflection and transmission coefficients with very small modification compared to normal incidence are obtained as follows where It is worth mentioning that since the effect of the radiation angle (θ 1 ) on reflectance and transmittance is considered analytically through Eqs. (29) to (31), there is no need to recalculate the eigenvector matrices (T Ei , T Hi ) in case of oblique incidence. The reflection and transmission coefficients of a free standing infinite graphene sheet under oblique incidence with θ 1 = 30 • are obtained using Eqs. (32) and (33) and compared with those of analytical solutions as shown in Fig. 5a. Again, we computed the relative error between our results and those of Multilayer PAMGS. In the analysis of a multilayer PAMGS as shown in Fig. 6, we use the reflection coefficient transformation approach to obtain the reflection and transmission coefficients of the whole periodic structure.
In this approach the output reflection coefficient (r out ) is assumed to be known and is to be transformed through layers and interfaces to obtain the input reflection coefficient (r in ).

By definition
where i is the number of layer. On the other hand, Ē f (z) = e −ŴzĒ f (0) and Ē b (z) = e ŴzĒ b (0) . Hence, according to Eq. (36) the reflection coefficient is transformed through the layers by For i = N we have r (i+1) − = r out which is set to zero in case of a semi-infinite layer. To calculate the reflection coefficient transformation formula through interfaces, the TGBCs have to be applied as follows Now we substitute Eq. (29) into Eq. (38) and rearrange the formulas to obtain r i − as a function of r i + as follows with Using the extracted recursive formulas for i = N : 1 , the reflection coefficient of the multilayer PAMGS as r in = r 1 − is obtained.
To obtain the transmission coefficient of the structure we should transform the fields in the opposite direction from the input to the output as follows By definition Hence, [B] = (I + r i + ), In short, we considered graphene as an especial impedance boundary condition at the interface of two adjacent layers due to its 2D structure and computed the reflection coefficient transformation formula through the interfaces, accordingly. Other 2D materials including Group-IV monoelemental honeycomb materials and binary compounds of group III-V elements can be treated similarly provided that we can model them as an appropriate boundary condition.

Numerical results of different PAMGSs
To validate the proposed method, the numerical results for two different PAMGSs are obtained and compared with those of other methods. As the first example, the Faraday rotation effect achieved by a PAMGS as shown in Fig. 7a is investigated using MoL and compared with Fourier-Floquet plane wave expansion method 36 . Then, a multilayer PAMGS as a tunable THz absorber shown in Fig. 7b is analyzed. In this case, the obtained results are compared with those of the circuit model theory method. In both cases, the two sets of the results are in good accordance. Fig. 7a is illuminated with a x-directed wave at normal incidence.

Faraday rotation. The PAMGS shown in
In presence of a magnetic field bias, the direction of the transmitted electric field is deviated from that of the incident wave by angle θ F called as Faraday rotation angle. By definition 36 where t cross = t yx and t co = t xx are the transmission coefficients in y-and x-direction due to an x-directed incident wave, respectively.
On the other hand, the transmission efficiency is defined by 36 in which n 1 and n 2 are the refractive indices of incident and transmitted regions, respectively. The Faraday rotation angle and the transmission efficiency of the structure shown in Fig. 7a with D = 2W g = 1µm and ε r = 4 as a function of frequency are obtained using MoL and compared with those of Ref. 36 and a commercial software as shown in Fig. 8a,b, respectively. In order to compare our numerical results with those obtained from commercial software, the structure shown in Fig. 7a is simulated using High Frequency Structural Simulator (HFSS) with its modal solution solver. It is worth mentioning that for this structure the number of discretization lines is set to N x = 50 for the MoL results to be converged. Moreover, to match the HFSS results to that of MoL, maximum Delta S is set to 0.0005 and took more than 40 minutes to run on an Intel Core i7 8700K CPU, 32 GB RAM desktop computer while the MoL took less than 4 seconds on the same computer. www.nature.com/scientificreports/ As it is seen the three sets of results agree very well. According to Fig. 8a,b the minimum transmittance occurs where the θ F changes sign. Moreover, to investigate the effect of the chemical potential on reflectance and transmittance, the co-and cross-polarized reflectance and transmittance for different µ c are obtained and shown in Fig. 9a,b, respectively. The resonance frequency of the PAMGS is blue-shifted as µ c increases.
Tunable THz absorber. The next example deals with a tunable THz absorber reported made of a multilayer PAMGS as shown in Fig.7b 6 . This absorber is composed of a perfect electric conductor (PEC), a continuous graphene sheet (CGS), and the periodic array of graphene strips (PAGS). The PEC, CGS and PAGS are separated with an intermediate dielectric layer. The PEC is used to suppress the transmission and increase the absorption.
To evaluate the efficiency and accuracy of the proposed MoL in analyzing this absorber with D = 23.6µm , W g = 18.28µm , d 1 = 300µm , d 2 = 0.3µm and ε r = 2.34 , the absorption of the structure as a function of frequency is obtained and compared with those of circuit model theory approach reported in literature 6 as shown in Fig. 10. The graphene sheet and strips are assumed to be isotropic ( B 0 = 0 ) and the structure is illuminated  www.nature.com/scientificreports/ by an x-directed wave with normal incidence. As illustrated in Fig. 10, the obtained results using MoL agree very well with those of circuit model theory approach. Next, the effect of incident angle on absorption for both TE and TM polarizations is investigated. As shown in Fig. 11a,b, increasing the incident angle causes a blue shift in absorption peaks. Besides, the absorption of the structure for both TE and TM polarizations are almost the same. As the result, this absorber can be referred to as a polarization independent absorber 6 .

Conclusion
In summary, we used the semi-numerical, semi-analytical method of lines (MoL) for the efficient analysis of the plane wave scattering in the most general case of oblique incidence for both TE and TM polarizations by multilayer periodic arrays of magnetically biased graphene strips (PAMGS). First, we discretized the governing equations in transverse directions concerning the periodic boundary condition using finite differences. Then, the reflection coefficient transformation approach was used to obtain the reflectance and transmittance of the multilayer structure analytically. Since the graphene parameters are discretized with finite differences in transverse directions, we set the surface conductivity to zero where graphene is absent. Hence, there is no need to introduce approximate boundary condition to model graphene. On the other hand, due to the semi-analytical nature of the method, there is no demand for the use of absorbing boundary conditions in the longitudinal direction unlike fully numerical methods. These advantages have made MoL a powerful computational tool for analyzing multilayered graphene-based periodic structures. To validate the proposed method, numerical results for different PAMGSs are obtained using the proposed method and compared with those obtained by other methods reported in the literature. In all cases, the two sets of the results agree very well, which confirms the validity of the method. Analyzing 3D graphene-based periodic structures with 2D periodicity will be the subject of our future work.

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