An Analytical and Rigorous Method for Analysis of an Array of Magnetically-Biased Graphene Ribbons

TL;DR

This paper introduces an accurate analytical method for analyzing magnetically-biased graphene ribbon arrays, enabling precise prediction of their electromagnetic behavior and Faraday rotation effects, surpassing existing semi-analytical approaches.

Contribution

It presents a novel, fast, and reliable integral-equation-based analytical technique for periodic graphene ribbon arrays under magnetic bias, including a closed-form solution for Faraday rotation.

Findings

Excellent agreement with full-wave simulations

Predicts resonant spectral effects not captured by previous methods

Provides a closed-form expression for Faraday rotation

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 normally incident plane waves, 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.

Equations

$${\overline{\overline{\sigma}}_{s}}(\omega,{E_{f}},\tau,T,{B_{0}})=\left[\begin{array}[]{*{20}{c}}{{\sigma_{xx}}}&{{\sigma_{xy}}}\\ {{\sigma_{yx}}}&{{\sigma_{yy}}}\end{array}\right]$$

$${\sigma }_{xx}={\sigma }_{yy}=\frac{e^2{E}_f\tau }{\pi {\hslash }^2}\frac{1+j\omega \tau }{{({\omega }_c\tau )}^2+{(1+j\omega \tau )}^2}$$

$${\sigma }_{xy}=-{\sigma }_{yx}=\frac{e^2{E}_f\tau }{\pi {\hslash }^2}\frac{{\omega }_c\tau }{{({\omega }_c\tau )}^2+{(1+j\omega \tau )}^2}$$

$${\overline{\overline{{\sigma_{s}}}}^{-1}}\left[\begin{array}[]{*{20}{c}}{{J_{x}}(x)}\\ {{J_{y}}(x)}\end{array}\right]=\left[\begin{array}[]{*{20}{c}}{E_{x}^{ext}(x)+E_{x}^{T}(x)}\\ {E_{y}^{ext}(x)+E_{y}^{T}(x)}\end{array}\right]$$

$$\begin{array}{rl}& E_x^T=\frac{1}{j\omega {\epsilon }_0}\frac{d}{dx}{\int }_{-w/2}^{w/2}{G}_0(x-x^{\prime })\frac{d{J}_x(x^{\prime })}{d{x}^{\prime }}d{x}^{\prime }-\\ & j\omega {\mu }_0\underset{-w/2}{\overset{w/2}{\int }}{G}_0(x-x^{\prime })J_x(x^{\prime })d{x}^{\prime }\end{array}$$

$$E_y^T=-j\omega {\mu }_0\underset{-w/2}{\overset{w/2}{\int }}{G}_0(x-x^{\prime })J_y(x^{\prime })d{x}^{\prime }$$

$$G_0(x-x^{\prime })=\frac{1}{4j}{H}_0^{(2)}(k_0|x-x^{\prime }|)$$

$$G_0(x-x^{\prime })\simeq -\frac{1}{2\pi }\ln \left(k_0|x-x^{\prime }|\right)$$

$$\begin{array}{rl}{J}_y(x)+& j\omega {\mu }_0{\sigma }_0{\int }_{-w/2}^{w/2}{G}_0(x-x^{\prime })J_y(x^{\prime })d{x}^{\prime }\\ & ={\sigma }_0\left[E_y^{ext}(x)-{\sigma }_a^{-1}{J}_x(x)\right]\end{array}$$

$${\sigma }_0^{-1}=\frac{{\sigma }_{xx}}{{\sigma }_{xx}^2+{\sigma }_{xy}^2},{\sigma }_a^{-1}=\frac{{\sigma }_{xy}}{{\sigma }_{xx}^2+{\sigma }_{xy}^2}$$

$$\left|\frac{j\omega {\mu }_{0}w{\sigma }_0}{4}\right|\ll 1$$

$$\begin{array}{r}{J}_y(x)={\sigma }_0\left[E_y^{ext}(x)-{\sigma }_a^{-1}{J}_x(x)\right]\end{array}$$

$$\frac{j\omega {\mu }_{0}w{\sigma }_0}{4}=(k_{0}w)\frac{k_0}{k_P}$$

