## Abstract

A sheet of graphene under magnetic bias attains anisotropic surface conductivity, opening the door for realizing compact devices such as Faraday rotators, isolators and circulators. In this paper, an accurate and analytical method is proposed for a periodic array of graphene ribbons under magnetic bias. The method is based on integral equations governing the induced surface currents on the coplanar array of graphene ribbons. For subwavelength size ribbons subjected to an incident plane wave, the current distribution is derived leading to analytical expressions for the reflection/transmission coefficients. The results obtained are in excellent agreement with full-wave simulations and predict resonant spectral effects that cannot be accounted for by existing semi-analytical methods. Finally, we extract an analytical, closed form solution for the Faraday rotation of magnetically-biased graphene ribbons. In contrast to previous studies, this paper presents a fast, precise and reliable technique for analyzing magnetically-biased array of graphene ribbons, which are one of the most popular graphene-based structures.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Graphene is a two-dimensional (2D) nanomaterial that is composed of a single layer of carbon atoms organized into a hexagonal lattice [1]. In the past decade, graphene-based devices have become a prime focus of research due to their unique electrical, optical, thermal, and mechanical properties [2],[3],[4]. The electrical conductivity of graphene can change by an external DC voltage or chemical doping. This tunability provides additional flexibility in the design of diverse functionalities. Therefore, numerous graphene-based applications have been proposed from microwave to terahertz regions [5–14]. Under an external magnetic bias, the surface conductivity of graphene becomes a tensor, and it will also have non-reciprocal properties [15]. These extraordinary properties make graphene a promising material for many novel applications, such as [16, 17] , nonreciprocal couplers [18], optical modulators [19], and Faraday rotators [20], in the terahertz regime.

Due to strong interaction of graphene with electromagnetic (EM) fields, patterned magnetically-biased graphene has been assessed up to now for several applications [21–23]. For example, in [23] it has been shown that the Faraday rotation of graphene micro ribbons can be controlled dynamically and that a giant Faraday rotation can be obtained, the frequency of which corresponds to the width of ribbons. But a closed-form relationship between the width of the graphene ribbons and the resonance frequency was not presented.

Consequently, the analytical and numerical study of continuous and patterned anisotropic graphene is currently of high interest. A continuous graphene sheet can be easily analyzed by applying appropriate boundary conditions. Sounas and Caloz studied the transmission and reflection properties of magnetically-biased graphene sheet using the general anisotropic conductivity tensor [24]. To the best of the authors’ knowledge, no rigorous and analytical solution has yet been presented on patterned anisotropic graphene under magnetic bias. Although certain numerical methods, including the Fourier modal-based [23], the finite difference time-domain method [25], and discontinuous Galerkin time domain (DGTD) [26] have been used to numerically solve this problem, they suffer from a slow convergence rate and more importantly, are not analytic results. Furthermore, some efforts have been made to decrease the computation time [27, 28]; however, they still require higher acceleration. In addition, a simple method for the analysis of magnetic-biased graphene ribbons based on the effective medium approach was proposed in [29]; however, as it will be shown, this method does not have accuracy and is only valid for the extreme subwavelength regime. Hence a rigorous and analytical analysis of magnetically-biased Graphene-based meta-surfaces is of paramount importance.

In this paper, a rigorous analysis is performed through solving the integral equations governing the induced surface currents on the graphene ribbons to accurately predict the electromagnetic performance of an array of magnetically-biased graphene ribbons. The paper is organized as follows: In Section (2) we shall first study the scattering of a TE/TM electromagnetic wave by a single magnetically-biased graphene ribbon. Then an analytical expression is extracted for the surface current density induced within a periodic array of magnetically-biased ribbons in the subwavelength regime. In Section (3), the reflection/transmission coefficient of the structure is calculated and the proposed method is validated against full-wave simulations. We show how plasmonic resonances across the graphene ribbons enter the reflection and transmission coefficients of such a structure. The limitations of our proposed method is discussed. Moreover, an analytical expression is presented for Faraday rotation of magnetically-biased graphene ribbons, as one of the important phenomena in optics and photonics.

## 2. Surface currents on a perodic array of anisotropic graphene ribbons

