Extraordinary optical transmission through periodic Drude-like graphene sheets using FDTD algorithms and its unconditionally stable approximate Crank–Nicolson implementation

Based upon the approximate Crank–Nicolson (CN) finite-difference time-domain method implementation, the unconditionally stable algorithm is proposed to investigate the wave propagation and transmission through extremely thin graphene layers. More precisely, by incorporating the CN Douglas–Gunn algorithm, the piecewise linear recursive convolution method and the auxiliary differential equation method, the analytical model is proposed for Drude-like graphene model. To obtain the solution of the governing equations, the perfectly matched layer and the periodic boundary condition are applied to the graphene structure with two dimensional nano-materials. Numerical examples are carried out for further investigation. During the simulation, the influences of the parameters such as the grating slit and its thickness on the wave transmission are investigated and discussed. The result shows that not only the graphene grating has high transmission performance but also the proposed methods have considerable performance and accuracy.


Theoretical approach
The Drude dispersion model can be applied to GSs 22 . According to the Drude model, the relative dielectric function ε r (ω) is given by where ω is the angular frequency of the impinging light, ε ∞ is the relative permittivity at infinite frequency, ν is the damping constant in the graphene sheet, and ω p is the plasma frequency given by where n is the density of electrons, m * is the electron effective mass, e is the electron charge, and ε 0 is the permittivity of the vacuum 23,24 . The source-free and normalized frequency-domain modified Maxwell's curl equations in two dimensional Drude media can be written as where S η , η = x, z is the stretched coordinate variables in the higher order PML regions, given as where κ ηn is assumed to be positive integer, σ ηn and α ηn are assumed to be positive, n = 1, 2 . By employing the ADE approach, S −1 η can be rewritten as where the coefficients can be given as a ηn = α ηn ε 0 , b ηn = α ηn ε 0 + σ ηn κ ηn ε 0 , k η = 1 κ η1 · 1 κ η2 . By substituting (7) into (3)- (5), the equations can be rewritten as Scientific Reports | (2020) 10:17462 | https://doi.org/10.1038/s41598-020-74552-5 www.nature.com/scientificreports/ Introducing the auxiliary variables, they can be given as the following forms where η is the rest part of the field components, for example, when calculating H y , η = x while η = z . Thus, the Maxwell's equations in the higher order PML regions can be expressed by the auxiliary variables in the timeharmonic domain as By substituting (11)- (14) into (15)- (17), the equations can be rewritten as Based on the CN scheme and PLRC methods, we obtain the FDTD difference equations as follows 25,26 where the coefficients can be expressed as And the operator δ η represents the CN scheme, for example, The updated equation of auxiliary field obtained by PLRC method can be calculated as By substituting (21) and (22) into (23), the updated equation of magnetic component in the y-direction can be given as where D 2η = p 3η p 6η δ η δ η . By employing the CNDG algorithm, adding D 2x D 2z H n+1 y and D 2x D 2z H n y at both sides of the equations and split the resultants, one obtains where A n is the other terms at right sides of (27). The equations can be solved by two steps as To obtain the solution of Maxwell's equations for the periodic structure, the PBC is introduced. Figure 1 shows the sketch picture of electric and the magnetic fields in a periodic structure. The PBC is used at upper and bottom of the x-direction. The boundary of z-direction in computational area is truncated by the proposed PML scheme. The boundary of the x-direction in the computational area is surrounded by PBCs due to its periodicity It can be observed that the matrices are no longer tri-diagonal matrices which can not be solve by employing the Thomas algorithm. Thus, an alternative method, the RCM method is employed to decrease the dimension of the matrices and improve the computational efficiency 27 .