$$\begin{array}{rl}\frac{J_x(x)}{{\sigma }_{xx}}-& \frac{1}{j\omega {\epsilon }_0}\frac{d}{dx}{\int }_{-w/2}^{w/2}{G}_0(x-x^{\prime })\frac{d{J}_x(x^{\prime })}{d{x}^{\prime }}d{x}^{\prime }\\ =& E_x^{ext}(x)+\frac{{\sigma }_{xy}}{{\sigma }_{xx}}{E}_y^{ext}(x)\end{array}$$

$$G\left(x-x^{\prime }\right)=\frac{1}{4j}\sum _{l=-\infty }^{\infty }{H}_0^{(2)}\left(k_0\left|x-x^{\prime }-lD\right|\right)$$

$$\begin{array}{r}G(x)=\frac{1}{2j{k}_{0}D}+\frac{1}{2\pi }\underset{0}{\overset{\infty }{\int }}\left[\frac{\cosh \left(\alpha u\right)}{\sinh u}-\frac{1}{u}\right]\frac{du}{u}+O\left(k_{0}D\right)\end{array}$$

$$J_y(x)={\sigma }_0\left[E_y^{ext}(x)-{\sigma }_a^{-1}{J}_x(x)\right]-\frac{\omega {\mu }_0{\sigma }_0}{2k_{0}D}{\int }_{-w/2}^{w/2}{J}_y(x^{\prime })d{x}^{\prime }$$

$$J_y(x)={\sigma }_0\left[\frac{E_y^{ext}}{1+\gamma }+\frac{\gamma }{1+\gamma }\frac{I_x}{{\sigma }_{a}w}-\frac{J_x(x)}{{\sigma }_a}\right]$$

$$I_x={\int }_{-w/2}^{w/2}{J}_x(x^{\prime })d{x}^{\prime }$$

$$\gamma =\frac{\omega {\mu }_0{\sigma }_{0}w}{2k_{0}D}=\frac{{\eta }_0{\sigma }_{0}w}{2D}$$

$$\begin{array}{rl}& \frac{J_x(x)}{{\sigma }_{xx}}=F+\\ & \frac{1}{j\omega {\epsilon }_0}\sum _{l=-\infty }^{\infty }\frac{d}{dx}{\int }_{-w/2}^{w/2}{G}_0(x-x^{\prime }-lD)\frac{d{J}_x(x^{\prime })}{d{x}^{\prime }}d{x}^{\prime }\end{array}$$

$$\begin{gathered} F=E_x^{ext}+\frac{{\sigma }_{xy}}{{\sigma }_{xx}}\left[\frac{E_y^{ext}}{1+\gamma }+\frac{\gamma }{1+\gamma }\frac{I_x}{{\sigma }_{a}w}\right] \\ -\frac{\omega {\mu }_0{I}_x}{2k_{0}D} \end{gathered}$$

$$\frac{dG(x)}{dx}=\sum _{l=-\infty }^{\infty }\frac{d{G}_0(x-lD)}{dx}+O\left(k_{0}D\right)$$

$$J_x(x)=F\sum _{n=1}^{\infty }{\mathcal{Y}}_n{S}_n{\psi }_n(x)$$

$${\mathcal{Y}}_n=\frac{{\sigma }_{xx}(2j\omega {\epsilon }_0/q_n)}{{\sigma }_{xx}+2j\omega {\epsilon }_0/q_n}$$

$$S_n=\underset{-w/2}{\overset{w/2}{\int }}{\psi }_n(x^{\prime })d{x}^{\prime }$$

$$\frac{1}{\pi }{\int }_{-w/2}^{w/2}\frac{1}{x-x^{\prime }}\frac{d{\psi }_n(x^{\prime })}{d{x}^{\prime }}d{x}^{\prime }=q_n{\psi }_n(x)$$

$${\int }_{-w/2}^{w/2}{\psi }_n^2(x)dx=1$$

$$I_x=\frac{Y}{1+\frac{1}{2}{\zeta }_{0}Y}D{F}_0$$

$$Y=\frac{1}{D}\sum _{m=1}^{\infty }{\mathcal{Y}}_m{S}_m^2$$

$$F_0=E_x^{ext}+\frac{{\sigma }_{xy}}{{\sigma }_{xx}}\frac{E_y^{ext}}{1+\gamma }$$

$${\zeta }_0={\eta }_0\left[1-\frac{{\sigma }_{xy}^2}{{\sigma }_{xx}^2(1+\gamma )}\right]$$

$$F=\frac{F_0}{1+\frac{1}{2}{\zeta }_{0}Y}$$