We shall consider a periodic array of 1D graphene ribbons under magnetic bias as shown in Fig. 1. Each ribbon has a width of *w*. The periodicity of the array is *D* in the x-direction. The graphene ribbons are placed on the *x* − *y* plane (infinite along *y*) and a static magnetic field *B*_{0} is applied perpendicular to the ribbons. A normal EM plane wave is incident on the arrangement. Due to the incident wave, a surface current density with components *J _{x}* and

*J*will be induced onthe ribbons. Since there is no variation in

_{y}*y*-direction, these induced currents are functions of

*x*only (i.e. ${\partial}_{y}{J}_{x,y}=0$).

#### 2.1. Electromagnetic parameters of magnetically-biased graphene

As noted earlier, graphene is a 2-D nanomaterial consisting of a single layer of carbon
atoms organized within a hexagonal lattice. Assume that the graphene
sheet is on the *x* − *y* plane
and that a magnetic field is also established perpendicular to the
graphene sheet ($\overrightarrow{B}={B}_{0}\widehat{z}$). The surface conductivity of this
graphene sheet is represented by a tensor (${\overline{\overline{\sigma}}}_{s}$) that is obtained from microscopic,
semi-classical, and quantum mechanical considerations as
[15]

*E*is the Fermi level energy and

_{f}*ω*is the angular frequency,

*τ*is the relaxation time, and

*T*is the Kelvin temperature ($T=300K$). For highly doped graphene, i.e. ${E}_{f}>>\hslash \omega $ and ${E}_{f}>>{k}_{B}T$ (

*k*= Boltzmann constant, ℏ = reduced Planck’s constant), (1) can be simplified by a Drude-like model as [30]

_{B}*μ*is the carrier mobility and ${n}_{s}={E}_{f}^{2}/\pi {\hslash}^{2}{v}_{f}^{2}$ is the carrier density [24]. Carrier mobility is assumed to be within the range of $0.1$ to $6\text{}{m}^{2}/V$, which is practical and achievable for graphene with different substrates at room temperature [31, 32]. It should also be noted that the time harmonic of the form $\mathrm{exp}\text{}(j\omega t)$ is assumed throughout this paper.

#### 2.2. Scattering from a single ribbon

Let us first consider a single anisotropic graphene ribbon in free space. The integral equation governing the induced surface currents on the single ribbon can be expressed as [33]

*x*- and

*y*-components of the incident (external) electric field, and

*J*and

_{x}*J*. In these equations

_{y}Assuming the strip to be narrow (${k}_{0}w\ll 1$), we next adopt the quasi-static approximation [34]. As a result, the second term on the right hand side of Eq. (3) becomes negligible compared to the first term. Moreover, the Green’s function Eq. (4) can be approximated by

From Eq. (3) and Eq. (3) one obtains

It must be noted that *σ*_{0} is the graphene conductivity in the absence of the external magnetic field.

Assume, for the moment, that $J\left(x\right)$ is known. Equation (6) can then be viewed as can integral equation for the function ${J}_{y}\left(x\right)$. This integral equation has a complicated, but exact solution. However, in the limit where

the second term on the left hand side of Eq. (6) becomes negligible and the solution is trivially given byTo see why Eq. (8) is satisfied, note that

where ${k}_{P}=-2j\omega {\epsilon}_{0}/{\sigma}_{0}$ is the wave number of graphene plasmons in the absence of an external magnetic field. Since ${k}_{0}w\ll 1$ and ${k}_{0}\ll {k}_{P}$ at frequencies of interest, Eq. (8) is satisfied.Let us now return to Eq. (3) and combine it with Eq. (3) to arrive at an equation for ${J}_{x}\left(x\right)$. However, one may simplify the resulting equation by neglecting the second term on the right hand side of Eq. (3). This term is small compared to the first term for a narrow ribbon. The integral equation governing *J _{x}* thus becomes

The solution to this equation has been described in detail in [34]. The only difference with the work in [34] is the appearance of the second term on the right hand side of Eq. (11).

#### 2.3. Scattering from an infinite array of ribbons

