Highly improved convergence approach incorporating edge conditions for scattering analysis of graphene gratings

This research developed an effective and efficient approach for improving the slow convergence in the scattering analysis of a one-dimensional graphene grating, made of a periodic array of parallel graphene strips, illuminated by a TM-polarized plane wave. Specifically, the electric fields over the graphene strips and slit regions in a unit cell are individually expressed as an expansion of local basis functions inherently satisfying edge conditions. Interestingly, convergence rate is highly improved compared to the customary and modified Fourier modal method. Additionally, with the aid of local basis functions, the Gibbs phenomenon occurring at both edges of graphene strip can be removed.

The plane wave scattering by a one-dimensional (1D) grating made up of dielectric or metallic mediums has been intensively and extensively studied. Some numerical methods such as the rigorous coupled-wave analysis (RCWA) 1 , the modal theory for dielectric and finitely conducting gratings 2 , and the modal transmission-line method [3][4][5] were developed to accurately calculate scattering characteristics of gratings. Moreover, graphene-based grating composed of graphene sheet has been a continued research interest [6][7][8][9][10][11][12][13][14] in both theoretical studies and practical applications in recent years.
In the RCWA method, both the permittivity function of a periodic medium and electromagnetic fields are expanded into Fourier series and Floquet-Fourier series, respectively. Therefore, the electromagnetic boundaryvalue problem can be converted into an eigenvalue problem. Such an approach is efficient in handling the grating with an arbitrary profile and a finite stack of multiple gratings, as well. However, the RCWA method 1 is known to be slowly converging for 1D metallic gratings in TM polarization (magnetic-field vector parallel to the grating vector). Fortunately, the inverse rule, by invoking adequate Fourier series of the permittivity and reciprocal permittivity functions of a periodic medium to reformulate the eigenvalue problem, was developed to achieved a highly improved convergence rate for the scattering analysis of a metallic grating in TM polarization [15][16][17][18] . Additionally, the Floquet modes in a periodic medium with the unit cell composed of a dielectric slab and a metal layer having finite conductivity can be determined by solving the dispersion equation 19 . However, finding their complex roots is a difficult task. Consequently, use of Fourier series expansion to calculate Floquet modes in a periodic medium, in general, is the most reliable and effective approach in handling a diffraction grating problem.
In this paper, we aimed at studying the numerical convergence of plane wave scattering by a graphene grating in TM polarization. Here, the structure under study is a 1D periodic array of graphene strips (ribbons) deposited on a dielectric substrate. The graphene sheet is assumed to be near-zero thickness (the thickness of mono-layer graphene is 0.335nm); therefore, the electromagnetic fields are merely in the upper and lower homogeneous mediums, as shown in Fig. 1. Due to the periodicity along the x-axis, electric and magnetic fields can be presented in the standard form of Floquet-Fourier series (or Rayleigh expansions). Moreover, the electrical property of graphene strips can be modeled with a surface conductivity ( σ g ). Consequently, the graphene conductivity function in a unit cell is σ (x) = σ g on graphene strip and σ (x) = 0 otherwise, which can be further expressed as a Fourier series. Furthermore, two electromagnetic boundary conditions including (1) the continuous of tangential electric fields across the graphene grating, and (2) the discontinuity of tangential magnetic fields across the graphene grating caused by the conduction current induced on graphene strips should be applied at the interface between two adjacent homogeneous mediums. Alternatively, such a problem amounts to imposing a Scientific RepoRtS | (2020) 10:12855 | https://doi.org/10.1038/s41598-020-69827-w www.nature.com/scientificreports/ periodic boundary condition on the tangential electric-and magnetic-fields at an interface between two uniform mediums. Moreover, the periodic boundary condition is obtained by expanding the graphene conductivity into a Fourier series expansion. The Laurent's rule 16 then can be applied for Fourier factorization of the conduction current σ (x)E x (x, z = 0) . By matching the Fourier coefficient corresponding to the same harmonic, an infinite set of linear equations for the input-output relation of the diffraction-order amplitudes are determined 5,7,20 . Notably, such an approach is termed as Fourier modal method (FMM) throughout this paper. Unfortunately, the poor convergence occurs in conventional FMM for the scattering analysis of periodic arrays of graphene ribbons reported in the literature 13,20 . Furthermore, the inverse rule by including the reciprocal function 1/σ (x) into the FMM is not applicable because σ (x) is zero outside the graphene ribbon. To resolve this problem, the author 20 proposed an approximate boundary condition (ABC) that takes into account the effective conductivity due to the displacement current in slit region without graphene. The effective conductivity is non-zero everywhere; therefore, the inverse rule can be successfully applied. Although FMM with ABC convincingly achieves numerical convergence, the convergent value varies in accordance with the enclosed-loop height (h), shown in Fig. 2, used to model the effective conductivity in slit region. Unfortunately, it is difficult to give a general criterion for h.
In this research, the respective local basis functions (LBFs) taking into account the electric-field edge conditions over the graphene-strip and slit regions in a unit cell are developed to replace the global basis functions (Floquet-Fourier series in a homogeneous medium). As will become clear later on, the tangential electric-field expanded by the local basis functions exhibits the fidelity of discontinuous behaviour, which is not seen in the conventional FMM, enabling a fast convergence in the scattering analysis of graphen-strip gratings in TM polarization.
This paper is organized as follows. We begin with the mathematical formulation using the conventional FMM that is taken as a general framework for theoretical analysis. Moreover, ABC is employed to reformulate the FMM by invoking the inverse rule. Additionally, the present approach, namely, FMM incorporating LBFs, will be comprehensively elaborated. Finally, convergence behaviour will be examined for the three aforementioned methods. Specifically, electric field distribution on graphene grating surface will be demonstrated for various incident conditions. Figure 1 shows structure configuration of a 1D graphene grating consisting of a periodic array of parallel graphene strips. The strip is infinite in extent along the y-axis and the electromagnetic fields have no variation along that direction. The strip width and period along the x-axis are denoted as w g and d x , respectively. The graphene layer is assumed to be zero thickness and characterized by a surface conductivity. Additionally, the 1D graphene strip array is sandwiched by two semi-infinite homogeneous mediums designed as ε + for z ≥ 0 and ε − for z < 0 , respectively. A plane wave with magnetic field H(x, z) =ŷH y is impinging on the array; its incident angle is designated as θ.