$$\begin{array}{rl}\left(E_x^{ext}\hat{\mathbf{x}}+E_y^{ext}\hat{\mathbf{y}}\right)& e^{-j{k}_{0}z}\\ +\sum _{n=-\infty }^{\infty }& \left(E_{x,n}^r\hat{\mathbf{x}}+E_{y,n}^r\hat{\mathbf{y}}\right)e^{-j{k}_{x,n}x+j{k}_{z,n}z}\end{array}$$

Full text

I 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, 6, 7, 8, 9, 10, 11, 12] Under an external magnetic bias, the surface conductivity of graphene becomes a tensor, and it will also have non-reciprocal properties [13]. These extraordinary properties make graphene a promising material for many novel applications, such as isolators [14], nonreciprocal couplers [15], optical modulators [16], and Faraday rotators [17], 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 [18, 19, 20]. For example, in [20] 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.

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. Saunas and Caloz studied the transmission and reflection properties of magnetically-biased graphene sheet using the general anisotropic conductivity tensor [21]. 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 [20], the finite difference time-domain method [22], and discontinuous Galerkin time domain (DGTD) [23] 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 [24, 25]; 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 [26]; 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 Sec. II 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 Sec. III, the reflection/transmission coefficient of the structure is calculated and the proposed method is validated against full-wave simulations. 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.

II SURFACE CURRENTS ON A PERIODIC 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_{y}$ will be induced on the ribbons. Since there is no variation in $y$-direction, these induced currents are functions of $x$ only (i.e. $\partial_{y}{J_{x,y}}=0$).