Now consider a periodic array of coplanar graphene ribbons with period *D* as shown in Fig. 1. The ribbon array is subject to a plane electromagnetic wave of arbitrary polarization that is normally incident on the array plane. Equation (3) that locally relates the ribbon surface currents to the total electric field is still applicable. However, on each ribbon, ${E}_{x,y}^{T}$ must now account for the total electric field that is produced by all ribbons. Fortunately, in case of normal incidence on an infinite number of array elements, the distributions of surface currents on all ribbons are identical. As a result, on each individual ribbon, Eq. (3) and Eq. (3) remain provided that the Green’s function Eq. (4) is replaced by

As in the case of a single ribbon, we first treat the equation for *J _{y}*. However, the interaction between ribbons, manifested in the periodic Green’s function Eq. (12), leads to an important modification of the approximation Eq. (5). For sub-wavelength arrays where ${k}_{0}D\ll 1$, a lengthy calculation presented in the appendix shows that

The interaction between sub-wavelength ribbons leads to a Green’s function whose leading term is of the order ${\left({k}_{0}D\right)}^{-1}$. Unlike the case of a single ribbon, the second term in the integral equation Eq. (6) cannot be neglected anymore. By replacing *G*_{0} with Eq. (13) in Eq. (6), and keeping the leading term alone, one ends up with the equation

Note that, since the external field is now a normally incident plane wave, ${E}_{x}^{ext}$,${E}_{y}^{ext}$ do not depend on *x*. This simple equation is solved by

Next, consider the equation for *J _{x}*. Upon replacing

*G*

_{0}by Eq. (13) in Eq. (3), keeping the leading order term, and using Eq. (3), one arrives at

This equation is identical to the one describing an array of graphene ribbons in the absence of a magnetic field [34], except for *F* replacing the external field ${E}_{x}^{ext}$. Using perturbation theory, the solution of Eq. (18) can be written as

Here ${\psi}_{n}\left(x\right)$ and *q _{n}* satisfy the eigenvalue equation [34, 35]

Note that ${\mathcal{Y}}_{n}$ in Eq. (22) equals the admittance of a capacitance $2{\epsilon}_{0}/{q}_{n}$ in series with a conductance ${\sigma}_{xx}$.

The constant *F* self-consistently depends on *J _{x}* through

*I*[see Eqs. (16) and (19)]. To calculate

_{x}*I*, an integration is carried out over Eq. (21), and Eq. (19) is used. After basic calculations one obtains

_{x}Substitution in Eq. (19) then yields

In all the aforementioned equations, the graphene ribbons are assumed to be suspended in free space. In a different medium with a relative dielectric constant of *ε _{r}*, all the expressions obtains remain valid provided that

*ε*

_{0}and

*k*

_{0}are replaced by ${\epsilon}_{0}{\epsilon}_{r}$ and ${k}_{0}\sqrt{{\epsilon}_{r}}$, respectively.

These calculations can be readily generalized to the cases of oblique incidence. This generalization can be obtained by the correction of $G\left(x\right)$. Consider an array of graphene ribbons when illuminated by an oblique incidence angle (*θ _{i}*). Equation (13) can be re-expressed as

After some straightforward mathematical manipulations it can be shown that Eqs. (21)–(28) are valid with this modification

It should be noted that in deriving Eq. (32) and Eq. (33), we neglected the spatial variations of ${E}_{x/y}^{ext}$ along the *x*-direction.

#### 2.4. Scattering parameters

The distribution of surface current on graphene ribbons is next used to calculate the reflection and transmission coefficients of the structure. To that end, the electric field in the region $z<0$ is expressed as the Rayleigh expansion

Similarly, the electric field in the region $z>0$ is written as

In these relationships

In the sub-wavelength regime where ${k}_{0}D\ll 2\pi $ only the zero-order mode (*n* = 0) is propagating. All other modes are evanescent. For sub-wavelength ribbon arrays the reflection and transmission coefficient are, therefore, defined using the zero-order mode amplitudes alone, according to the relationships

In order to compute the reflection and transmission matrices, the usual boundary conditions on electric and magnetic fields is applied at *z* = 0, and an integration is carried out over the resulting equations on a single period, i.e., from *x* = 0 to $x=D$. It then follows that

