Theory for polaritons in graphene photonic crystals in an applied magnetic field

One-dimensional (1D) graphene photonic crystals (PCs) are multilayered systems where one of the constituent materials is graphene, the well-known one-atom-thick metallic film discovered in 2004 [1]. These graphene PCs form an emergent research area of intense interest [2–8]. Systems such as PCs [9] or quasicrystals [10] that incorporate the atomically-thin graphene material have been studied because they exhibit novel optical properties and may be valuable for modern applications in the field of tunable plasmonic and graphene electronics at far-infrared frequencies, for example. Also, these structures can support highly confined plasmon polaritons with very low loss [11] due to graphene's unique electronic dispersion with gapless features near the so-called Dirac points in the Brillouin zone [12].

Polaritons are hybrid (or mixed) modes that occur when a photon strongly couples to an excitation in the material [13, 14]. Examples include exciton-polaritons, when a photon couples to an electron–hole pair in a semiconductor, or plasmon-polaritons, when the coupling is with the electron plasma in a metal or doped semiconductor. On the other hand, graphene constitutes a very special kind of 2D electron gas system. For example, the capability to change its carrier concentration using an applied voltage-gate opens the way to precise control of the frequency range for excitations like surface plasmon-polaritons (SPP) [15]. Also long propagation lengths can be achieved [16] compared with conventional SPPs.

On the other hand, it is known that when an external magnetic field is applied perpendicular to a 2D electron gas [16] a mix between cyclotron excitations and the electronic oscillations (or plasmons) occurs in that system, leading to magnetoplasmon-polaritons [17]. An effect of the magnetic field is to cause an anisotropy in the optical conductivity tensor and strong resonant absorptions in the dielectric functions. Consequently the dispersion properties and frequencies of the polaritons become very sensitive to the applied magnetic field, especially in systems that incorporate graphene sheets [18–20]. For example, Ferreira et al [19] reported an exotic polariton mode with weak damping in a single-layer graphene film with an applied static magnetic field; they termed it a quasi-transverse-electric mode because of the resemblance to conventional transverse electric (TE) modes. Also Melo [21] studied the excitation of surface polaritons in magnetically biased graphene at the interface separating two dielectric media. The focus was on the influence of the magnetostatic field for the surface polaritons in graphene in the 500 GHz to 5 THz frequency range, where the author calculated the reflectance in an attenuated total reflection (ATR) set up [14] and associated the reflectance dip with the excitation of the surface polaritons. Michalski [22] used the modal transmission line theory to calculate the reflectivity and transmissivity rate in a multilayered system with multiple anisotropic conductive graphene sheets at the interfaces. He noted that the analysis of such structures is rendered non-trivial, because of the coupling between the transverse magnetic (TM) and TE polarization field caused by the graphene sheets. It was concluded that a TE plane wave should also excite SPP modes, due to this coupling of the TE and TM waves, and this coupling becomes enhanced with increased magnetic field bias.

In this work we emphasize the role (and its interesting consequences) of an applied magnetic field on the surface and bulk plasmon-polariton modes in a photonic crystal with embedded graphene sheets at the interfaces, henceforth referred to as a graphene photonic crystal (GPC). The plasmon-polaritons are studied here by using Maxwell's equations in a macroscopic formalism, where graphene is characterized by a frequency-dependent optical conductivity tensor. In section 2 we present a modification of the plasmon-polariton theory for dielectric photonic crystal with graphene sheets [9] to relate to our present case where we have an external magnetic field applied. Here the graphene sheets are characterized by a modified optical conductivity, which includes the Landau levels' quantization due to the applied magnetic field, and we concentrate initially on the bulk polaritons. Next the localized surface polaritons of a semi-infinite GPC are considered in section 3. Numerical results are then presented in section 4 when the dielectric material is taken to be SiO₂ in the various layers of the GPC. These results illustrate the important role of the applied magnetic field in GPCs in modifying the polariton dispersion curves, especially when we have a magnetic field of order 1 T. Also we report a sensitive dependence on the polariton dispersion relation curves on the electronic scattering rate parameter and on the intensity of the applied magnetic field in the graphene sheets. In the last section we present the conclusions of this work.

We consider a photonic crystal (PC) consisting of alternating unit cells AB juxtaposed to form the periodic ABABAB... array, with graphene sheets at each of the interfaces and a static magnetic field $\vec {B}_0$ applied in the direction perpendicular to the layers (see figure 1). The building blocks A and B are layers of a dielectric material with thicknesses a and b, respectively.

Zoom In Zoom Out Reset image size

2.1. Field-dependent conductivity of graphene

When a static magnetic field B₀ is applied perpendicular to the graphene sheets, it is known that the conductivity becomes a tensor with the components in the xy-plane satisfying [23]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \widetilde{\sigma}=\left(\begin{array}{@{}rr@{}} \sigma_{xx} & \sigma_{xy} \nonumber \\ \sigma_{yx} & \sigma_{yy} \end{array}\right). \nonumber \end{align} \tag{ 1 }$$

Here, the general expressions for the diagonal and off-diagonal terms (sometimes called Hall conductivities) can be obtained using the Kubo formula [18, 24] as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq1a} \sigma_{xx}&= \sigma_{yy}= \frac{e^{2}v^{2}_F|eB_0|(\hbar{\omega}+2i\Gamma)}{\pi i} \nonumber \\ &\times \sum_{n=0}^{\infty} \Big\{\frac{[n_{F}(M_{\nu})-n_{F}(M_{\nu+1})]-[n_{F}(-M_{\nu})-n_{F}(-M_{\nu+1})]}{(M_{\nu+1}-M_{\nu})^{3}-(\hbar \omega+2i\Gamma)^{2}(M_{\nu+1}-M_{\nu})} \nonumber \\ & \qquad +\frac{[n_{F}(-M_{\nu})-n_{F}(M_{\nu+1})]-[n_{F}(M_{\nu})-n_{F}(-M_{\nu+1})]}{(M_{\nu+1}+M_{\nu})^{3}-(\hbar \omega-2i\Gamma)^{2}(M_{\nu+1}+M_{\nu})} \Big\},\nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq2a} \sigma_{xy}&= -\sigma_{yx}=-\xi \frac{e^{2}v^{2}_F |eB_0|}{\pi } \nonumber \\ &\times \sum_{n=0}^{\infty}\Big\{[n_{F}(M_{\nu})-n_{F}(M_{\nu+1})]+[n_{F}(-M_{\nu})-n_{F}(-M_{\nu+1})] \Big\} \nonumber \\ & \qquad \times\frac{2[M_{\nu}^{2}+M_{\nu+1}^{2}-(\hbar \omega+2i\Gamma)^{2}]}{[M_{\nu}^{2}+M_{\nu+1}^{2}-(\hbar \omega+2i\Gamma)^{2}]^{2}-4M_{\nu}^{2}M_{\nu+1}^{2}}. \nonumber \end{align} \tag{ 3 }$$

