## Abstract

We propose a simple semi analytical model that allows to compute the transmittance and reflectance of a one dimensional subwavelength graphene strip grating under an external static magnetic field. In this model graphene is treated as an anisotropic layer with atomic thickness and a frequency dependent complex permittivity tensor. The model is based on an effective medium approach (EMA) and a rigorous phase correction. The scattering matrix approach is also used to take into account the different resonant phenomena occurring in the structure. The approach is validated against the Polynomial Modal Method (PMM) through numerical examples.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

## 1. Introduction

Graphene, a two dimensional (2D) material made of carbon atoms arranged in a honeycomb lattice [1], has recently emerged in the plasmonic field as a potential and promising plasmonic material. This 2D material which, for decades, was predicted to be experimentally unviable, is a gapless semi conductor with unique and unusual optical, mechanical and thermal properties [2]. In addition to the aforementioned properties and due to its unusual optical conductivity which can be dynamically tuned by electrostatic bias or chemical doping [3], doped graphene can support surface plasmons polaritons modes (GSPPs) in the mid-infrared and the terahertz regions with attractive and fascinating properties such as tunability, extreme confinement and relatively low losses. GSPPs have been investigated in several studies [4–7] and they have shown very important potential in various applications. In particular, recent studies have been focused on surface plasmons polaritons over graphene nanoribbons arrays. For example, in [4] the plasmonic coupling between two crossed layers of graphene nanoribbons has been exploited to achieve plasmonically induced transparency (PIT).

Moreover, the conductivity of graphene can be also tuned through a magnetic bias. In fact, when a perpendicular magnetic field is turned on, graphene becomes gyrotropic and exhibits non reciprocal properties [8]. This leads to many magneto-optical (MO) phenomena such as giant Faraday rotation [9] and Kerr effects [10]. Furthermore, in the presence of a static magnetic field, plasmons and cyclotron excitations hybridize which leads to the appearance of graphene surface magnetoplasmons polaritons (GSMPPs). One of the great advantages of graphene over conventional plasmonic materials, is the high sensitivity of magnetoplasmons to external magnetic fields because the cyclotron frequency is comparable to the plasmon frequency [11]. An additional advantage is that the cyclotron mass of the charge carriers is small and can be tuned by doping. Moreover while cryogenic temperatures are required to observe magneto-optical effects in 2DEGs (two-dimensional electron gases), in graphene, they already arise at room temperature. All this makes graphene one of the most promising candidates for use in tunable plasmonic nonreciprocal devices in the microwave and terahertz regimes such as optical isolators [12] and absorbers [13]. Recently, graphene surface magneto-plasmons have been widely studied in several theoretical and experimental works under various forms and configurations including extended layers of graphene [14], nano-ribbons [15,16] and antidots [17]. It has been shown that the excitation of GSMPs significantly enhances the MO effects in graphene structures and strong magnetoplasmonics resonances have been observed in a patterned arrays of graphene antidots. As a result, it is of great interest to analytically and numerically investigate the optical properties of patterned anisotropic graphene under a magnetic bias. Among these structures, periodic magnetically biased graphene ribbons arrays have attracted a great deal of research interest and became one of the most studied graphene based structures [4].

Several and various numerical methods have been employed for modelling this structure such as the Fourier Modal Method (FMM) [16], the Finite Difference Time-Domain (FDTD) method [18] and the Discontinuous Galerkin Time Domain (DGTD) approach [19]. However, it has been shown that due to the ultra-thin nature of graphene, these numerical methods suffer from slow convergence rates and require too much computational effort. In addition to these numerical methods, an analytical approach has been proposed in [20] to study the electromagnetic behavior of strongly anisotropic metasurfaces implemented by a subwavelength array of graphene strips. It is based on an Effective Medium Approach (EMA) applied to the conductivity tensor of graphene in the electrostatic limit. Nevertheless, this method provides accurate results only for structures in the extreme subwavelength regime (with infinitesimally small periodicity as compared to the wavelength). More recently, an analytical method has been proposed in [21] to analyse magnetically-biased arrays of graphene ribbons by solving the integral equations governing the induced surface currents on the graphene ribbons. The major drawback of this method is that it does not take into account the interaction between neighboring ribbons which makes it less accurate when one is dealing with strip gratings with small gaps.