The next step is to express ${I}_{x},{I}_{y}$ in terms of the incident field. To that end, we carry out an integration over Eq. (15) from $x=-w/2$ to $x=w/2$. This results in

For simplification, in all the equations of this subsection, we assume that graphene is suspended in free space. These equations can be generalized to the case where the ribbon is sandwiched between two homogenous media with permittivities *ε*_{1} and *ε*_{2}, respectively as shown in Fig. 1. However, this method has complicated calculations. In this case, we employ another simple method to solve the integral equations of Eq. (14) and Eq. (18). In this method we omit the term $-j\omega {\mu}_{0}{\int}_{-w/2}^{w/2}G\left(x-{x}^{\prime}\right){J}_{x/y}\left({x}^{\prime}\right)d{x}^{\prime}$ but account for its effect by including the reflected fields in the external fields [34, 36]. ${E}_{x}^{ext}$ and ${E}_{y}^{ext}$ now not only include the true external field but also the field reflected by the ribbon array. For example, for an incident TM wave with unity amplitude, the total amplitude of the electric field must be replaced by $1+{R}_{xx}$ and *R _{yx}* in the

*x*and

*y*directions, respectively.

## 3. Result, discussion and application

#### 3.1. Numerical results

In this section we validate and verify the accuracy of the proposed method through some examples. Consider a periodic array of graphene ribbons in free space with $D=4\mu m,w=2\mu m,{E}_{f}=0.5eV,\tau =1ps,{B}_{0}=10T$ for the first example. For a normally incident plane wave, the magnitude (absolute value) of reflection coefficients Eqs. (43), (44) and (45) is plotted as function of frequency in Fig. 2. For the sake of comparison the results obtained using effective medium theory [29] (dotted blue), and finite integration technique (FIT) (dashed black) are also plotted. The FIT results are obtained using CST Microwave Studio 2018. Excellent agreement can be observed between the FIT results and our analytical method. By contrast, the effective medium approach [29] can not reproduce the electromagnetic response of magnetically-biased graphene ribbons.

The resonance features observed in Fig. 2 can be explained as follows. Equations (21) and (22) imply that when ${q}_{n}{\sigma}_{xx}+2j\omega {\epsilon}_{0}$ becomes very small for some *n*, ${\mathcal{Y}}_{n}$ becomes largeand a resonance occurs. Near the corresponding resonance frequency, the distribution of ${J}_{x}\left(x\right)$ is predominantly given by ${\psi}_{n}\left(x\right)$. Apart from constant terms, the distribution of ${J}_{y}\left(x\right)$ is identical to ${J}_{x}\left(x\right)$ as can be seen from Eq. (15). Figure 3 shows the magnitude of ${J}_{x},{J}_{y}$ across a graphene ribbon at the first (*n* = 1) and second (*n* = 3) resonances, that occur at 9.78- and 19.13 THz, respectively. It must be noted that ${\psi}_{n}\left(x\right)$ is an odd function of *x* with respect to *x* = 0 (center of a ribbon) for even values of *n*. As a result *S _{n}* = 0 for even values of

*n*[see Eq. (23)] so that only odd values of

*n*enter the summation Eq. (21). The resonances observed are thus due to odd modes $n=1,3,\cdots $. Near these resonance frequencies, peaks in reflection coefficients are observed. Note, however, that

*R*does not show any resonance features in the absence of an external magnetic bias, i.e., when ${\sigma}_{xy}=0$. This can be seen from Eq. (45).

_{yy}As a second example, consider a periodic array of graphene ribbons that are illuminated by an oblique incident wave. The periodic graphene ribbon structure has $D=10\mu m,w=7.5\mu m$, ${E}_{f}=0.75eV$, $\tau =1ps$, ${B}_{0}=6T$. The result of the magnitude and phase and of *R _{yx}* vs frequency is presented in Fig. 4. For brevity, only

*R*has been plotted here. A very good agreement is seen between the analytical method and the FIT. It can be observed from Fig. 4 that oblique incidence does not change the resonance frequency, and this point is clear from Eqs. (21),(22),(32) and (33). Oblique incidence does not change the eigenfunction and eigenvalue of the magnetically-biased graphene ribbons, and as previously mentioned each resonance frequency corresponds to an eigenfunction and eigenvalue. It can be further seen from the figure that if the angle of incidence increases, the magnitude of