Structure configuration and surface conductivity model of graphene
Incidentally, graphene conductivity ( σ g = σ intra + σ inter ), having a close-form expression for the condition | µ c | ≫ k B T , consists of both the intraband ( σ intra ) and interband ( σ inter ) terms 11 :  where -e is the electron charge, ℏ is the reduced Planck constant, γ is a phenomenological carrier scattering rate, µ c is the chemical potential, k B is Boltzmann's constant, and T is the ambient temperature.

fourier modal method
In a homogeneous medium, the electric and magnetic field components of TM polarization in the presence of periodicity d x along the x-axis can be represented by Rayleigh expansions (or Floquet-Fourier series) in the form with ψ n (x) is named as the n th space harmonic given below with the propagation constant along the z-axis k zn = k 2 o ε s − k 2 xn and wave admittance Y n = ωε o ε s /k zn in an infinite medium with relative dielectric constant ε s . Parameters a n and b n usually are termed as the amplitudes of the forward-and backward-propagating waves of the n th diffraction order, respectively.
Owing to the zero thickness approximation of the graphene grating, discontinuity in the tangential component of magnetic fields at the interface (graphene grating surface) between two uniform mediums equals to the conduction current induced on the graphene strips array, yielding where graphene conductivity σ (x) is a periodic function that can be expressed as a Fourier series We first substitute magnetic field H y in Eq. (3) and conductivity function σ (x) in Eq. (10) into Eq. (9) together with the electric field E x on the graphene grating surface approximated by the Floquet-Fourier expansion in Eq. (4). Using Laurent's rule 16 and equaling the same Fourier coefficient corresponding to the same harmonic on both sides, one obtains the system of linear equations where index m is running from negative to positive infinity.
Equation (11) can be expressed in a compact matrix form, one obtains is the Toeplitz matrix with (m, n) entry σ m−n given in Eq. (10); the n th element in column vectors I(0 ± ) and V (0) respectively are I n (0 ± ) and V n (0) .
Eq. (12) establishes a relationship between voltage and current waves across the grating layer; it is the so-called admittance matrix in microwave engineering 5 . Equation (12) is also regarded as the (7) V n (z) = a n exp(−jk zn z) + b n exp(+jk zn z), (8) I n (z) =Y n a n exp(−jk zn z) − b n exp(+jk zn z) , (2020) 10:12855 | https://doi.org/10.1038/s41598-020-69827-w www.nature.com/scientificreports/ input-output relation among incident, reflected and transmitted diffraction-order amplitudes with respect to a graphene grating. Moreover, due to continuous of E x at z = 0 , Eq. (4) arrives at V n (0 − ) = V n (0 + ) ; its vector form can be written as follow: Furthermore, substitution of I(0 + ) = [[Y + ]]V (0 + ) (no downward waves in the upper medium) and Eq. (13) into Eq. (12), one obtains where the input admittance matrix looking into the interface in lower medium is defined as Here the incident and reflected wave vectors, defined at z = 0 − in the lower medium, are individually denoted as column vectors a and b whose elements are a n and b n , respectively. Since the voltage and current wave vectors at z = 0 − can be written as with their components given in Eqs. (7) and (8) Notably, for a graphene grating incident by a single plane wave, the wave vector a is known and usually defined as a 0 = 1 and a n = 0 for n = 0 . The amplitude of each reflected and transmitted diffraction order can then be obtained via Eqs. (15) and (18).
Moreover, the time average Poynting power along the z-axis over a grating period is defined as Therefore, the incident-, reflected-and transmitted-power are obtained as follows: Re n Y − n |b n | 2 , and P tx = 1 2 Re n Y + n |c n | 2 , respectively. The absorptance can be obtained by evaluating e abs = 1 − P ref /P inc − P tx /P inc .
Notably, the mathematical formulation in FMM is rigorous and the result is exact when the space harmonic index n runs from − ∞ to + ∞ ; however, they have to be truncate into [−N, +N] in numerical computation, where N is denoted as truncated order. The total number of space harmonics is then designated as N tot equal to 2N + 1.