To the best of our knowledge, no analytical or semi-analytical method exists for the prediction of the electromagnetic response of a one dimensional graphene strip grating in the presence of an external magnetic field when the graphene is modelled as a bulk material using an Effective Medium Approach. Accordingly, in the present work, we provide a simple and fast semi-analytical model for the investigation of the propagation properties of a 1D magnetically biased graphene strip grating. It is based on an EMA and a phase correction. In our model, graphene is regarded as an extremely thin anisotropic layer with a finite thickness and unlike the approach proposed in [20] the EMA will be applied to its permittivity anisotropic tensor instead of the conductivity tensor. Note that the present approach has two main advantages over the above reported methods: (i) although this method is developed to deal with the general case of a 1D magnetically biased graphene strip grating, it is still valid even for a non magnetized structure where the magnetic bias is absent, and (ii) this method can efficiently and accurately treat subwavelength graphene strip gratings with either small gaps or large gaps.

This paper is organized as follows: In section 2., we introduce the physical system. Then, in the following section, the proposed semi analytical method is be explained and described in detail. Finally, in the last section, some numerical results are given and discussed to illustrate the effectiveness of our approach. The accuracy of the proposed model will be validated through comparisons with both the Polynomial Modal Method (PMM) [22,–24] and the Effective Medium Approach proposed in [20].

## 2. Physical system

The structure under study is depicted in Fig. 1. It consists of a periodic 1D graphene ribbons array (along the $x$ direction with width $a$ and period $d=a+s$) biased with a magnetostatic field $\mathbf {B}$. The graphene ribbons array is suspended in free space and occupies the $xy$ plane. It is enlightened from the upper medium $(z>0)$ by a plane wave under normal incidence. Here we consider the Faraday geometry where the external static magnetic field $\mathbf {B}$ is perpendicular to the array. Under the influence of such a magnetostatic bias, graphene will exhibit an induced anisotropy leading to a surface magneto-optical conductivity tensor of the form [25]:

For highly doped graphene i.e $\mu _c >>\hbar \omega$ and $\mu _c >> k_BT$, the conductivity is dominated by the intraband contribution and thus the general quantum mechanical conductivity reduces to the semi-classical Drude model form [26]:

## 3. Theoretical Method

To simply model the electromagnetic response of the studied structure and to understand its physical properties we use the effective medium approach which is an intuitive and effective technique allowing to obtain the macroscopic averaged response of an array of subwavelength resonators [27–29]. Since, the period $d$ of the graphene strip grating is small compared to the wavelength of the incident wave ($d<<\lambda$), it can be approximated by a homogeneous anisotropic effective layer with an equivalent effective permittivity tensor $\tilde {\varepsilon }$. To find $\tilde {\varepsilon }$ we first consider the constitutive relations linking the displacement to the electric field components:

According to these relations, we can rewrite $D^y$ and $E_x$ in terms of $D^x$ and $E_y$ as:

Averaging all quantities involving local fields and densities over the period $d$, and assuming that the tangential component of the electromagnetic field $E_y$ and the normal component of the density $D^x$ are continuous at any coordinate curve $x=constant$, we obtain:

Treating the grating as an effective medium requires that these macroscopic fields satisfy the constitutive relations given by Eq. (5), such that:

Comparing with Eq. (8), we can deduce the effective medium permittivity tensor components :