_{yx}*R*can increase Eq. (44).

_{yx}With regards to the limitations of the presented theory, it should be noted that our analysis assumed operation in the sub-wavelength regime where ${k}_{0}D\ll 1$. As shown in Fig. 2, when the frequency increases, the accuracy of the proposed method decreases. For the position of the first resonance, we have a relative error of $0.5\%$ while the relative error for the fourth resonance is $2.71\%$. Moreover, solution Eq. (21) is based on perturbation theory in which the interaction between neighboring ribbons is assumed not to affect the eigenvalue Eq. (24) [34]. Increasing the fill factor ($w/D$), will lead to a smaller gap and a stronger interaction between neighboring ribbons. One would, therefore, expect the accuracy of our method to decrease. The frequency of the first resonance is plotted as function of fill factor in Fig. 5 for $D=10\mu m,{E}_{f}=0.3eV,\tau =1ps,{B}_{0}=7.5T$. It is observed that by increase in the fill factor, the relative error in predicting the resonance frequency is slightly increased. Finally, derivation of Eq. (9) was based on Eq. (8). According to Eqs. (2) and (7), when *E _{f}* is large, Eq. (8) cannot be satisfied. If the Fermi level energy increases, the accuracy of the proposed method can slightly decrease. The effect of this approximation, as well as perturbation theory, is small compared to the sub-wavelength approximation. The latter, however, does not pose a strict limitation, since anisotropic graphene ribbons have almost always been utilized in the sub-wavelength regime [21–23, 29]. Consequently, the proposed analysis is much more accurate and affordable than its counterparts already reported in the literature [23, 28, 29].

#### 3.2. Design of a Faraday rotator

Faraday rotation is one of the most attractive phenomena in optics, with many applications in optical diodes, sensing, magnetic microscopy, optical communications, data storage, optical isolators, and phase modulators [20, 23]. For the array of graphene ribbons under magnetic bias, Faraday rotation angle can be directly computed as follows [23]:

*T*and

_{xx}*T*can be calculated from Eqs. (43) and (44). Large Faraday rotation angles are observed at resonance. To achieve a giant Faraday rotation at 10 THz, the structure parameters are designed as $D=4.5\mu m,w=2.7\mu m,{E}_{f}=0.8eV,\tau =2ps,{B}_{0}=7T$. The results are shown in Fig. 6. A good agreement between the results of our method and numerical FIT is observed. Also, we compare our results with the semi-analytical results [23]. It is clear that our method yields much more accurate results for the Faraday rotation angle of the magnetically-biased graphene ribbons. This is because the analysis in [23] is based on the low-relaxation-time assumption. Besides, it assumes that period

_{yx}*D*is much smaller than the free-space wavelength. Neither condition is satisfied here. Finally, we investigate the effect of magnetic bias on performance of the Faraday rotator in Fig. 6(b). It can be seen that increasing the magnetic bias leads to larger Faraday rotation angles. Also, giant Faraday rotation occurs at higher frequencies.

## 4. Conclusion

Based on integral equations governing the field and the surface current density on magnetically-biased graphene ribbons, this paper proposed a analytical, fast, and accurate method for the analysis of such structures. Unlike previous studies, the proposed method does not suffer from extremely slow convergence rates of numerical methods, and more importantly is able to predict spectral resonance features, not possible with previous analytic solutions. The expressions were derived under a quasi-static approximation and extended to periodic arrays of graphene ribbons using perturbation theory. By Rayleigh expansion, we determined the reflection/transmission coefficient using a rigorous method. The proposed method was validated against commercial software equipped with finite integration technique (FIT), showing excellent agreement between the two solutions. In addition, a Faraday rotator was designed and investigated as an important application. The present study suggests a new method for the analysis of magnetically-biased graphene ribbons, which is more precise, reliable, and affordable than the methods already reported in the literature. The proposed method might also be useful for the analysis and development of different devices such as polarizer, isolators, and so on.

## 5. Appendix

Consider the Green function Eq. (12). The Hankel function is written as