fourier modal method with approximate boundary condition
As reported in the literature 15,16 , the slow convergence in the RCWA for TM polarization is not caused by the Fourier series expansion but the form where the Fourier series of the permittivity and the reciprocal permittivity functions are utilized. The same problem also occurs in graphene gratings with thickness of almost zero. Nevertheless, the inverse rule 15,16 is not applicable in the graphene grating because 1/σ (x) goes to infinity in the slit region. The approximate boundary condition on periodic arrays of graphene ribbons was proposed to replace σ (x) by the effective conductivity function: σ eff (x) incorporating the contribution of displacement current 20 in slit region. In doing so, the reciprocal function 1/σ eff (x) exists everywhere. More precisely, by Ampere's law with Maxwell's modification, the discontinuity in magnetic fields intensity (H) on the graphene grating between two uniform mediums equals to the sum of conduction-and displacement-current. At the left-hand side of Fig. 2, the line integral over the graphene strip equals to the conduction current flowing along the x-axis; while it equals to the displacement current filled in the rectangular box at the right-hand side figure. Combining these two results, Eq. (9) can be rewritten as 20 where the effective conductivity at z = 0 is www.nature.com/scientificreports/ Parameter σ (x) is the graphen conductivity function given in Eq. (10) and ε = (ε + + ε − )/2 . Parameter h is the height of the rectangular enclosed loop in Fig. 2. Now, the reciprocal function 1/σ eff (x) exists everywhere as well as its Fourier expansion. By the standard procedure of FMM in the previous section incorporating inverse rule 20 , we obtain the same input-output relation in Eq. (12)  Once the admittance matrix is obtained, the standard procedure of FMM in the previous section can be applied to calculate the reflection, transmission and absorption efficiencies.