Numerical results and discussion
Effectiveness of the proposed unconditionally stable implementation. Before the investigation of GSs, the effectiveness and efficiency of the proposed CNDG-HO-PML algorithm is validated through a full-filled domain with graphene with the parameters of T 0 = 300 K , ω p = 2.42 × 10 13 s −1 , ε ∝ = 7 and ν = 1.93 × 10 11 s −1 , where T 0 is the environment temperature. The sketch picture of the computational domain is shown in Fig. 2. The computational domain has dimensions of 70 x × 70 z . At the boundaries of the x-direc- The waveform observed at the observation point obtained by different CFLNs are shown in Fig. 3. All of the waveforms are almost overlapped. In addition, it can be observed that although the computational accuracy degenerates with the increment of CFLNs, the proposed implementations when CFLN = 20 shows higher accuracy.
The accuracy of the calculation also affects by the absorbing performance of the PML regions. The absorbing performance of the PML regions can be reflected by the relative reflection error in the time domain which can be defined as is the test solution which can be observed directly from the observation point and E r x (t) is the reference solution which can be obtained by enlarging the computational domain to 7000 x × 7000 z and terminating by 128-cell-PML without changing the relative position between the source and observation point. During the calculation of the reference solution, the reflection waves has no effected on the waveform. Figure 4 shows the relative reflection error obtained by different algorithms with different CFLNs. The absorbing performance degenerates with the increment of CFLNs. The maximum value of the relative reflection error of FDTD-HO-PML, CNDG-HO-PML CFLN = 1, 10 and 20 are − 109 dB, − 96 dB and − 81 dB, respectively. Although the absorbing performance decreases 21 dB when CFLN = 20, the proposed scheme can still be employed in the practical engineering (usually below − 40 dB) 29 . The computational efficiency, memory consumption and time reduction are shown in Table 1. As is shown that computational efficiency decreases when CFLN = 1. The reason is that the two matrices obtained by RCM method must be calculated at each time step resulting in such condition. The computational efficiency can be further improved by employing larger CFLNs to reduce the iteration steps. It can be observed that when CFLN = 10 and 20, the CPU time reduces by 75.0% and 87.3%, respectively. This indicating the computational efficiency can be obviously improved by employing the proposed scheme. In conclusion, the proposed scheme shows considerable performance and accuracy. Thus, the following investigations are based on CNDG-HO-PML when CFLN = 20. The FDTD-HO-PML chosen for comparison. It should be noticed that the mesh size and the time step is no longer limited by the CFL condition. Thus, the physical size can be employed during the simulation. To ensure the accuracy of calculation, the thickness of single-layer GS is divided into 10 cells as x = z = = 0.34 nm. The time step can be obtained from t = √ 2c as 8.02 × 10 −19 s. Firstly, the wave transmission through GSs can be investigated through the maximum value of the waveform obtained from the observation point. Figure 6 shows the normalized transmitted energy versus different thickness obtained by FDTD-HO-PML and CNDG-HO-PML CFLN = 20. It can be observed that the curves are almost overlapped indicating they hold same accuracy. The CPU time of FDTD-HO-PML and CNDG-HO-PML CFLN = 20 are 78.3 s and 10.6 s, respectively. By employing the proposed scheme, the CPU time can be decreased by 86.5%. Thus, the effectiveness and the accuracy of the proposed scheme have been demonstrated.
In addition, it can be observed that the normalized transmitted energy decreases with the increment of the thickness. This indicates that the thicker sheet transmits less energy. This phenomenon can be reflected by the normalized reflection energy, shown in Fig. 7. It can be observed that the reflection energy increases with the increment of thickness. Especially, it can be founded that by adding the normalized reflection energy and the normalized transmitted energy, the resultant is no longer 1 at some cases. This indicates that a part of energy is lost by the sheet. The normalized loss energy is shown in Fig. 8. The results prove that thicker sheets losses much energy during the wave transmission.   www.nature.com/scientificreports/ To demonstrate the wave transmission in the frequency domain, the S parameters in the microwave theory is introduced during the simulation 30 . As defined in the microwave theory, the transmission can be reflected by S 21 . Thus, the port 1 is defined at the left side of sheet with the incidence plane wave. Port 2 is defined at the right side of the sheet. In this simulation, the S 21 is obtained by CNDG-HO-PML when CFLN = 20, as shown in Fig. 9. It can be observed that S 21 decreases with the increment of the thickness. At the frequency of 6.2 PHz, the graphene sheets show the least S parameters. At the low-frequency, the graphene sheets have higher S parameters indicating the less loss during the wave transmission. Figure 10 show the wavelength versus S 21 parameters obtained by different thickness. It can be observed that the S 21 becomes lower in small wavelength.