At this stage, it is important to emphasize that, close to a single strip, the wave sees the latter as surrounded on both sides by an effective medium. Thus, the structure can be regarded as a strip grating inlayed into an effective medium. For that reason and in order to take into account the resonance phenomena occurring in the grating, the original periodic structure will be modelled as a subwavelength periodic graphene strip grating with the same dimensions but instead of vacuum the gaps between neighboring strips will be filled by the effective medium (see Fig. 2). Let us now take a look at the different modes living in the new approximate structure. From Fig. 2, one can distinguish two kinds of modes: the first is the most slowly decaying evanescent mode of the sub-wavelength periodic structure given by the couple $(\tilde {\alpha },\tilde {\gamma })$ and the second one corresponds to the surface magnetoplasmon polariton of graphene (GSMP) guided by the graphene strip $(\alpha,\gamma )$. The challenge now is to successfully introduce the coupling between these modes. One powerful and insightful way to accomplish this is to adopt an analysis based on a scattering parameters of a two-port network. In the proposed analysis, we will use the scattering matrix formalism to describe the optical characteristics of the structure (see Refs. [30,31] for more details). The total scattering matrix associated to one period can be written as:

As mentioned above, $\widetilde {\alpha }$ denotes the in-plane effective index of the most slowly decaying evanescent mode, of the sub-wavelength periodic structure, and can be approximated by:

In the case of a periodic subwavelength strip grating with a small gap size ($s<<d$), we have $\widetilde {\varepsilon }=\widetilde {\varepsilon }_{xx}>> \varepsilon _s$ which yields:

Thus, in this particular case the effective index $\nu _{\text {eff}}$ can be straightforwardly given by:

In the proposed model, we assume that the phase term $\phi$ appearing in Eq. (13) is composed of two contributions:

The first term $\phi _a$ handles the propagation along the $a-width$ graphene strip. While the second term $\phi _c$ is related to the contribution of the periodicity of the strips array. In this term $\theta _c$ is a corrective phase required to take into account the phase delay. To get an idea about the behavior of the phase of the electric field in the considered structure, we plot, in Fig. 3, its spatial distribution (calculated with the PMM) at the resonance wavelength. One can clearly see that the phase of the electric field reveals a jump at the ends of the graphene strips.

Note that, unlike in the phase term $\phi _a$, the phase delay $\theta _c$ cannot be provided analytically. However it may be estimated thanks to a rigorous simulation. In this context, several numerical simulations showed that this phase depends only on the geometric parameters, through the ratio $a/s$. For example, for $s/a=1/4$, it is found that $\theta _c(s/a)=\pi /2$ whatever the value of $a$ and the magnitude of $B$. In the following, we seek a robust estimation of this phase. To this end, we write $\theta _c$ as a series of the ration $\zeta =s/a$ as follows :

Then, the weighted coefficients of the above series are obtained numerically as:

To show the reliability and the efficiency of the proposed approximate model, we plot in Fig. 4 the behaviour of $\theta _c$ as a function of the ratio $s/a$ for the following numerical parameters: $\lambda =30.35\mu m$, $a=0.735\mu m$, $B=0T$, $\mu _c=0.3eV$. The numerical results are compared to those computed by the approximate model. The red solid line shows the results obtained numerically and the results obtained using the approximate model are depicted by the blue dashed and the black dotted lines. One can see that the blue line representing $\theta _c^{(1)}$ fits perfectly the red line for $\zeta \geq 0.5$ while the black one corresponding to $\theta _c^{(2)}$ matches the latter for $\zeta >0.5$. Hence, we conclude that in this case, the behavior of the phase delay $\theta _c$ can be efficiently predicted by the above proposed approximate model. Note that although we have presented here only the case of the structure without a magnetic bias, the model has been also verified for $B\neq 0$ and it is found to remains valid and efficient for all values of $B$ less than $7T$.

After having defined all the parameters required for the determination the S-matrix, let us now find the resonance modes of the system which are obtained by setting $det(S^{-1})=0$, since the S matrix is calculated for a symmetric structure (the input and the output media are the same), it is found that the zeros of $det(S^{-1})$ are the same as those of $det(S)$:

We get then two categories of solutions:

To find out which one corresponds to the resonant modes of the system, we compare each of them with the curves computed by the PMM method. As seen from Fig. 5, the resonance frequency defined by the class of solutions satisfying to $\Delta _+=S_{11}(\omega )+S_{12}(\omega )$ matches with the resonant dip in the transmission spectrum of the structure. Thus, we can set:

where $t$ and $r$ denote the transmissions and reflection coefficients resulting from the resonance phenomena.In addition to this contribution, we should add a second term that takes into account the contribution of a single strip of graphene. Therefore the transmission and the reflection spectra of the system can take the following forms:

where $T_s$ and $R_s$ are respectively the transmission and reflection coefficients through a graphene mono-layer taken from [32]. Now, our semi-analytical model is complete and in the next section we are going to test and validate it.## 4. Numerical Results and discussion

In this section, some numerical examples are provided to validate the proposed semi analytical model and demonstrate its viability. For this purpose, we investigate the transmission and reflection properties for two different cases: with an external magnetic field and without. First, let us consider a graphene strip grating with the following parameters: $a=0.9\mu$m, $s=0.45 \mu$m, $B=3$T, $\mu _c=0.6$eV, $\tau =0.6$ps. The reflection and transmission spectra are shown in Fig. 6. Their real and imaginary parts are also plotted in Fig. 7. These results are compared with the effective medium approach given in Ref. [20] and the PMM method. We can see a good agreement between the PMM results and those obtained with our proposed model. By contrast, the comparison with the effective conductivity model shows that the latter cannot reproduce the magneto-optical properties of the studied graphene grating.

As a second example, we consider the case of a graphene strip grating without an external magnetic field. The strip width is $a=0.25 \mu$m and the periodicity is assumed to be $d=0.35 \mu$m. The chemical potential of the graphene strips is set to $\mu _c=0.4 e$V and the relaxation time $\tau$ is taken to be $0.4$ps. Figure 8 shows the reflection and transmission spectra computed through the proposed semi-analytical model (red solid lines), the PMM model (blue solid lines) and the effective conductivity approach proposed in [20] (dashed black lines). Similarly to the first example, an excellent agreement is observed between our model and the PMM results while the effective conductivity approach [20] is not able to predict the magneto-optical response of this structure. From these results, we can conclude that our semi-analytical model can efficiently and successfully treat the general case of a graphene strip grating biased with an external magnetic field as well as the particular case where no magnetic field is applied.

In the following, we are going to test and check the performance and the accuracy of our approach by examining it as a function of grating and graphene parameters including the geometric parameter $a/d$, the chemical potential and the magnetic bias. Let us begin a test for different filling factors $a/d$. To do so, we set $d=1.35\mu m$ and plot the wavelength of the first resonance as a function of filling factor for two values of the magnetic field $B=0 T$ and $B=3 T$ in Fig. 9. The results obtained with our approach (red pluses) are compared with those performed by the PMM method (blue asterisks). One can clearly see that the results correlate well with each other for all filling factors. In particular, an excellent agreement is observed for $a/d$ smaller than $0.4$ and for $a/d$ greater than $0.7$. For the remaining filling factors $0.4<a/d<0.7$ the relative error does not exceed $4\%$. For example for $a/d=0.5736$, the relative error is equal to $2.04\%$ for $B=0 T$ and $3.12\%$ for $B=3T$. Hence, the proposed model provides accurate results for all filling factors and can efficiently deal with gratings with a small gap as well as a large gap. It is important to note that the grating with small gaps constitute a critical case because of the strong interaction between neighboring strips and that has to be taken into account to ensure an accurate model.