fourier modal method incorporating local basis functions
Although FMM with ABC can improve the convergence rate, the parameter h is still a core factor affecting the value of convergence, as will be demonstrated in the next section. Moreover, the criterion for determining the appropriate h remains to be studied in detail. So far, we have implemented a computer program based on FMM with ABC for calculating the scattering characteristics of a graphene grating as well as the electric-field ( E x ) distribution on the graphene grating surface. Let us skip ahead to the numerical result concerning a graphenestrip grating normally incident by a TM-polarized plane wave. Figure 5 shows the distribution of E x (red dotted curve obtained using FMM with ABC) over the graphene-strip ( x ∈ [0, 20] ) and slit ( x ∈ [20, 70] ) regions. It is obvious to see jump discontinuities occurring around x = 0 and x = 20µm ; wherein E x (or J x /σ g ) vanishes close to strip edges. In fact, induced current vanishes at strip edges can be explained physically as follows.
Referring to Fig. 1, at x = w + g in the slit region, magnetic field component H y must be continuous across the interface due to J Additionally, the vanishing current at metal-strip edges was reported in the research of a metal-strip grating illuminated by a plane wave in TM polarization 19 , but unfortunately the distribution of E x in the slit region was not shown in that paper. Regarding the electric field in the slit region, an exponential growth of E x around the slit edges can be observed (red dotted curve) in Fig. 5. Alternatively, as is well know, TM-polarized electric field near the edge of a thin sheet is proportional to ρ −1/2 and becomes singular as ρ approaches zero 21 , where ρ is defined as the radius (in the polar coordinate system) with its original point locating at the edge.
Nevertheless, the Gibbs phenomenon taking place near the discontinuities reveals that the customary global basis of Floquet-Fourier series is inappropriate for expanding the field directly on the graphene grating surface. In view of that, it is essential to construct the local basis functions inherently satisfying the field nature in respective regions. Furthermore, the criterion for choosing basis functions contains: (1) to use only a few basis functions to approach the correct solution, and (2) to have closed forms in the overlap integral between the local basis functions and the space harmonic in Eq. (5).
More specifically, the sin-based local basis function vanishing at its both ends for any harmonic order is used to expand E x over the graphene strip, which is given as ; index n is ranging from 1 to N g . On the other hand, in the slit region, we have the singular basis functions with singularities at its two edges, which are commonly used to approximate the current parallel to the edges in a micro-strip line 22  In a unit cell on the graphene grating surface (at z = 0 ), E x (x) can be written as Scientific RepoRtS | (2020) 10:12855 | https://doi.org/10.1038/s41598-020-69827-w www.nature.com/scientificreports/ Parameters N g and N s represent the number of basis in the graphene strip and slit regions, respectively. Due to electromagnetic boundary condition, E x must be continuous at the interface between graphene grating and uniform medium at z = 0 . Therefore, equality of Eqs. (4) and (25) gives By multiplying the complex conjugate of ψ m (x) on both sides of Eq. (26) and taking integration over one period, we obtain where integer m is ranging from -N to +N. The notation of �a(x)|b(x)� =    Here, we obtain a totally different admittance matrix while the input-output relation in Eq. (12) remains the same. The same procedure in FMM can be applied to calculate the reflect and transmit amplitudes of each diffraction order. Once the voltage V (0) = [[T]]a is obtained, the expansion coefficients of local basis functions, p and q in Eq. (31) can be readily determined, as well as the distribution of E x on the graphene strip and slit.
Furthermore, the power dissipated on the graphene strips array can be directly determined by