*θ*plane. Since $\left|x\right|<D$, one has:

Equation (48) can be rewritten as

The integrand in Eq. (49) is devoid of singularities arising at $\theta =0,\pi $ so that the integration contour is deformed into the one shown on page 916 of [37]. The resulting integral is written as

Since ${k}_{0}D$ is small (and so is ${k}_{0}x$) the integrand in the first integral can be expanded in powers of ${k}_{0}D$,${k}_{0}\left|x\right|$. The result can be easily shown to be of the order ${k}_{0}D$. To compute the second integral it is rewritten as

## References

**1. **A. K. Geim and K. S. Novoselov, “The rise of graphene,” in *Nanoscience and Technology: A Collection of Reviews from Nature Journals*, (World Scientific, 2010), pp. 11–19.

**2. **K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric field effect in atomically thin carbon films,” science **306**, 666–669 (2004). [CrossRef] [PubMed]

**3. **A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, “Superior thermal conductivity of single-layer graphene,” Nano letters **8**, 902–907 (2008). [CrossRef] [PubMed]

**4. **Q. Bao and K. P. Loh, “Graphene photonics, plasmonics, and broadband optoelectronic devices,” ACS nano **6**, 3677–3694 (2012). [CrossRef] [PubMed]

**5. **M. Rahmanzadeh, H. Rajabalipanah, and A. Abdolali, “Multilayer graphene-based metasurfaces: robust design method for extremely broadband, wide-angle, and polarization-insensitive terahertz absorbers,” Appl. optics **57**, 959–968 (2018). [CrossRef]

**6. **A. Momeni, K. Rouhi, H. Rajabalipanah, and A. Abdolali, “An information theory-inspired strategy for design of re-programmable encrypted graphene-based coding metasurfaces at terahertz frequencies,” Sci. reports **8**, 6200(2018). [CrossRef]

**7. **P.-Y. Chen, J. Soric, Y. R. Padooru, H. M. Bernety, A. B. Yakovlev, and A. Alù, “Nanostructured graphene metasurface for tunable terahertz cloaking,” New J. Phys. **15**, 123029 (2013). [CrossRef]

**8. **M. Rahmanzadeh, H. Rajabalipanah, and A. Abdolali, “Analytical investigation of ultrabroadband plasma–graphene radar absorbing structures,” IEEE Transactions on Plasma Sci. **45**, 945–954 (2017). [CrossRef]

**9. **J. T. Kim and C.-G. Choi, “Graphene-based polymer waveguide polarizer,” Opt. express **20**, 3556–3562 (2012). [CrossRef] [PubMed]

**10. **S. AbdollahRamezani, K. Arik, A. Khavasi, and Z. Kavehvash, “Analog computing using graphene-based metalines,” Opt. letters **40**, 5239–5242 (2015). [CrossRef]

**11. **R. Filter, M. Farhat, M. Steglich, R. Alaee, C. Rockstuhl, and F. Lederer, “Tunable graphene antennas for selective enhancement of thz-emission,” Opt. express **21**, 3737–3745 (2013). [CrossRef] [PubMed]

**12. **K. Rouhi, H. Rajabalipanah, and A. Abdolali, “Multi-bit graphene-based bias-encoded metasurfaces for real-time terahertz wavefront shaping: From controllable orbital angular momentum generation toward arbitrary beam tailoring,” Carbon **149**, 125–138 (2019). [CrossRef]

**13. **B. Sensale-Rodriguez, S. Rafique, R. Yan, M. Zhu, V. Protasenko, D. Jena, L. Liu, and H. G. Xing, “Terahertz imaging employing graphene modulator arrays,” Opt. express **21**, 2324–2330 (2013). [CrossRef] [PubMed]

**14. **C.-C. Lee, S. Suzuki, W. Xie, and T. Schibli, “Broadband graphene electro-optic modulators with sub-wavelength thickness,” Opt. express **20**, 5264–5269 (2012). [CrossRef] [PubMed]

**15. **G. W. Hanson, “Dyadic green’s functions for an anisotropic, non-local model of biased graphene,” IEEE Transactions on Antennas Propag. **56**, 747–757 (2008). [CrossRef]