We look now at the accuracy of the proposed approach as a function of the magnetic field. Figure 10 shows the wavelength of the first resonance as a function of the magnetic field for three values of the filling factor $a/d=0.2$, $a/d=0.5$ and $a/d=0.8$. It is shown that for $a/d=0.2$ the results obtained from our model match well with those computed with the PMM method for all magnetic intensities between [0,10]T. While for $a/d=0.5$ and $a/d=0.8$ the results correlate well up to $7$T . For $7 T<B<10 T$, they shift away from those of the PMM method. This shift increases slightly while increasing the magnitude of $B$. For instance for $B=8.42$T, the relative error is equal to $2\% (a/d=0.9)$ and it increases for $B=10$T, to reach $3.41\%$. This could be explained by the fact that the model proposed in (20) and (21) to approximate the phase $\theta _c$ is only valid for $B \leq 7T$. So, to remove or minimize the difference between the two methods, we should numerically adjust the coefficients $a^{(1)}_n$ and $a^{(2)}_n$ to find the best fit model for the phase. Obviously this will affect the accuracy of our model which becomes less accurate for $B > 7T$. But, it does not raise a considerable and real problem since, from a practical point of view, the magnetic bias usually used does not exceed $7T$ [33]. We end this section by examining the stability of the proposed model when varying the chemical potential of the graphene strips. Figure 11 shows the wavelength of the first resonance as a function of the chemical potential. The reflection and transmission spectra for different chemical potential values namely for $\mu _c=0.4 ,0.6$ and $0.8$eV are also plotted in Fig. 12. In these figures , we compare the results computed by our model with those obtained from the PMM Method. We can see from Fig. 11 that for all chemical potentials the first resonance is successfully predicted by the proposed method. From Fig. 12, it is shown that the transmittance and the reflectance obtained with the PMM are well fitted by the proposed approach. One can also see that by decreasing $\mu _c$, the accuracy of the proposed method with respect to the amplitude of the transmittance and specifically of the reflectance is slightly decreased. This last behavior is expected. In fact as seen from Eq. (25), in our calculations, we do not take into account the absorption features within the structure. Therefore, when the chemical potential decreases, the losses in graphene increase and then the model becomes less accurate. However, it should be noted that the difference between the two curves is relatively small and thus this cannot be considered as a significant limitation of the proposed semi analytical approach.

## 5. Conclusion

In this paper, we presented a simple, accurate and fast semi analytical model to compute the transmittance and reflectance of a magnetically-biased subwavelength graphene strip grating. It is based on the effective medium approach when the graphene is modelled as an anisotropic layer with atomic thickness and a frequency-dependent and complex permittivity tensor. The proposed model is validated against the PMM method by numerical examples. The performance and the accuracy of this method has been also tested and checked as a function of the graphene and grating parameters. We first showed that this method is a general approach that can efficiently treat the structure when the magnetic field is applied and still valid even in the particular case of a structure without an external magnetic field which is not possible with the analytical methods already encountered in the literature. In addition, the study of the accuracy of the method as a function of the filling factor of the grating revealed that the proposed model can efficiently deal with gratings with small gaps as well as with large gaps. Specifically, this method provides precise and reliable results in the case of gratings with small gaps which constitute a complicated case to deal with. Through this last structure it is possible to achieve strong reflection ( near $90\%$) which could be of interest in devices that require high reflectivity.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **A. K. Geim and K. S. Novoselov, “The Rise of Graphene,” Nat. Mater. **6**(3), 183–191 (2007). [CrossRef]

**2. **F. Bonaccorso and Z. Sun and T. Hasan and A. Ferrari, “Graphene photonics and optoelectronics,” Nat. Mater. **4**(9), 611–622 (2010). [CrossRef]

**3. **O. J. Sik and K. K. Nam and Y. G. Young, “Graphene Doping Methods and Device Applications,” J. Nanosci. Nanotechnol. **14**(2), 1120–1133 (2014). [CrossRef]