numerical results and discussion
A free-standing graphene-strip grating is taken as an example to examine the convergence behaviour for the three approaches. The parameters of graphene are µ c = 0.39 eV , T = 300 K , and ℏγ = 0.658 meV (relaxation time of charge carriers τ = 1/2γ = 0.5 ps ). The period and strip width are d x = 70 µm and w g = 20 µm , respectively. The upper and lower semi-infinite mediums are free space with ε + = ε − = ε o . We first calculate the absorptance against the truncated order N running from 1 to 200. In Fig. 3, the three methods, including the conventional FMM, FMM with ABC, and the present approach FMM incorporating LBFs, were employed to carry out the convergence test. Three different enclosed-loop heights (h) in FMM with ABC are considered: w g /200 (magenta dashed curve), w g /2000 (blue dashed curve) and w g /20000 (red dashed curve). As was widely reported in literature [15][16][17][18]20 , the conventional FMM (green dotted curve) indeed runs into the serious problem of oscillating convergence behaviour. On the other hand, the FMM with ABC can improve the convergence rate; however, the convergent value changes with h accordingly. Interestingly, the result with h = w g /2000 approaches to that of our method. Apparently, the convergence rate of our approach (black dashed curve) is superior to the other two methods; even only a few number of truncated order is needed to achieve the numerical convergence.
Additionally, the absorption versus frequency for both convincing methods including the FMM with ABC and our approach are demonstrated in Fig. 4. Because of oscillating convergence in the conventional FMM, its result is unreliable and was neglected here. The aforementioned graphene and structure parameters are used in this example. The graphene grating is obliquely incident by a TM-polarized wave with incident angle 60°. The enclosed-loop height ( h = w g /2000 ) is chosen since it shares approximately the same convergence value with that of our approach in Fig. 3. Two truncated orders (N) are used to examine the performance of numerical convergence. It is obvious to see that FMM with ABC (yellow solid curve) and our approach (purple solid curve) agree very well for the case of N = 100 . Specifically, the result of our approach with N = 35 (red solid curve) coincides with those of N = 100 . However, the result obtained by FMM with ABC for N = 35 (blue solid curve) shows apparent discrepancy, particularly in the higher frequency range. Additionally, the Wood's anomalous taking place at 2.3 THz can be observed in both approaches. We may conclude that compared to FMM with ABC, our approach can achieve numerical convergence even if only a few truncated orders is used. www.nature.com/scientificreports/ Figure 5 depicts the absolute value of E x (x, z = 0) , which was normalized to the incident E x , versus x in a unit cell on the graphene grating surface. The red dotted curve is calculated based on FMM with ABC, while the solid curve in blue colour is obtained using FMM incorporating LBFs. The region for x ∈ [0, 20µm] is the graphene strip, and is otherwise the slit region. Although the Gibbs phenomenon, an overshoot (oscillating) of a Fourier series occurring at jump discontinuities, is obvious in red dotted curve based on FMM with ABC, the vanishing current density ( σ g E x ) at the strip edges and exponentially growth in E x around the slit edges can still be clearly observed. Contrarily, owing to the two local basis functions, given in Eqs. (23) and (24), inherently satisfy the individual edge condition, the Gibbs phenomenon is removed, as shown in the blue solid curve. Although not shown here, the case of N = 35 shares almost the same profile with the case of N = 100 ; this means that only a few LBFs is needed to expand E x . Notably, in this case the FMM with ABC can achieve almost the same result of absorptance for the case of N = 100 ; however, it can not reflect the essence of field nature around the graphene strip edges. Additionally, the operating frequency 2.5 THz is near the absorption peak; the incident wave is resonant with the graphene current J x along the x-direction. Therefore, the induced current exhibits the first normal mode (standing wave) pattern. It is very similar to the current induced on a radio-frequency (RF) dipole antenna excited at its first resonant frequency.
Since the approach of FMM incorporating LBFs can effectively remove the Gibbs phenomenon, it should be in a good position to observe the effect of incident angles on the electric field ( E x ) distribution over the  Comparison of |E x (x, z = 0)| , normalized to the incident E x in a unit cell on the graphene grating surface using the FMM with ABC (dotted curve in red colour) and our approach incorporating LBFs in FMM (solid curve in blue colour). The operating frequency and incident angle are 2.5 THz and θ = 0 • , respectively. The truncated order is N = 100 . The region for x ∈ [0, 20 µm] is in the graphene strip; otherwise is in the slit region.
Scientific RepoRtS | (2020) 10:12855 | https://doi.org/10.1038/s41598-020-69827-w www.nature.com/scientificreports/ graphene-strip grating surface. In Fig. 6a and b, the absolute value of E x (x, z = 0) normalized to the incident E x on the graphene grating surface against x-axis was demonstrated for the five cases with different incident angles including 0°, 15°, 30°, 45°, and 60°, respectively. Fig. 6a shows the electric field strength over the graphene strip, while Fig. 6b shows that in the slit region at z = 0 . As shown in Fig. 6a, the induced E x (or J x /σ g ) on the graphene strip changes insignificantly as θ is increasing from 0° to 30°, while its peak is decreasing as the incident angle increases up to 45° and 60°. Interestingly to find that the current distribution is almost symmetric with respect to the graphene strip centered at x = 10 µm even for the oblique incidence. As was reported in literature 13 , the resonance effects are in connection with the leaky plasmonic modes existing in individual graphene strip but with weak coupling between strips. Such a normal mode is a source-free solution; therefore, its mode pattern is almost independent of the incident angle. Contrarily, the asymmetric distribution is observed in the slit region for oblique incidence, which maybe caused by unbound waves with continuous spectrum.

conclusion
In this research, the three approaches including the conventional FMM, FMM with ABC, and our approach incorporating LBFs in FMM were implemented to examine the convergence behaviour of absorptance with respect to a periodic array of parallel graphene strips obliquely incident by a TM-polarized plane wave. Because of the individual local basis functions inherently satisfying the electric-field edge condition at graphene-strip and slit edges, the Gibbs phenomenon due to the Fourier expansion of global basis functions (space harmonic) in conventional FMM disappears. Furthermore, the convergence rate of the present approach is superior to the other two methods. Additionally, a new admittance matrix is obtained in the present approach while the whole formulation can still fit into the standard procedure of FMM. Significantly, the inverse rule and ABC for FMM are no longer needed. Such an approach can drastically reduce the required number of space harmonics and is more efficient for scattering analysis of stacked multiple graphene gratings.