**16. **M. Tamagnone, C. Moldovan, J.-M. Poumirol, A. B. Kuzmenko, A. M. Ionescu, J. R. Mosig, and J. Perruisseau-Carrier, “Near optimal graphene terahertz non-reciprocal isolator,” Nat. communications **7**, 11216 (2016). [CrossRef]

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

**18. **N. Chamanara, D. Sounas, and C. Caloz, “Non-reciprocal magnetoplasmon graphene coupler,” Opt. express **21**, 11248–11256 (2013). [CrossRef] [PubMed]

**19. **M. Tamagnone, A. Fallahi, J. R. Mosig, and J. Perruisseau-Carrier, “Fundamental limits and near-optimal design of graphene modulators and non-reciprocal devices,” Nat. photonics **8**, 556 (2014). [CrossRef]

**20. **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**, 48 (2011). [CrossRef]

**21. **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. communications **8**, 14626 (2017). [CrossRef]

**22. **A. Fallahi and J. Perruisseau-Carrier, “Manipulation of giant faraday rotation in graphene metasurfaces,” Appl. Phys. Lett. **101**, 231605 (2012). [CrossRef]

**23. **M. Tymchenko, A. Y. Nikitin, and L. Martín-Moreno, “Faraday rotation due to excitation of magnetoplasmons in graphene microribbons,” ACS nano **7**, 9780–9787 (2013). [CrossRef] [PubMed]

**24. **D. L. Sounas and C. Caloz, “Electromagnetic nonreciprocity and gyrotropy of graphene,” Appl. Phys. Lett. **98**, 021911 (2011). [CrossRef]

**25. **X.-H. Wang, W.-Y. Yin, and Z. Chen, “Matrix exponential fdtd modeling of magnetized graphene sheet,” IEEE Antennas Wirel. Propag. Lett. **12**, 1129–1132 (2013). [CrossRef]

**26. **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 Transactions on Antennas Propag. **63**, 4458–4467 (2015). [CrossRef]

**27. **P. K. Khoozani, M. Maddahali, M. Shahabadi, and A. Bakhtafrouz, “Analysis of magnetically biased graphene-based periodic structures using a transmission-line formulation,” JOSA B **33**, 2566–2576 (2016). [CrossRef]

**28. **M. Feizi, V. Nayyeri, and O. M. Ramahi, “Modeling magnetized graphene in the finite-difference time-domain method using an anisotropic surface boundary condition,” IEEE Transactions on Antennas Propag. **66**, 233–241 (2018). [CrossRef]

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

**30. **A. Ferreira, J. Viana-Gomes, Y. V. Bludov, V. Pereira, N. Peres, and A. C. 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**, 235410 (2011). [CrossRef]

**31. **L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. Liang, A. Zettl, Y. R. Shen, *et al.*, “Graphene
plasmonics for tunable terahertz
metamaterials,” Nat.
Nanotechnology **6**, 630
(2011). [CrossRef]

**32. **C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, *et al.*, “Boron nitride
substrates for high-quality graphene
electronics,” Nat.
Nanotechnology **5**, 722
(2010). [CrossRef]

**33. **C.-T. Tai, *Dyadic Green Functions in
Electromagnetic Theory* (Institute of
Electrical & Electronics Engineers (IEEE),
1994).

**34. **A. Khavasi and B. Rejaei, “Analytical modeling of graphene ribbons as optical circuit elements,” IEEE J. Of Quantum Electron. **50**, 397–403 (2014). [CrossRef]

**35. **M. Rahmanzadeh, A. Abdolali, A. Khavasi, and H. Rajabalipanah, “Adopting image theorem for rigorous analysis of a perfect electric conductor–backed array of graphene ribbons,” JOSA B **35**, 1836–1844 (2018). [CrossRef]

**36. **S. Barzegar-Parizi, B. Rejaei, and A. Khavasi, “Analytical circuit model for periodic arrays of graphene disks,” IEEE J. Quantum Electron. **51**, 1–7 (2015). [CrossRef]

**37. **I. S. Gradshteyn and I. M. Ryzhik, *Table of Integrals, Series, and
Products* (Academic press,
2014).