II-A 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}\hat{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 [13]

[TABLE]

where $E_{f}$ is the Fermi level energy and $\omega$ is the angular frequency, $\tau$ is the relaxation time, and $T$ is the Kelvin temperature ($T=300K$). For highly doped graphene, i.e. $E_{f}>>\hbar\omega$ and $E_{f}>>k_{B}T$ ($k_{B}$ = Boltzmann constant, $\hbar$ = reduced Planck’s constant), (1) can be simplified by a Drude-like model as [27]

[TABLE]

[TABLE]

where $\omega_{c}={e}{B_{0}}{f^{2}_{F}}/{E_{f}}$ is the cyclotron frequency with $v_{f}=10^{6}m/s$ denoting the Fermi velocity. The relaxation time can be described by $\tau=\mu\hbar\sqrt{{n_{s}}\pi}/e{v_{F}}$ where $\mu$ is the carrier mobility and ${n_{s}}=E_{F}^{2}/\pi{\hbar^{2}}v_{F}^{2}$ is the carrier density [21]. Carrier mobility is assumed to be within the range of $0.1$ to $6$ $m^{2}/V$, which is practical and achievable for graphene with different substrates at room temperature [28, 29]. It should also be noted that the time harmonic of the form $\exp(j{\omega}t)$ is assumed throughout this paper.

II-B 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 [30]

[TABLE]

represent the electric field generated by the surface current $J_{x}$ and $J_{y}$. In these equations

[TABLE]

designates the Green’s function of the 2D Helmholtz equation. Here $H_{0}^{(2)}$ is the Hankel function of the second kind and $k_{0}=\omega(\varepsilon_{0}\mu_{0})^{1/2}$ is the vacuum wave number.

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

[TABLE]

From (3a) and (3c) one obtains

[TABLE]

where

[TABLE]

It must be noted that $\sigma_{0}$ is the graphene conductivity in the absence of the external magnetic field.

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

[TABLE]

the second term on the left hand side of (6) becomes negligible and the solution is trivially given by

[TABLE]

To see why (8) is satisfied, note that

[TABLE]

where $k_{P}=-2j\omega\varepsilon_{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, (8) is satisfied.

Let us now return to (3a) and combine it with (3b) to arrive at an equation for $J_{x}(x)$. However, one may simplify the resulting equation by neglecting the second term on the right hand side of (3b). This term is small compared to the first term for a narrow ribbon. The integral equation governing $J_{x}$ thus becomes

[TABLE]

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

II-C 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 (3a) 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, (3c) and (3c) remain provided that the Green’s function (4) is replaced by

[TABLE]

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 (12), leads to an important modification of the approximation (5). For sub-wavelength arrays where $k_{0}D\ll 1$, a lengthy calculation presented in the appendix shows that

[TABLE]

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 (6) cannot be neglected anymore. By replacing $G_{0}$ with (13) in (6), and keeping the leading term alone, one ends up with the equation

[TABLE]

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

[TABLE]

where

[TABLE]

where $\eta_{0}=(\mu_{0}/\varepsilon_{0})^{1/2}$.

Next, consider the equation for $J_{x}$. Upon replacing $G_{0}$ by (13) in (3c), keeping the leading order term, and using (3a), one arrives at

[TABLE]

where

[TABLE]

and we have used

[TABLE]

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

[TABLE]

Here $\psi_{n}(x)$ and $q_{n}$ satisfy the eigenvalue equation[31, 32]

[TABLE]

where $\psi_{n}(x)$ is subject to the normalization condition

[TABLE]

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

The constant $F$ self-consistently depends on $J_{x}$ through $I_{x}$ [see (16),(19)] . To calculate $I_{x}$, an integration is carried out over (21), and (19) is used. After basic calculations one obtains

[TABLE]

Substitution in (19) then yields

[TABLE]

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 $\varepsilon_{r}$, all the expressions obtains remain valid provided that $\varepsilon_{0}$ and $k_{0}$ are replaced by $\varepsilon_{0}\varepsilon_{r}$ and $k_{0}\sqrt{\varepsilon_{r}}$, respectively.

II-D 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

[TABLE]

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

[TABLE]

In these relationships

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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 (15) from $x=-w/2$ to $x=w/2$. This results in

[TABLE]

where we have used (7). Combination of (26),(28), (51), (48)- (40) then yields

[TABLE]

For simplification, in all the equations of this subsection, we assume that graphene is suspended in free space. These equations can be readily generalized to the case where the ribbon is sandwiched between two homogenous media with permittivities $\varepsilon_{1}$ and $\varepsilon_{2}$, respectively as shown in Fig.1.

III RESULT, DISCUSSION AND APPLICATION

III-A *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 a normally incident plane wave, the magnitude of reflection coefficients (52), (53) and (54) is plotted as function of frequency in Fig.2. For the sake of comparison the results obtained using effective medium theory [26] (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 [26] 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),(22) imply that when $q_{n}\sigma_{xx}+2j\omega\varepsilon_{0}$ becomes very small for some $n$, $\mathcal{Y}_{n}$ becomes large and a resonance occurs. Near the corresponding resonance frequency, the distribution of $J_{x}(x)$ is predominantly given by $\psi_{n}(x)$. Apart from constant terms, the distribution of $J_{y}(x)$ is identical to $J_{x}(x)$ as can be seen from (15). Fig. 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}(x)$ 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 (23)] so that only odd values of $n$ enter the summation (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_{yy}$ 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 ((54)).

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 (21) is based on perturbation theory in which the interaction between neighboring ribbons is assumed not to affect the eigenvalue equation (24) [31]. 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.4 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 (9) was based on (8). According to (2) and (7), when $E_{f}$ is large, (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 [18, 19, 20, 26]. Consequently, the proposed analysis is much more accurate and affordable than its counterparts already reported in the literature [20, 25, 26].

III-B *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 [17, 20]. For the array of graphene ribbons under magnetic bias, Faraday rotation angle can be directly computed as follows[20]:

[TABLE]

where $T_{xx}$ and $T_{yx}$ can be calculated from (52) and (53). 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.5. A good agreement between the results of our method and numerical FIT is observed. Also, we compare our results with the semi-analytical results [20]. 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 [20] is based on the low-relaxation-time assumption. Besides, it assumes that period $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.5b. It can be seen that increasing the magnetic bias leads to larger Faraday rotation angles. Also, giant Faraday rotation occurs at higher frequencies.

IV 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.

Appendix A

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

[TABLE]

were the integration is carried out over the path shown on page 916 of [33] in the complex $\theta$ plane. Since $\left|x\right|<D$, one has:

[TABLE]

where it is assumed that the integration contour does not pass through points where $\sin(\theta)=0$.

(A.2) can be rewritten as

[TABLE]

where we have used

[TABLE]

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

[TABLE]

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}|x|$. The result can be easily shown to be of the order $k_{0}D$. To compute the second integral it is rewritten as

[TABLE]

where $\alpha=1-2|x|/D$ and $u=(k_{0}D/2)\sinh t$. Inserting this result in (A.3) leads to (13). Take note that

[TABLE]

which is identical to (20).