Extraordinary optical transmission through graphene sheet.
To investigate the extraordinary optical transmission phenomenon through the graphene sheet, a periodic graphene structure with slit is proposed. Figure 11 demonstrates the constitution of the whole structure and its single unit.  www.nature.com/scientificreports/ P G , D G , S G represent the length of a single unit in the z-direction, the width of graphene sheet in x-direction and the length of slit in the z-direction, respectively. The whole computational domain of single unit is divided into two regions including the total field (TF) and scattering field (SF) regions. The TF region is located on the right side of the incidence plane and the SF is on the left of the incidence plane. The length of the SF in the z-direction is 10 cells. The distance between incidence plane and sheet is 50 cells. The distance between observation plane and GS is 50 cells as well. P G is 100 cells in this simulation. The width of the graphene S G changes from 6 to 36 cells. The thickness of the graphene D G changes from 10 to 80 cells. The excitation source, mesh size and time step are the same as the numerical example above. The resultants are obtained by CNDG-HO-PML with CFLN = 20.
Compared with the graphene sheet without slit in the previously numerical example, due to the existence of slit, the interference of wave transmission results in the generation of the crest and trough in the observation plane. Thus, it should be noticed that the energy from each mesh grid in the observation plane should be employed during the calculation. Figures 12, 13 and 14 show the normalized transmitted energy, normalized reflected energy, and normalized lost energy versus different thickness with the slit width of 6 cells. It can be observed that the transmitted energy decreases with the increment of thickness. The normalized reflected energy increases at the same time. Meanwhile, compared with the normalized lost energy in the graphene sheet, it decreases in this numerical example. The reason is that due to the existence of the slit, less energy is consumed during the wave transmission of the graphene. In addition, compared with the single graphene sheet, the normalized transmitted energy through the graphene grating becomes larger. The electric current that pass through  www.nature.com/scientificreports/ the slit causes such phenomenon. Figures 15, 16 and 17 show the normalized transmitted energy, normalized reflected energy, and normalized lost energy versus different slit width with the 8 layers thick sheet. It can be observed that the transmitted energy increases with the increment of the slit width. The reason is that large slit width can transmit much power compared with the narrow slit. Meanwhile, the consumption energy of graphene decreases with the increment of width slit indicating that less energy is cost during the wave transmission. Following the same step as the numerical example above, the S 21 parameters can be employed in evaluating the performance and investigating the properties of graphene. Figures 18 and 19 show the S 21 parameters versus frequency and wavelength with different slit width under 8-layer-sheet, respectively. It can be observed that the S parameters increases obviously with the increment of slit width. This indicates that the wave is much easier to transmit with the increment of slit width before the 6 PHz frequency.
As is demonstrate in the previously that the wave is no longer uniform at the observation plane. The reason is that the interference of electric current occurs when wave transmitting through the slit. The crest and the trough are formed at the observation point. Thus, it is important to investigate the crest and the trough of the waveform.
In this simulation, the parameter η = |max {H t (t)}| 2 |max {H r (t)}| 2 is defined as the ratio of the field intensity at the observation plane, where H t (t) is the magnetic field different situations, H r (t) is the magnetic field under the circumstance of single sheet [31][32][33][34] . Figures 20 and 21 demonstrate the ratio versus thickness and silt width, respectively. As is shown in Fig. 20 that the influences of graphene gratings on the magnitude are significantly placid when the graphene thickness is larger than 4 layers. As shown in Fig. 21, the ratio of the wave intensity increases linearly with slit widths. This can be explained that when the thickness of the GS increases, part of the electric current pass through the slit. The portion which traverse the GS suffers more loss. Therefore, the effect of wave interference that causes by the  www.nature.com/scientificreports/ electric current becomes weak. In conclusion, the ratio between crest and trough is coincidence not only with the trend of energy but also with the S 21 parameters.

Conclusions
Based on the CNDG algorithm, the Drude-like model and higher order concept, an unconditionally stable CNDG-HO-PML algorithm is presented to investigate the wave propagation through graphene gratings. The Maxwell's equations are applied to the model to study the transmission characteristics. The ADE-and PLRC-FDTD methods with faster calculation and less memory is used to solve the partial differential equations for time stepping of the electromagnetic fields. The effectiveness and efficiency are testified through a full-filled structure, a single sheet and the sheet with slit. Based on numerical analysis, the obtained result shows that the