The summations over n are for the energies $M_{\nu}=(2\hbar \nu|eB_{0}|v_{F}^2){}^{1/2}$ of the Landau level, n_F is the Fermi–Dirac distribution function at temperature T, $v_{F}$ is the Fermi velocity of the electrons in graphene, $\Gamma$ is the electron scattering rate, $\xi =$ sgn(B₀) is the sign function of $\vec {B}_0$, and e is the electron charge. Here we have considered the excitonic energy to be equal to zero, i.e. $\Delta=0$ in [24].

When $B_{0} \rightarrow 0$, the sums over n for $\sigma_{xy}$ go to zero and $\sigma_{xx}$ converges very slowly. In this case (not studied here) we could alternatively use the Drude form given in [21, 25]. To illustrate the behaviour of the optical conductivity for graphene for the field values of interest, we consider T  =  300 K, chemical potential $\mu_c= 0.2$ eV for the graphene sheets [26–28], and electron scattering rate $\Gamma =0.11$ meV [21]. Figure 2 shows $\sigma_{xx}$ and $\sigma_{xy}$ plotted versus frequency for the values of the applied magnetic field, B₀  =  0.5, 1.0, and 1.5 T, where the solid (dashed) lines correspond to the real (imaginary) part of $\sigma_{xx}$ (in figure 2(a)) and $\sigma_{xy}$ (in figure 2(b)). The vertical lines (located at the frequencies 0.390 and 0.558 THz) characterize the maximum and minimum value for the imaginary part of $\sigma_{xx}$ (real part of $\sigma_{xy}$), respectively, when B₀  =  1.5 T. These values will turn out to be important later for interpreting the bulk polariton dispersion relations.

Zoom In Zoom Out Reset image size

2.2. Bulk polariton modes

In order to study the bulk polariton modes for an infinite GPC, we need to consider solutions of Maxwell's equations within each dielectric medium, subject to the standard electromagnetic boundary conditions at the interfaces with graphene. As a consequence of the Landau levels leading to a non-diagonal conductivity tensor for the graphene sheets, we need to consider the electric and magnetic field solutions within the constituent media A and B of the lth layer in a general way. Specifically, we do not make a separation into particular polarizations like TM or TE. Instead, in both media we write $\vec {E_j}=(E_{xj},E_{yj},E_{zj})$ for the electric fields, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ex} E_{xj} = \left[A_{xj}^{(l)}e^{ik_{zj}(z-z_l)}+ B_{xj}^{(l)}e^{-ik_{zj}(z-z_l)}\right] e^{i(k_xx-\omega t)}, \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ey} E_{yj} = \left[A_{yj}^{(l)}e^{ik_{zj}(z-z_l)} + B_{yj}^{(l)}e^{-ik_{zj}(z-z_l)}\right]e^{i(k_xx-\omega t)}, \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ez} E_{zj} = -(k_x/ k_j)\left[A_{xj}^{(l)}e^{ik_{zj}(z-z_l)} - B_{xj}^{(l)}e^{-ik_{zj}(z-z_l)} \right]e^{i(k_xx-\omega t)}. \nonumber \end{align} \tag{ 6 }$$

The factors $A_{xj}^{(l)}$, $B_{xj}^{(l)}$, $A_{yj}^{(l)}$, and $B_{yj}^{(l)}$ (with j   =  A or B; $l= 1, 2, 3, \dots$) are the amplitudes for the forward- and backward-travelling waves, respectively, and L  =  a  +  b is the periodicity length. Without loss of generality, we have taken the in-plane wave vector to be along the x direction, and so $\newcommand{\e}{{\rm e}} k_{zj}=(\epsilon_j\omega^ 2/c^ 2-k_x^2){}^{1/2}$ for the (complex) wave vector in the perpendicular direction.

Correspondingly, for the magnetic field components $\vec {H_j}=(H_{xj},H_{yj},H_{zj})$ in medium A (or B), we deduce from $\newcommand{\e}{{\rm e}} \vec {H}=(\epsilon_0 c^ 2/i\omega) \vec {\nabla} \times\vec {E}$ that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Hx} H_{xj}=-\frac{k_{zj}\epsilon_0c^2}{\omega} \left [A_{yj}^{(l)}e^{ik_{zj}(z-z_l)} - B_{yj}^{(l)}e^{-ik_{zj}(z-z_l)}\right] e^{i(k_xx-\omega t)} , \nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Hy} H_{yj} =(\frac{\epsilon_0 \epsilon_j \omega}{k_{zj}}) \left [A_{xj}^{(l)}e^{ik_{zj}(z-z_l)} - B_{xj}^{(l)}e^{-ik_{zj}(z-z_l)}\right] e^{i(k_xx-\omega t)}, \nonumber \end{align} \tag{ 8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Hz} H_{zj}=\frac{k_{x}\epsilon_0c^2}{\omega} \left [A_{yj}^{(l)}e^{ik_{zj}(z-z_l)} - B_{yj}^{(l)}e^{-ik_{zj}(z-z_l)}\right] e^{i(k_xx-\omega t)} . \nonumber \end{align} \tag{ 9 }$$

From the electromagnetic boundary conditions at any A  −  B interface we require

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Bound1} \vec {E}_{yA}=\vec {E}_{yB}, \ \ \ \ \ \vec {E}_{xA}=\vec {E}_{xB}, \nonumber \end{align} \tag{ 10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -({H}_{yB}-{H}_{yA})=\sigma_{xx}{E}_{xB} +\sigma_{xy}{E}_{yB}, \label{Bound2} \nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {H}_{xB}-{H}_{xA}=\sigma_{yx}{E}_{xB} +\sigma_{yy}{E}_{yB}. \label{Bound3} \nonumber \end{align} \tag{ 12 }$$

Then, by considering these condition at the interfaces z  =  nL  +  a within a cell and at $z=(n+1)L$ between cells (see figure 1), we can relate the electromagnetic fields for cell n to those for cell n  +  1. From equation (10) we deduce

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app1} A_{xA}^{(n)}e^{ik_{zA}a}+B_{xA}^{(n)}e^{-ik_{zA}a}= A_{xB}^{(n)}+ B_{xB}^{(n)}, \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{yA}^{(n)}e^{ik_{zA}a}+B_{yA}^{(n)}e^{-ik_{zA}a}= A_{yB}^{(n)}+ B_{yB}^{(n)} .\label{app2} \nonumber \end{align} \tag{ 14 }$$