**4. **S.-X. Xia, X. Zhai, L.-L. Wang, and S-C. Wen, “Plasmonically induced transparency in double-layered graphene nanoribbons,” Photonics Res. **6**(7), 692–702 (2018). [CrossRef]

**5. **S.-X. Xia, X. Zhai, L.-L. Wang, and S.-C. Wen, “Plasmonically induced transparency in in-plane isotropic and anisotropic 2D materials,” Opt. Express **28**(6), 7980–8002 (2020). [CrossRef]

**6. **M. Ben Rhouma, M. Oueslati, and B. Guizal, “Surface plasmons on a doped graphene sheet with periodically modulated conductivity,” Superlattices Microstruct. **96**, 212–219 (2016). [CrossRef]

**7. **M. B. Rhouma, O. Meherzi, and B. Guizal, “Strong coupling between a plasmonic waveguide and graphene surface plasmons,” J. Opt. Soc. Am. B **34**(4), 884–890 (2017). [CrossRef]

**8. **D. L. Sounas, H. S. Skulason, H. V. Nguyen, A. Guermoune, M. Siaj, T. Szkopek, and C. Caloz, “Faraday rotation in magnetically biased graphene at microwave frequencies,” Appl. Phys. Lett. **102**(19), 191901 (2013). [CrossRef]

**9. **I. Crassee, J. Levallois, A. L. Walter, M. Ostler, A. Bostwick, E. Rotenberg, T. Seyller, D. van der Marel, and A. B. Kuzmenko, “Giant Faraday rotation in single- and multilayer graphene,” Nat. Phys. **7**(1), 48–51 (2011). [CrossRef]

**10. **T. Yoshino, “Giant Faraday rotation in single- and multilayer graphene,” J. Opt. Soc. Am. B **30**(5), 1085–1091 (2013). [CrossRef]

**11. **T. Yoshino, “Quasiclassical cyclotron resonance of Dirac fermions in highly doped graphene,” Phys. Rev. B **82**(16), 165305 (2010). [CrossRef]

**12. **X. Lin, Z. Wang, F. Gao, B. Zhang, and H. Chen, “Atomically thin nonreciprocal optical isolation,” Sci. Rep. **4**(1), 4190 (2015). [CrossRef]

**13. **M. Wang, Y. Wang, M. Pu, C. Hu, X. Wu, Z. Zhao, and X. Luo, “Circular dichroism of graphene-based absorber in static magnetic field,” J. Appl. Phys. **115**(15), 154312 (2014). [CrossRef]

**14. **B. Hu, J. Tao, Y. Zhang, and Q. J. Wang, “Magneto-plasmonics in graphene-dielectric sandwich,” J. Appl. Phys. **22**(18), 21727–21738 (2014). [CrossRef]

**15. **J. M. Poumirol, W. Yu, X. Chen, C. Berger, W. A. de Heer, M. L. Smith, T. Ohta, W. Pan, M. O. Goerbig, D. Smirnov, and Z. Jiang, “Magnetoplasmons in Quasineutral Epitaxial Graphene Nanoribbons,” Phys. Rev. Lett. **110**(24), 246803 (2013). [CrossRef]

**16. **M. Tymchenko, A. Yu. Nikitin, and L. Martin-Moreno, “Faraday Rotation Due to Excitation of Magnetoplasmons in Graphene Microribbons,” ACS Nano **7**(11), 9780–9787 (2013). [CrossRef]

**17. **J.-M. Poumirol, P. Q. Liu, T. M. Slipchenko, A. Y. Nikitin, L. Martin-Moreno, J. Faist, and A. B. Kuzmenko, “Electrically controlled terahertz magneto-optical phenomena in continuous and patterned graphene,” Nat. Commun. **8**(1), 14626 (2017). [CrossRef]