Next from equations (11) and (12) we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app3} &A_{xA}^{(n)}e^{ik_{zA}a}- B_{xA}^{(n)}e^{-ik_{zA}a}= (\overline{\sigma}_{xx}+\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}})A_{xB}^{(n)}+ (\overline{\sigma}_{xx}-\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}}) B_{xB}^{(n)} \nonumber \\ &\qquad\qquad\qquad\qquad \qquad + \overline{\sigma}_{xy}(A_{yB}^{(n)}+B_{yB}^{(n)}), \nonumber \end{align} \tag{ 15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{app4} &A_{yA}^{(n)}e^{ik_{zA}a}- B_{yA}^{(n)}e^{-ik_{zA}a}= (\overline{\sigma}_{yy}+\frac{k_{zB}}{k_{zA}})A_{yB}^{(n)}+ (\overline{\sigma}_{yy}-\frac{k_{zB}}{k_{zA}})B_{yB}^{(n)} \nonumber \\ &\qquad\qquad\qquad\qquad \qquad + \overline{\sigma}_{yx}(A_{xB}^{(n)}+B_{xB}^{(n)}), \nonumber \end{align} \tag{ 16 }$$

where we define $\newcommand{\e}{{\rm e}} \overline{\sigma}_{xx}={\sigma_{xx}k_{zA}}/{\epsilon_A\epsilon_0\omega}$, $\newcommand{\e}{{\rm e}} \overline{\sigma}_{xy}={\sigma_{xy}k_{zA}}/{\epsilon_A\epsilon_0\omega}$, $\newcommand{\e}{{\rm e}} \overline{\sigma}_{yy}={\sigma_{yy}\omega}/{\epsilon_0 c^2 k_{zA}}$, and $\newcommand{\e}{{\rm e}} \overline{\sigma}_{yx}={\sigma_{yx} \omega}/{\epsilon_0 c^2 k_{zA}}$. The above equations can be conveniently re-expressed in a matrix form as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{T1} M_{A} \left(\begin{array}{@{}c@{}} A_{xA}^{(n)} \nonumber \\ B_{xA}^{(n)} \nonumber \\ A_{yA}^{(n)} \nonumber \\ B_{yA}^{(n)} \end{array}\right)=M_{AB} \left(\begin{array}{@{}c@{}} A_{xB}^{(n)} \nonumber \\ B_{xB}^{(n)} \nonumber \\ A_{yB}^{(n)} \nonumber \\ B_{yB}^{(n)} \end{array}\right). \nonumber \end{align} \tag{ 17 }$$

Here the $4\times4$ propagation matrix M_A is simply given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{matrix1} M_{A}= \left(\begin{array}{@{}cccc@{}} e^{ik_{zA}a} & 0 & 0 & 0 \nonumber \\ 0 & e^{-ik_{zA}a} & 0 & 0 \nonumber \\ 0 & 0 & e^{ik_{zA}a} & 0 \nonumber \\ 0 & 0 & 0 & e^{-ik_{zA}a} \end{array}\right), \nonumber \end{align} \tag{ 18 }$$

while the more complicated $4\times4$ interface matrix due to the graphene sheet is defined in block form as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{M-AB} M_{AB} = \left(\begin{array}{@{}cc@{}} M^{11}_{AB} & M^{12}_{AB} \nonumber \\ M^{21}_{AB} & M^{22}_{AB} \end{array} \right), \nonumber \end{align} \tag{ 19 }$$

where the forms of the $2\times 2$ matrices $M_{AB}^{i'j'}$ (with $i'$ and $j'$ taking the values 1 or 2) are specified in the appendix.

An analogous matrix equation arises when we apply the same boundary conditions at the interfaces where $z=(n+1)L$. In this case it relates the amplitudes $A_{pB}^{(n)}$ and $B_{pB}^{(n)}$ ($p=x,y$) to $A_{pA}^{(n+1)}$ and $B_{pA}^{(n+1)}$, giving

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{T2} M_B \left(\begin{array}{@{}c@{}} A_{xB}^{(n)} \nonumber \\ B_{xB}^{(n)} \nonumber \\ A_{yB}^{(n)} \nonumber \\ B_{yB}^{(n)} \end{array}\right)=M_{BA} \left(\begin{array}{@{}c@{}} A_{xA}^{(n+1)} \nonumber \\ B_{xA}^{(n+1)} \nonumber \\ A_{yA}^{(n+1)} \nonumber \\ B_{yA}^{(n+1)} \end{array}\right). \nonumber \end{align} \tag{ 20 }$$

Here M_B and M_BA are the same as M_A and M_AB, respectively, provided the substitutions $A \rightarrow B$, $B \rightarrow A$, and $a \rightarrow b$ in the matrices (19) and (18) are made. On defining the kets formed by the undetermined set of coefficients as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |A_j^{(n)}\rangle= \left(\begin{array}{@{}c@{}} A_{xj}^{(n)} \nonumber \\ B_{xj}^{(n)} \nonumber \\ A_{yj}^{(n)} \nonumber \\ B_{yj}^{(n)} \end{array}\right), \nonumber \end{align} \tag{ 21 }$$

we find that equations (17) and (20) can be re-expressed as the matrix equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_A|A_A^{(n)}\rangle=M_{AB}|A_B^{(n)} \rangle, \nonumber \end{align} \tag{ 22 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_B|A_B^{(n)} \rangle=M_{BA}|A_A^{(n+1)} \rangle. \nonumber \end{align} \tag{ 23 }$$

It is easy to show that, using these equations, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{rel1} |A_A^{(n+1)}\rangle = ~T |A_A^{(n)}\rangle=\exp{(iQL)}|A_A^{(n)}\rangle, \nonumber \end{align} \tag{ 24 }$$

where in the last step Bloch's theorem (see [14]) was employed, with Q being the modulus of the Bloch wave vector. Here the so-called transfer matrix T is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T =M_{BA}^{-1}M_BM_{AB}^{-1}M_A, \nonumber \end{align} \tag{ 25 }$$

and therefore

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{d1} [T -\exp{(iQL)}I]|A_A^{(n)} \rangle = 0, \nonumber \end{align} \tag{ 26 }$$

where I is the $4\times 4$ unit matrix. The equation above is the secular equation that defines the bulk polariton dispersion relation. From the property that T is an unimodular matrix (det T  =  1) it can be shown [29] that its eigenvalues $\{s_{1},s_{2},s_{3},s_{4}\}$ occur in reciprocal pairs, allowing us to write $s_3 = s_1^{-1}$ and $s_4 = s_2^{-1}$. Therefore, as we conclude below, it is necessary only to find two eigenvalues, which are related to the Bloch wave vectors through [29]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{disp_ge} \exp(iQ_{j}L)=s_{j},\ j=1,\ 2. \nonumber \end{align} \tag{ 27 }$$

By solving equation (26) we find that the eigenvalues satisfy the polynomial equation $s^{4}+\alpha s^{3}+\beta s^{2}+\gamma s+1=0$, where the coefficients are related to the elements of the transfer matrix by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha=\gamma=-Tr(T), \nonumber \end{align} \tag{ 28 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\beta=T_{1}{}_{1}T_{22}+T_{33}T_{44}+(T_{11}+T_{22})(T_{33}+T_{44})-T_{24}T_{42} \nonumber \\ &\qquad -T_{34}T_{43}-T_{23}T_{32} -T_{12}T_{21}-T_{13}T_{31}-T_{14}T_{41}. \nonumber \end{align} \tag{ 29 }$$

After some straightforward algebra the above polynomial can be rearranged as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H^{2}+\alpha H+(\beta-2)=0 \nonumber \end{align} \tag{ 30 }$$

with $H_{i} = s_{i}+s_{i}^{-1}$. This roots of this equation, together equation (24), gives us the dispersion relations for the bulk bands in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{QL12} \cos(Q_{1}L)=H_{1}/2 \quad {\rm and} \quad \cos(Q_{2}L)=H_{2}/2, \nonumber \end{align} \tag{ 31 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{1,2} = \ \frac{-\alpha \pm \sqrt{\alpha^{2}-4(\beta-2)}}{2}. \nonumber \end{align} \tag{ 32 }$$

In this case we have two independent solutions for the Bloch wave vector Q_i, as illustrated later by numerical examples.

It is worth noting that, as a special case of the above theory, we may use equation (17) to consider the polaritons propagating at a single interface (at z  =  a) between two different semi-infinite media (A and B). The general result reduces to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\begin{array}{@{}c@{}} 0 \nonumber \\ B_{xA}^{(1)} \nonumber \\ 0 \nonumber \\ B_{yA}^{(1)} \end{array}\right) = M_{A}^{-1}M_{AB} \left(\begin{array}{@{}c@{}} A_{xB}^{(1)} \nonumber \\ 0 \nonumber \\ A_{yB}^{(1)} \nonumber \\ 0 \end{array}\right), \nonumber \end{align} \tag{ 33 }$$

which leads to a $2\times 2$ matrix secular equation for the polariton modes. After some algebra, we can recover the result given in [19] for the single graphene sheet at the interface between two media, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\frac{\epsilon_{A}}{q_A}+\frac{\epsilon_{B}}{q_{B}}+\frac{i\sigma_{xx}}{\epsilon_A\epsilon_0\omega}) (q_A+q_B-\frac{i\sigma_{yy}\omega}{\epsilon_0 c^2})= \frac{\sigma_{yx}\sigma_{xy}}{(\epsilon_0c)^2} \, , \nonumber \end{align} \tag{ 34 }$$

where we define $q_A =-ik_{zA}$ and $q_B =-ik_{zB}$.

For a study of the surface modes we consider a semi-infinite GPC with the geometry as in figure 1, where the region z  <  0 is taken to be a dielectric medium C. In practice, a surface mode in a photonic crystal is excited by an incident electromagnetic wave either in p-polarization (a TM mode) or s-polarization (a TE mode), e.g. as in an ATR experiment. Here we will consider both the cases of p-polarization and s-polarization in the external medium C.

3.1. Case of p-polarization in the external medium

We consider now a p-polarized wave incident on the PC sample from the external medium C. The electric and magnetic fields are given by $\vec {E_C}=(E_{xC},0,E_{zC})$ and $\vec {H_C}=(0,H_{yC},0)$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{xC} = A_{xC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}, \label{Exc} \nonumber \end{align} \tag{ 35 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{zC} = -\frac{k_x}{k_{zC}}A_{xC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}, \label{Ezc} \nonumber \end{align} \tag{ 36 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{yC} =-\frac{\epsilon_0\omega\epsilon_C}{ik_{zC}} A_{xC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}. \label{Hyc} \nonumber \end{align} \tag{ 37 }$$

Here, $\newcommand{\e}{{\rm e}} k_{zC}=(k_x^2-\epsilon_C\omega^2/c^2){}^{1/2}$ and we assume $\newcommand{\e}{{\rm e}} k_x^ 2>\epsilon_C\omega^2 /c^ 2$ so that k_zC takes only real and positive values. To obtain equations (36) and (37), we used $\newcommand{\e}{{\rm e}} \vec {H}=(\epsilon_0 c^ 2/i\omega) \vec {\nabla} \times\vec {E}$ together with $\vec {\nabla} \cdot \vec {D}=0$. By applying the previous electromagnetic boundary conditions as in equations (10)–(12) to the interface at z  =  0, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{xC}^{(0)}=A_{xA}^{(1)}+B_{xA}^{(1)}, \label{AxBx} \nonumber \end{align} \tag{ 38 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0=A_{yA}^{(1)}+B_{yA}^{(1)}, \label{AyBy} \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\epsilon_0\omega\epsilon_A}{k_{zA}}(A_{xA}^{(1)}-B_{xA}^{(1)})+ \frac{\epsilon_0\omega\epsilon_C}{ik_{zC}}A_{xC}^{(0)}= -\sigma_{xx}(A_{xA}^{(1)}+B_{xA}^{(1)})-\sigma_{xy}(A_{yA}^{(1)}+B_{yA}^{(1)}), \ \ \ \ \ \ \label{r1} \nonumber \end{align} \tag{ 40 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\frac{k_{zA}\epsilon_0c^2}{\omega}(A_{yA}^{(1)}-B_{yA}^{(1)})= \sigma_{yx}(A_{xA}^{(1)}+B_{xA}^{(1)})+\sigma_{yy}(A_{yA}^{(1)}+B_{yA}^{(1)}). \label{r2} \nonumber \end{align} \tag{ 41 }$$

From these coupled equation we can deduce a relationship between the $A_{xA}^{(1)}$ and $B_{xA}^{(1)}$ amplitudes. The result is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{xA}^{(1)}=\lambda B_{xA}^{(1)}, \label{r7} \nonumber \end{align} \tag{ 42 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda= \Big(1-\frac{k_{zA}\epsilon_C}{ik_{zC}\epsilon_A}-\frac{\sigma_{xx}k_{zA}}{\epsilon_0\omega\epsilon_A}\Big) / \Big(1+\frac{k_{zA}\epsilon_C}{ik_{zC}\epsilon_A}+\frac{\sigma_{xx}k_{zA}}{\epsilon_0\omega\epsilon_A}\Big). \label{r8} \nonumber \end{align} \tag{ 43 }$$

Then from equations (40)–(42) we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle B_{yA}^{(1)}=\Delta B_{xA}^{(1)} \quad {\rm and}\quad A_{yA}^{(1)}=-\Delta B_{xA}^{(1)}, \label{r6a} \nonumber \end{align} \tag{ 44 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta=\frac{\sigma_{yx}\omega}{2k_{zA}\epsilon_0c^2}(\lambda+1).\label{r8b} \nonumber \end{align} \tag{ 45 }$$

For the semi-infinite GPC in figure 1, instead of using Bloch's theorem as in the bulk-mode case, we now seek solutions for localized SPP that have decaying amplitudes for the electromagnetic fields with respect to distance from one unit cell to the next interior cell, starting from the surface plane z  =  0. Therefore, taking equation (24) with $Q=i\beta'$, we must have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T |A_A^{(n)}\rangle=\exp{(-\beta'L)}|A_A^{(n)}\rangle, \qquad (z>0) \nonumber \end{align} \tag{ 46 }$$

with Re$(\beta')>0$ (or equivalently $\newcommand{\e}{{\rm e}} |\exp(-\beta'L)| < 1$) as a requirement for localization. On applying this equation to the first layer, we have (following [13])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{r9} T \left(\begin{array}{@{}c@{}} A_{xA}^{(1)} \nonumber \\ B_{xA}^{(1)} \nonumber \\ A_{yA}^{(1)} \nonumber \\ B_{yA}^{(1)} \end{array}\right)=\exp{(-\beta'L)} \left(\begin{array}{@{}c@{}} A_{xA}^{(1)} \nonumber \\ B_{xA}^{(1)} \nonumber \\ A_{yA}^{(1)} \nonumber \\ B_{yA}^{(1)} \end{array}\right). \nonumber \end{align} \tag{ 47 }$$

From this equation, we can find two conditions for localized surface modes. One condition comes from writing out the first and second lines of equation (47) in full, and leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{11} +T_{12}/\lambda -T_{13}\Delta/\lambda +T_{14}\Delta /\lambda=T_{21}\lambda +T_{22} -T_{23}\Delta +T_{24}\Delta, \label{surfe1} \nonumber \end{align} \tag{ 48 }$$

together the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |T_{11} +T_{12}/\lambda -T_{13}\Delta/\lambda +T_{14}\Delta/\lambda|<1. \nonumber \end{align} \tag{ 49 }$$

Also, from the third and fourth lines of (47), we obtain a second implicit dispersion relation for surface modes as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{41}\lambda +T_{42} -T_{43}\Delta +T_{44}\Delta=-T_{31}\lambda -T_{32} +T_{33}\Delta-T_{34}\Delta, \label{surfe2} \nonumber \end{align} \tag{ 50 }$$

together with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |T_{31}\lambda/\Delta +T_{32}/\Delta -T_{33}+T_{34}|<1. \label{surfe2b} \nonumber \end{align} \tag{ 51 }$$

Two brief comments on the above results are useful. First, if we put $\sigma_{xy}=\sigma_{yx}=0$, we just recover the polariton results in our previous paper [9] where there was no applied magnetic field. Secondly, our present results for the surface modes simplify in the regime where retardation effects are small (as at larger in-plane wavenumbers). In this case, choosing $k_{zC} \approx k_x$ and $k_{zA}=i k_x$, following [30], the parameters in equations (43) and (45) reduce to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda= (1-\frac{\epsilon_C}{\epsilon_A}-\frac{\sigma_{xx}ik_x}{\epsilon_0\omega\epsilon_A}) / (1+\frac{\epsilon_C}{\epsilon_A}+\frac{\sigma_{xx}ik_{x}}{\epsilon_0\omega\epsilon_A}), \nonumber \end{align} \tag{ 52 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta=\frac{\sigma_{yx}\omega}{2ik_{x}\epsilon_0c^2}(\lambda+1). \nonumber \end{align} \tag{ 53 }$$

3.2. Case of s-polarization in the external medium

To obtain the surface modes in a semi-infinite GPC in the case of an incident wave with s-polarization (TE) we write $\vec {E_C}=(0,E_{yC},0)$ and $\vec {H_C}=(H_{xC},0,H_{yC})$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{yC} = A_{yC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}, \label{Eyc} \nonumber \end{align} \tag{ 54 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{xC} = -\frac{k_{zC}\epsilon_0 c^2}{i\omega}A_{yC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}, \label{Hxc} \nonumber \end{align} \tag{ 55 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{zC} =\frac{k_x\epsilon_0c^2}{\omega} A_{yC}^{(0)}e^{k_{zC}z}e^{i(k_xx-\omega t)}. \label{Hzc} \nonumber \end{align} \tag{ 56 }$$

To obtain equations (55) and (56) we used $\newcommand{\e}{{\rm e}} \vec {H}=$ $\newcommand{\e}{{\rm e}} (\epsilon_0 c^ 2/i\omega) \vec {\nabla} \times\vec {E}$. Now, by applying the same boundary conditions, (10)–(12), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0=A_{xA}^{(1)}+B_{xA}^{(1)}, \label{eq1} \nonumber \end{align} \tag{ 57 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{yC}^{(0)}=A_{yA}^{(1)}+B_{yA}^{(1)}, \label{eq2} \nonumber \end{align} \tag{ 58 }$$

for the electrical fields amplitudes, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\epsilon_0\omega\epsilon_A}{k_{zA}}(A_{xA}^{(1)}-B_{xA}^{(1)})= -\sigma_{xx}(A_{xA}^{(1)}+B_{xA}^{(1)})-\sigma_{xy}(A_{yA}^{(1)}+B_{yA}^{(1)}), \label{eq3} \nonumber \end{align} \tag{ 59 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\frac{k_{zA}\epsilon_0c^2}{\omega}(A_{yA}^{(1)}-B_{yA}^{(1)})+\frac{k_{zC}\epsilon_0c^2}{i\omega}A_{yC}^{(0)} \nonumber \\ =\sigma_{yx}(A_{xA}^{(1)}+B_{xA}^{(1)})+\sigma_{yy}(A_{yA}^{(1)}+B_{yA}^{(1)}), \label{eq4} \nonumber \end{align} \tag{ 60 }$$

for the magnetic fields amplitudes. In similar way as for p-polarization, we have from (57), $B_{xA}^{(1)}=-A_{xA}^{(1)}$. Therefore, applying it in (59) and (60) we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{2\epsilon_0\omega\epsilon_A}{k_{zA}}(A_{xA}^{(1)})= -\sigma_{xy}(A_{yA}^{(1)}+B_{yA}^{(1)}), \label{eq5} \nonumber \end{align} \tag{ 61 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\frac{k_{zA}\epsilon_0c^2}{\omega}(A_{yA}^{(1)}-B_{yA}^{(1)})+\frac{k_{zC}\epsilon_0c^2}{i\omega}A_{yC}^{(0)}= \sigma_{yy}(A_{yA}^{(1)}+B_{yA}^{(1)}). \nonumber \end{align} \tag{ 62 }$$

Now, using (58), we can rewrite (3.2) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\frac{k_{zA}\epsilon_0c^2}{\omega}(A_{yA}^{(1)}-B_{yA}^{(1)})+\frac{k_{zC}\epsilon_0c^2}{i\omega}(A_{yA}^{(1)}+B_{yA}^{(1)})= \sigma_{yy}(A_{yA}^{(1)}+B_{yA}^{(1)}). \nonumber \end{align} \tag{ 63 }$$

After some algebra we can find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{yA}^{(1)}(1-\frac{k_{zC}}{ik_{zA}}+\frac{\sigma_{yy}\omega}{k_{zA}\epsilon_0 c^2})=B_{yA}^{(1)}(1+\frac{k_{zC}}{ik_{zA}}-\frac{\sigma_{yy}\omega}{k_{zA}\epsilon_0 c^2}), \nonumber \end{align} \tag{ 64 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{yA}^{(1)}=\lambda' B_{yA}^{(1)}, \nonumber \end{align} \tag{ 65 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Lambl} \lambda'=\frac{(1+\frac{k_{zC}}{ik_{zA}}-\frac{\sigma_{yy}\omega}{k_{zA}\epsilon_0 c^2})}{(1-\frac{k_{zC}}{ik_{zA}}+\frac{\sigma_{yy}\omega}{k_{zA}\epsilon_0 c^2})}. \nonumber \end{align} \tag{ 66 }$$

By using (57) and (58) we can find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{xA}^{(1)}=-\Delta' B_{yA}^{(1)}, \nonumber \end{align} \tag{ 67 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle B_{xA}^{(1)}=\Delta' B_{yA}^{(1)}, \nonumber \end{align} \tag{ 68 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Deltal} \Delta'=\frac{\sigma_{xy}k_{zA}}{2\epsilon_0\omega\epsilon_A}(\lambda'+1). \nonumber \end{align} \tag{ 69 }$$

The conditions for having localized s-polarized surfaces modes can be obtained using the transfer matrix in a similar fashion to the analysis in the previous subsection for the case of p-polarization. The analogous equations to (48)–(51) are found to be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (T_{11} -T_{12})\Delta' -T_{13}\lambda' -T_{14}= (T_{22}-T_{21})\Delta' +T_{23}\lambda' +T_{24}, \nonumber \end{align} \tag{ 70 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (T_{32} -T_{31})\Delta'/\lambda' +T_{33}+T_{34}/\lambda'=(T_{42} -T_{41})\Delta' +T_{43}\lambda' +T_{44}, \nonumber \end{align} \tag{ 71 }$$

with the localization conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |(T_{11} -T_{12}) -T_{13}\lambda'/\Delta' -T_{14}/\Delta'|>1, \nonumber \end{align} \tag{ 72 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |(T_{32} -T_{31})\Delta'/\lambda' +T_{33}+T_{34}/\lambda'|>1. \nonumber \end{align} \tag{ 73 }$$

Again there are simplifications in the case of weak retardation. Specifically, equations (66) and (69) reduce to give

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda'= \frac{2ik_{x}\epsilon_0 c^2}{\sigma_{yy}\omega} -1 \nonumber \end{align} \tag{ 74 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta'=-\frac{\sigma_{xy}}{\sigma_{yy}}\left(\frac{k_x^2}{\epsilon_A(\omega/c)^2} \right). \nonumber \end{align} \tag{ 75 }$$

In this section we apply the previous theory to a semi-infinite GPC to illustrate the results for the bulk and surface polaritons. We take both of the constituent materials, A and B, to be silica (SiO₂), with a dielectric constant $\newcommand{\e}{{\rm e}} \epsilon_A=\epsilon_B \equiv \epsilon =2.2$, as appropriate for THz frequencies [31]. We mention here, as a first approximation, our results will not include the effects of absorption near 4 THz and 8 THz due to SiO₂, since our goal here is to simulate a lossless material and focus on the role of the graphene layers in a magnetic field. Indeed, absorption processes will modify the properties of polaritons most sensitively when the frequencies lie close to a resonance. Generally, absorption leads to attenuation of the polaritons (surface and bulk modes) as they propagate through the crystal (see, e.g. [32]). The absorption effects could be studied in future work.

We choose vacuum ($\newcommand{\e}{{\rm e}} \epsilon_ {C}= 1$) for the external medium. We ignore the role of absorption on the silica since it is typically small for other distinct crystalline forms and composition [33] in the temperature and frequencies of interest in this work. Also we choose the temperature T  =  300 K and a chemical potential $\mu_c= 0.2$ eV for the graphene sheets [26–28]. The other physical parameters used here are a  =  b  =  10 $\mu$m and applied magnetic field B₀  =  1.5 T. We have put our expressions for the optical conductivity in terms of the fine structure constant $\newcommand{\e}{{\rm e}} e^2/4\pi \epsilon_0 \hbar c$ ($\approx$1/137), where c is the light speed in a vacuum. This constant arises naturally from our equations. We consider the electron scattering rate (the inverse of the intrinsic relaxation time) to be $\Gamma =0.11$ meV [21].

In the following figures, we present numerical dispersion relations for plasmon polaritons in the 1D GPC, focusing on the effects of the applied magnetic field. The bulk modes occur here as bands separated by gaps and they are shown as the shaded areas (brown areas in the online version). These bulk bands are bounded by the curves for Q_jL  =  0 and $Q_jL = \pi$, where Q_jL is given in equation (31) with L  =  a  +  b. The surface modes, characterized by black lines, lie between and above the bulk bands, subject to the constraint $\newcommand{\e}{{\rm e}} k_x^ 2>\epsilon_C\omega^2 /c^ 2$. Consequently any surface mode in these plots lies to the right of the vacuum light line (VLL), shown as a dashed straight line, given by $\newcommand{\e}{{\rm e}} \omega=ck_x/\sqrt{\epsilon_ {0}}$. For reference, we also show the dielectric light line (DLL) for the PC given by $\newcommand{\e}{{\rm e}} \omega=ck_x/\sqrt{\epsilon}$. We have chosen to plot the frequency in THz units and k_xa is the scaled 2D in-plane wave vector.

The bulk and surface plasmon-polariton dispersion relations for the frequency (f ) region up to 10 THz are presented in figure 3 for the in-plane wave vector ranging from k_xa  =  −5.0 to k_xa  =  5.0. A reciprocal propagation effect is found for the bulk bands, as expected, with respect to the sign of the wave-vector component, but this is not the case for the surface modes where we have a non-reciprocity effect. In describing the behaviour of these modes we conveniently separate this spectrum into three distinct regions: (a) to the right (or left) side of the DLL for k_x  >  0 (or k_x  <  0) and frequencies greater than 1.0 THz; (b) inside the region bounded by the the DLLs for k_x  >  0 and k_x  <  0 with f   >  1.0 THz; and (c) the region where the frequencies are less than 1 THz.

Zoom In Zoom Out Reset image size

In region (a) in figure 3 we have two bulk bands that start at the DLL and merge into a narrow band approximating a straight line with a constant group velocity for high values of k_xa. Those bulk bands show a behaviour that is reciprocal with respect to k_x. However, we also have five surface modes, labeled by S1, S2, S3, S4 and S5, that are non-reciprocal with respect to the substitution $k_x \rightarrow -k_x$. The explanation for this behaviour is related to the symmetry properties of the optical conductivity tensor. The origin is similar to the well-known non-reciprocity in surface modes [34] and for waves in electromagnetic systems in general [35]. The modes S3, S4 and S5 are p-polarized modes in the external medium (vacuum), while S2 and S3 are s-polarized modes in this external medium. Note that these latter two modes have a particular limit when k_x (or  −k_x) goes to infinity (minus infinity): specifically they tend to f ₃  =  8.32 THz (long-short-dashed line in figure 3) draw attention here to the fact that we have s-polarized (TE) plane surface polariton modes, for which there is no counterpart in non-graphene PCs (or dielectric superlattices) [13, 14]. We also point out the existence of two other very interesting surface modes, labeled by S6 and S7, which can more easily be studied by expanding the frequency region from zero to 1 THz as done in figure 4.

Zoom In Zoom Out Reset image size

In region (b) we have a bulk region terminating at the DLL, for k_x  >  0 or k_x  <  0. In figures 3 and 4 we can also see the so-called graphene induced gap, by analogy with results in [8, 9], but in the present case the gap starts at about 0.32 THz and finishes at about 2.48 THz for k_xa  =  0. Naturally, this gap depends on the GPC parameters (such as layer thicknesses, dielectric constants, etc), but here we have the added feature that it depends on the applied magnetic field as well. For the case where there is zero applied magnetic field, this gap starts at the origin [9], i.e. where k_xa  =  0 and f   =  0.

The low-frequency region (c) is shown in figure 4. Here we can see the specific behaviour for the bulk modes close to the frequency values f ₁  =  0.39 THz and f ₂  =  0.55 THz. The bulk bands, labeled by B1 and B2 are reciprocal and tend to f ₂ when k_x goes to infinity. This value represents a resonance frequency corresponding to the maximum in Im(${\sigma_{xx}}$), or the maximum in Re(${\sigma_{xy}}$). The other bulk bands, labeled by B3 and B4 are also reciprocal and they tend to f ₁ at larger magnitude of k_xa. Here f ₁ is another resonance frequency that corresponds to the minimum in Im${\sigma_{xx}}$ (or the minimum in Re${\sigma_{xy}}$). Clearly, these bulk bands depend on the resonant frequencies f ₁ and f ₂, and consequently on the intensity of the applied magnetic field. These bulk modes are very different from the bulk modes in dielectric (semiconductor) PCs, as well as those reported in previous works [13, 14]. We infer that these bulk modes have this different behaviour due to the absorption, as characterized by the scattering rate $\Gamma$ in equation (2), in the graphene sheets. Also, we see that the applied magnetic field is a very important parameter for tuning the bulk (and surface) polariton modes in the GPC. We note the small gap region that starts at 0.32 THz and finishes at 0.51 THz, for k_xa  =  0, and it is limited by the lowest bulk band and the DLLs. It depends sensitively on the applied magnetic field (see figure 5). On the other hand, for low frequencies, below to 0.05 THz (or 50 GHz), we have more two bulk bands, labeled by B5 and B6 (which are reciprocal), and tend to 0.045 THz when |k_xa| becomes large. The surface modes in the region (c) are also quite different. We have three surface modes, labeled by S8, S7, and S6. The mode S8 starts at f   =  0.38 THz for k_xa  =  0.11 and finishes at f   =  0.49 THz for k_xa  =  0.133. The surface mode S7 starts at f   =  0.51 THz for k_xa  =  0.139 and finishes at the VLL. These modes, S8 and S7, have a behaviour that is typical for polariton surface modes when we consider the damping (or absorption) in the dielectric function [36]. Indeed, here we have taken account of an equivalent damping effect through the scattering rate $\Gamma$ in the optical conductivity tensor (see equations (2) and (3)). We should mention here that these surface modes are very sensitive to this parameter. For example, we will not have surface modes if $\Gamma =1$ meV.

Zoom In Zoom Out Reset image size

The non-reciprocity aspects of the surface modes, as described above, do not change with the increasing of the magnetic field. However, these modes will be shifted (except for the surface modes S7 and S8 that will disappear). We can interpret this using figure 2, because, essentially, we have a shift in the optical conductivity (including the peaks and dips in the imaginary and real parts of $\sigma$). Therefore, we present now a study to find the behaviour of these modes for bulk bands with the applied magnetic field. This behaviour will be similar for all the spectra.

Finally, in figure 5 we show the bulk polariton modes (and the band gap) for the region below 1 THz, for k_xa  =  0 plotted as a function of the applied magnetic field. Here, we can see that it is possible to tune the bulk bands (and gap) by changing the applied magnetic field intensity. Note that the frequencies of these bulk bands will increase approximately linearly as a function of the applied magnetic field when k_xa  =  0. In general terms, we comment that the absence of an applied magnetic field in this system would make this spectrum more simple, mainly in the region around f   <  0.6 THz, in figure 4. Further, we would not have surface modes for s-polarization [9] in that case, as discussed also by Albuquerque and Cottam [13, 14]. On the other hand, the bulk spectrum with applied magnetic field is rather similar to the case in the absence of an applied magnetic field for p-polarization (see [9] for more detail).

In summary, we have presented a general theory for the plasmon-polaritons in graphene PCs subjected to an external applied magnetic field. Specifically we have considered a periodic photonic crystal consisting of alternating unit cells arranged in a periodic array with embedded graphene sheets at the interfaces and a static applied magnetic field in the direction perpendicular to the layers. On comparing our results with previous work we find that the inclusion of an external magnetic field in the GPC will affect strongly the dispersion relations of the polaritons, mainly at low frequencies (below about 1 THz, for B  =  1.5 T). Also, we find a novel nonreciprocity effect for the propagation of the surface modes. Further, there is the possibility of having s-polarized (TE) surface polariton modes that no have counterpart in the dispersion relations in semiconductor (or dielectric) PCs [13].

In the context of applications, we mention that surface magnetoplasmons (or magnetoplasmon polaritons) in multilayer structures are typically difficult to realize because they require very high magnetic fields. The explanation is that in order to observe the effects due to an external magnetic field for metals, for example, the frequencies would need to be of the order of 10¹⁶ Hz, which is the order of the characteristic cyclotron frequency in metals [37]. The characteristic cyclotron frequency is a function of the applied magnetic field, and, therefore, in the case of metals, it needs a magnetic field of order 10³ T, which is exceptionally difficult to realize in laboratories. One solution to avoid this generally unreachable magnetic field strength is to use doped semiconductors (instead of metals) where the characteristic cyclotron frequency is in the THz regime (for a review see [38]). In this case, a more modest external magnetic field of less than 2 T will suffice. From our results, by using GPCs instead semiconductor PCs, we can see that there is the possibility for the excitation of surface and bulk magnetoplasmon polaritons in this same frequency range and the same order for the specified external magnetic field.

M S Vasconcelos thanks the Department of Physics and Astronomy at the University of Western Ontario for hospitality during his sabbatical as visiting professor. This study was financed in part by the Coordenação de Aperfeioamento de Pessoal de Nível Superior (CAPES) of Brasil (Finance Code 88881.172293/2018-01) and the Natural Sciences and Engineering Research Council (NSERC) of Canada (Grant RGPIN-2017-04429).

Here, we provide the definitions for the matrix blocks constituting the interface matrix in equation (19). First, $M^{11}_{AB}$ and $M^{22}_{AB}$, which involve the diagonal elements of the graphene conductivity tensor, are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{11}_{AB} = \frac{1}{2}\left(\begin{array}{@{}cc@{}} 1+\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}}+ \frac{\sigma_{xx}k_{zA}}{\epsilon_A\epsilon_0\omega} & 1-\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}}+\frac{\sigma_{xx}k_{zA}}{\epsilon_A\epsilon_0\omega} \nonumber \\ 1-\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}}-\frac{\sigma_{xx}k_{zA}}{\epsilon_A\epsilon_0\omega} & 1+\frac{\epsilon_{B}k_{zA}}{\epsilon_{A}k_{zB}}-\frac{\sigma_{xx}k_{zA}}{\epsilon_A\epsilon_0\omega} \end{array}\right), \nonumber \end{align} \tag{ A.1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{22}_{AB} = \frac{1}{2} \left(\begin{array}{@{}cc@{}} 1+\frac{k_{zB}}{k_{zA}}+\frac{\sigma_{yy}\omega}{\epsilon_0 c^2 k_{zA}} & 1-\frac{k_{zB}}{k_{zA}}+\frac{\sigma_{yy}\omega}{\epsilon_0c^2k_{zA}} \nonumber \\ 1-\frac{k_{zB}}{k_{zA}}-\frac{\sigma_{yy}\omega}{\epsilon_0 c^2 k_{zA}} & 1+\frac{k_{zB}}{k_{zA}}-\frac{\sigma_{yy}\omega}{\epsilon_0c^2k_{zA}} \end{array}\right). \nonumber \end{align} \tag{ A.2 }$$

We note that these matrices are similar to those in [9] for p-polarization and s-polarization, respectively. The other matrices $M^{12}_{AB}$ and $M^{21}_{AB}$, which involve the off-diagonal elements of the graphene conductivity tensor, are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{12}_{AB} = \frac{1}{2} \left(\begin{array}{@{}cc@{}} \frac{\sigma_{xy}k_{zA}}{\epsilon_A\epsilon_0\omega} & \frac{\sigma_{xy}k_{zA}}{\epsilon_A\epsilon_0\omega} \nonumber \\ -\frac{\sigma_{xy}k_{zA}}{\epsilon_A\epsilon_0\omega} & -\frac{\sigma_{xy}k_{zA}}{\epsilon_A\epsilon_0\omega} \nonumber \\ \end{array}\right), \nonumber \end{align} \tag{ A.3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{21}_{AB} = \frac{1}{2} \left(\begin{array}{@{}cc@{}} \frac{\sigma_{yx} \omega}{\epsilon_0 c^2 k_{zA}} & \frac{\sigma_{yx} \omega}{\epsilon_0 c^2 k_{zA}} \nonumber \\ -\frac{\sigma_{yx} \omega}{\epsilon_0 c^2 k_{zA}} & -\frac{\sigma_{yx} \omega}{\epsilon_0 c^2 k_{zA}} \end{array}\right). \nonumber \end{align} \tag{ A.4 }$$