**18. **X. Wang, W. Yin, and Z. Chen, “Matrix Exponential FDTD Modeling of Magnetized Graphene Sheet,” Antennas Wirel. Propag. Lett. **12**, 1129–1132 (2013). [CrossRef]

**19. **P. Li and L. J. Jiang, “Modeling of Magnetized Graphene From Microwave to THz Range by DGTD With a Scalar RBC and an ADE,” IEEE Trans. Antennas Propag. **63**(10), 4458–4467 (2015). [CrossRef]

**20. **J. S. Gomez-Diaz and A. Alu, “Magnetically-biased graphene-based hyperbolic metasurfaces,” 2016 IEEE International Symposium on Antennas and Propagation (APSURSI), 359–360 (2016).

**21. **M. Rahmanzadeh, B. Rejaei, M. Memarian, and A. Khavasi, “Analytical and rigorous method for analysis of an array of magnetically-biased graphene ribbons,” Opt. Express **27**(20), 28395–28409 (2019). [CrossRef]

**22. **K. Edee, “Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings,” J. Opt. Soc. Am. A **28**(10), 2006–2013 (2011). [CrossRef]

**23. **K. Edee, I. Fenniche, G. Granet, and B. Guizal, “Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings: weighting function, convergence and stability,” PIER **133**, 17–35 (2013). [CrossRef]

**24. **K. Edee and J.-P. Plumey, “Numerical scheme for the modal method based on subsectional Gegenbauer polynomial expansion: application to biperiodic binary grating,” J. Opt. Soc. Am. A **32**(3), 402–410 (2015). [CrossRef]

**25. **G. W. Hanson, “Dyadic Green functions for an anisotropic, nonlocal model of biased graphenne,” IEEE Trans. Antennas Propagat. **56**(3), 747–757 (2008). [CrossRef]

**26. **A. Ferreira, J. Viana-Gomes, Yu. V. Bludov, V. Pereira, N. M. R. Peres, and A. H. Castro Neto, “Faraday effect in graphene enclosed in an optical cavity and the equation of motion method for the study of magneto-optical transport in solids,” Phys. Rev. B **84**(23), 235410 (2011). [CrossRef]

**27. **J. M. Bell and G. H. Derrick and R. C. McPhedran, “Diffraction Gratings in the Quasi-static Limit,” Opt. Acta **29**(11), 1475–1489 (1982). [CrossRef]

**28. **P. Lalanne and D. Lemercier-lalanne, “On the effective medium theory of subwavelength periodic structures,” J. Mod. Opt. **43**(10), 2063–2085 (1996). [CrossRef]

**29. **P. Lalanne and D. Lemercier-lalanne, “Depth dependence of the effective properties of subwavelength gratings,” J. Opt. Soc. Am. A **14**(2), 450–459 (1997). [CrossRef]

**30. **K. Edee, “Understanding the plane wave excitation of the metal-insulator-metal gap plasmon mode of a nanoribbons periodic array: role of insulator-metal-insulator lattice mode,” OSA Continuum **2**(2), 389–399 (2019). [CrossRef]

**31. **K. Edee, M. Benrhouma, M. Antezza, J. A. Fan, and B. Guizal, “Coupling between subwavelength nano-slit lattice modes and metal-insulator-graphene cavity modes: a semi-analytical model,” OSA Continuum **2**(4), 1296–1309 (2019). [CrossRef]

**32. **D. L. Sounas and C. Caloz, “Gyrotropy and Nonreciprocity of Graphene for Microwave Applications,” IEEE Trans. Microwave Theory Tech. **60**(4), 901–914 (2012). [CrossRef]

**33. **I. Crassee, M. Orlita, M. Potemski, A. L. Walter, M. Ostler, Th. Seyller, I. Gaponenko, J. Chen, and A. B. Kuzmenko, “Intrinsic Terahertz Plasmons and Magnetoplasmons in Large Scale Monolayer Graphene,” Nano Lett. **12**(5), 2470–2474 (2012). [CrossRef]