The feasibility of composite right/left-handed (CRLH) metamaterial waveguides based upon graphene plasmons is demonstrated via numerical simulation. Designs are presented that operate in the terahertz frequency range along with their various dimensions. Dispersion relations, radiative and free-carrier losses, and free-carrier based tunability are characterized. Finally, the radiative characteristics are evaluated, along with its feasibility for use as a leaky-wave antenna. While CRLH waveguides are feasible in the terahertz range, their ultimate utility will require precise nanofabrication, and excellent quality graphene to mitigate free-carrier losses.

Graphene shows considerable promise as an optical material in the terahertz and mid-infrared frequency ranges as a result of the strong optical response associated with intraband conduction of Dirac fermions. Two-dimensional plasmons in single and multi-layer graphene structures have been extensively described theoretically and via numerical simulations.1–5 Furthermore, direct evidence of their existence has been shown in the THz and mid-IR using near-field probe techniques,6–9 as well as using infrared absorption spectroscopy in patterned structures (micro- and nano-ribbons, discs, dipoles, etc.) to excite dipolar modes associated with standing-wave plasmons.10–12 For sheet carrier concentrations of 10121013 cm−2, THz resonances can be achieved for structures with widths of several microns;10,11 if the sizes are reduced to tens of nanometers, resonances in the 3–12 μm range can be obtained.12–14 By leveraging these effects, graphene has been demonstrated as a THz modulator15–17 and has been proposed for strong light-matter coupling3,12 as well as for resonant and leaky-wave antennas (LWAs).18–21 

In this paper, we explore the feasibility of graphene for composite right/left-handed (CRLH) metamaterial waveguides in the terahertz and far-IR frequency range (1–15 THz). The CRLH metamaterial concept is most easily understood using transmission-line (TL) theory. A conventional (right-handed) transmission-line has a series inductance (LR) and a shunt capacitance (CR); by loading the line with distributed series capacitance (CL) and shunt inductance (LL), the line becomes highly dispersive and exhibits backward wave (left-handed) propagation.22–24 The dispersion relation of the angular frequency ω vs. the propagation constant β is shown for a typical CRLH line in Fig. 1. A stopband is evident at β = 0, where the band edge frequencies are given by the shunt resonance ωsh=(LLCR)12 and series resonance ωse=(LRCL)12. If the line is properly designed, these resonance frequencies can be made equal so that the stopband closes; this is known as the balanced condition. A balanced CRLH line is particularly useful for traveling wave devices since all modes exhibit non-zero group velocity, and the characteristic impedance of the line is independent of frequency. One- and two-dimensional metamaterials have been widely explored in the microwave frequency range and have been used to demonstrate a variety of guided-wave devices (e.g., multi-band and enhanced bandwidth components, power combiners/splitters, compact resonators, phase shifters, and phased array feed lines), as well as radiating devices (e.g., 1D and 2D resonant and leaky-wave antennas).23,25 More recently, CRLH metamaterial waveguides and metasurfaces have been demonstrated in the THz frequency range, both for passive waveguides and in active devices integrated with THz quantum-cascade laser gain material.26,27

FIG. 1.

CRLH dispersion relationship which exhibits both right- and left-handed characteristics. A stopband exists between ωsh and ωse—the shunt and series resonance frequencies, respectively. The inset shows the characteristic CRLH transmission line unit cell, which contains series inductance LR, shunt capacitance CR, series capacitance CL, and shunt inductance LL.

FIG. 1.

CRLH dispersion relationship which exhibits both right- and left-handed characteristics. A stopband exists between ωsh and ωse—the shunt and series resonance frequencies, respectively. The inset shows the characteristic CRLH transmission line unit cell, which contains series inductance LR, shunt capacitance CR, series capacitance CL, and shunt inductance LL.

Close modal

We begin by considering the propagation of plasmons by infinite, unpatterned graphene sheets on a dielectric half-space. If coupling with free-space propagating electromagnetic waves is neglected, graphene plasmons are described by a dispersion relation for the in-plane wavenumber q

q=i(ε1+ε2)ωσ(ω),
(1)

where ε1 and ε2 are the permittivities of the two dielectric regions, and σ(ω) is the ac sheet conductivity of the graphene. Throughout this paper, we model the graphene conductivity using the Drude model as given by28 

σ2D(ω)=e2EFπ2τ1iωτ,
(2)

where EF is the Fermi level, vF=106 m/s is the Fermi velocity, e is the fundamental charge, is the reduced Planck's constant, and τ is the Drude momentum relaxation time. The Fermi level is related to the sheet carrier density n by

EF=vFπ|n|.
(3)

Inserting the Drude conductivity into (1) gives the dispersion relation

q(ω)=π2(ϵ1+ϵ2)e2EF(1+iτω)ω2,
(4)

which exhibits the characteristic ωq1/2 behavior associated with 2D plasmons and the ωn1/4 behavior associated with a Dirac bandstructure.2 

It has been shown that the Drude model for conductivity does an excellent job of capturing the essential behavior of graphene plasmons, particularly when the wavenumber q is modest compared to the Fermi wavenumber (qEF/vF),29 and at sufficiently low temperatures such that kBTEF. Both conditions are well satisfied for the doping levels of graphene considered here where n1012 cm−2 (EF0.13 eV). Since we will consider plasmons with frequencies of at most 15 THz (ω62 meV), interband transitions will be forbidden since ω2EF. More accurate conductivity expressions that include interband transitions and comparisons with random-phase-approximation (RPA) based models are given in Ref. 29; treatments of nonlocal conductivity (i.e., spatial dispersion) are given in Refs. 30 and 31. It has been shown that classical local models for the conductivity are generally accurate for graphene nanostructures larger than ∼20 nm.32 Our confidence in using the local Drude conductivity is further bolstered by reports which show that even for 60 nm wide graphene nanoribbons, use of a nonlocal conductivity model (where conductivity σ(ω,q) is a function of both frequency and wavenumber) only shifts resonances by 5% or less for our considered carrier concentrations.33 For larger nanostructures at lower frequencies, the shifts are smaller, since the in-plane wavenumbers q involved are smaller. As a final check, we have performed a post-computational analysis on our simulation results by using Fourier analysis to extract the spectral power of the current density distribution. For the frequencies and dimensions considered here, we confirmed that 95% or more of the spectral power is at wavenumbers sufficiently small that corrections due to spatial dispersion are a few percent or less.30,31

Since CRLH metamaterials are typically described using transmission-line (TL) models, it is useful to begin our discussion by describing the fundamental lateral plasmon mode (m = 0) on a microribbon of width w using a transmission line (TL) model (e.g., as in Refs. 31 and 34). This mode is schematically illustrated in Fig. 2(a) along with its dispersion relation. In the low-damping limit (ωτ1), the kinetic inductance of the electrons dominates the graphene impedance; this can be represented by an equivalent series inductor LR in the circuit diagram in Fig. 1. The fringing E-field longitudinally couples charges along the graphene ribbon and plays the role of the shunt capacitor CR. Effective values for these circuit elements are given in the  Appendix.

FIG. 2.

Conceptual design progression for CRLH graphene plasmonic waveguide along with their associated circuit models (for one unit cell) and dispersion relations. (a) and (b) The fundamental m = 0 even mode with and without series gaps. (c) and (d) The higher order m = 1 odd mode with and without series gaps. The latter case shown in (d) exhibits CRLH behavior and is the basis for the designs presented here.

FIG. 2.

Conceptual design progression for CRLH graphene plasmonic waveguide along with their associated circuit models (for one unit cell) and dispersion relations. (a) and (b) The fundamental m = 0 even mode with and without series gaps. (c) and (d) The higher order m = 1 odd mode with and without series gaps. The latter case shown in (d) exhibits CRLH behavior and is the basis for the designs presented here.

Close modal

In order to realize a CRLH line, a shunt inductance LL and a series capacitance CL must be added. It is straightforward to create a distributed series capacitance by introducing repeated gaps of size a and period p along the graphene microribbon. This adds a high-pass characteristic to the dispersion relation for the m = 0 mode, with a cutoff frequency equal to the series resonance ωse (see Fig. 2(b)).

In microstrip implementations, the shunt inductance LL is typically achieved using a stub and/or via to the ground plane. At THz frequencies, such techniques are often lossy and difficult to fabricate; furthermore, they are incompatible with the geometry of a single graphene microribbon. However, an effective shunt inductance arises naturally if we consider propagation of the first higher-order lateral mode (labeled m = 1) of the microribbon (also known as an edge mode), which is associated with transverse conduction currents. This can be modeled by two transmission lines in parallel, coupled together by an inductance 2LL (see Fig. 2(c)). The effective shunt resonant frequency ωsh is given by the cutoff frequency for this higher order mode, which occurs when the width w is approximately equal to one half of the plasmon wavelength, i.e., q(ωsh)π/w. This approach to obtaining LL has been used for passive and active THz microstrip-type waveguides.26,27,35

When gaps are introduced, the dispersion of the m = 1 mode takes on a CRLH character, with both right-handed and left-handed propagating characteristics. The double transmission line model is shown in Fig. 2(d). When this line is excited in its odd mode, the central symmetry plane acts as a virtual ground, so that each branch is equivalent to the circuit in Fig. 1. With proper design, one can obtain the condition ωse=ωsh, and we say, the structure is “balanced”; the stopband closes and propagating waves exist with non-zero group velocity even at β = 0.

Numerical simulations were performed using a commercial finite element full-wave electromagnetic solver (Ansys HFSS). The eigenfrequency solver was applied to a single unit cell of length p with periodic Floquet-Bloch boundary conditions in the axial direction for a fixed propagation constant β; this was used to extract the complex eigenfrequency ω=ω+iω. This procedure allowed the retrieval of the dispersion relationships ω(β).

The graphene is modeled as a 2D impedance boundary condition of Zs=σ(ω)1 at the dielectric halfspace interface, where the conductivity σ is taken from Eq. (2).36 The simulation space consists of a cylindrical “airbox” of 15 μm radius surrounding the graphene microribbon. The upper semi-cylinder has permittivity ε1 (assumed to be vacuum throughout), and the lower semi-cylinder has permittivity ε2 which represents a semi-infinite substrate. For radiating modes, it is important to suppress reflections at the transverse airbox boundaries; this is accomplished by using a layered impedance boundary condition. Various airbox sizes were considered to verify that their dimensions did not affect the complex eigenfrequency solution.

Since the conductivity of graphene is dependent on frequency, an iterative approach was used to achieve a self-consistent solution. The algorithm is as follows: first, for a given value of β, an initial radian frequency ω0 is guessed using the lumped element circuit model; this value is used to calculate σ(ω0) for graphene which is input into the simulation, and the complex eigenfrequency ωi is solved numerically for each i-th iteration. Second, the conductivity σ(ωi) is used to re-solve for the eigenfrequency ωi+1. The process is iterated until the real part of the complex eigenfrequency converges within 20 GHz. As verification, this method is seen to accurately reproduce the analytic dispersion relation for an unpatterned graphene microribbon in its fundamental m = 0 mode in Fig. 3.

FIG. 3.

Simulated dispersion relationship for Design A CRLH structure (m = 1 with series gaps) with a center frequency of 2.8 THz (red) and dimensions shown. Also shown is the dispersion for the m = 0 plasmonic mode on a graphene microribbon of identical width (green). Shown in the insets are plots of the graphene charge density ρ at the shunt and series resonance. Light lines (ω=cβ/n) for vacuum (n = 1) and silicon (n = 3.41) are also plotted in black and gray, respectively.

FIG. 3.

Simulated dispersion relationship for Design A CRLH structure (m = 1 with series gaps) with a center frequency of 2.8 THz (red) and dimensions shown. Also shown is the dispersion for the m = 0 plasmonic mode on a graphene microribbon of identical width (green). Shown in the insets are plots of the graphene charge density ρ at the shunt and series resonance. Light lines (ω=cβ/n) for vacuum (n = 1) and silicon (n = 3.41) are also plotted in black and gray, respectively.

Close modal

The imaginary part of the eigenfrequency gives information both about Drude loss and radiative loss (for modes within the light line (β<nω/c)). The modal quality factor Q is obtained from the complex eigenfrequency according to Q=ω/2ω, and with knowledge of the group velocity vg=dω/dβ of the mode, the power loss coefficient per unit length is given by α=ω/Qvg.

However, a complication arises if the finite element eigenfrequency solver is used to directly extract the Drude losses because the graphene impedance is treated as frequency independent during any given iteration. Essentially, the imaginary part of the impedance is not allowed to vary with frequency during the eigenfrequency solution process (equivalent to finding a root in complex frequency space). This results in the eigenfrequency solver underestimating the Drude contribution to the losses by as much as a factor of two (depending upon the value of β). For example, for the m = 0 guided mode, the HFSS solver calculates Q=ωτ/2, instead of the expected Q=ωτ that can be calculated analytically from Eq. (4). This issue is addressed in detail in Sec. IV A via a hybrid method, where only the radiative losses are extracted numerically, and an analytic expression is used for the Drude losses based upon a transmission line model.

An exemplar CRLH waveguide was designed to operate with a center frequency of 2.8 THz, with dimensions as follows: unit cell length p = 2.5 μm, capacitive gap size a = 100 nm, and microribbon width w = 4.55 μm. Design proceeds as follows. For a given carrier density, the width w of the microribbon is chosen first to set ωsh to the frequency of choice. The circuit model is used to give a coarse estimate (see the  Appendix); the value is then refined using a 2D full-wave simulation for a finite width ribbon. Then, gaps are added to the design, with an initial guess for the period p and gap size a again provided by the circuit model. Then, the dimensions are fine tuned using full-wave simulations so that the CRLH structure is nearly balanced (i.e., ωsh=ωse). For this baseline design, which we will refer to as “Design A,” the carrier density was taken as n=3×1013 cm−2 (EF=0.64 eV), a value which is readily experimentally achievable. We also assumed that the graphene structure was at the interface of a dielectric half-space between vacuum and a lossless silicon substrate (εr=11.6). This was chosen somewhat arbitrarily, and our results should be understood to be applicable to a variety of dielectric substrates (i.e., SiO2, boron nitride (BN), etc.), provided the appropriate modifications are made to the dimensions, plasmon damping rates, and substrate dielectric and free-carrier losses. Similarly, we have neglected any coupling with optical phonon Reststrahlen bands in the substrates;12 such coupling is negligible at THz frequencies but is important in the mid-IR (i.e., λ10μm for SiO2 and BN). Unless otherwise specified, momentum relaxation time was set to be τ = 1 ps. However, simulations showed that provided ωτ>1, the dispersion relation is not sensitive to the particular value of τ. The simulated dispersion relation is presented in Fig. 3 and exhibits the characteristics shown in Fig. 1. Left-handed propagation is exhibited between 1.8–2.8 THz. The structure is nearly balanced; however, if closely examined, there is a slight ∼20 GHz gap that remains at β = 0. At this point, distinct series and shunt resonance modes can be observed, and their charge densities are shown in the inset of Fig. 3. However, when one considers non-zero values of β, the propagating modes are hybrids of these two basis states. For comparison, also shown in this figure is the dispersion of the fundamental m = 0 mode for a ribbon without gaps, which exhibits the characteristic ωβ1/2 behavior. As expected, the shunt resonance frequency occurs close to the half-wavelength condition (i.e., β(ωsh)w0.9π for the m = 0 mode). The residual 20 GHz gap is likely due to the imbalances in radiative losses associated with the shunt and series branches of the transmission line.37,38

Several factors drive the choice of p. While ideally p should be as small as possible to obey the metamaterial homogenization condition, in practice as long as the period satisfies, the condition p<(2π/q)/2 no higher order Bragg scattering will contribute to propagation. The finite-sized unit cell also gives rise to a high-pass filtering effect, where the series capacitance and shunt inductance behave as a high pass filter with a cutoff of ωcL=1/2LLCL, which is seen to be approximately 1.8 THz in this design. Reducing the unit cell size p results in a larger value of LL and CL (required to keep ωsh and ωse constant) which reduces the cutoff frequency. This can be used to engineer the group velocity and leaky-wave bandwidth of the CRLH waveguide. From a practical standpoint, reducing the unit cell is challenging, since the series capacitance has only a logarithmic dependence on the gap size a. For example, our simulations showed that modifying Design A by reducing the period from 2.5 μm to 1 μm requires reducing the gap size from 100 nm to 1.3 nm to maintain the balanced condition.

One of the attractive features of graphene is the ability to tune the plasmon dispersion via electrostatic gating, chemical doping, or optical pumping. This property can be exploited to tune CRLH waveguides while approximately maintaining their balanced character. Simulated dispersion relations for various concentrations for our Design A structure are shown in Fig. 4. Based on the graphene dispersion relationship presented in (4), we expect a fourth root scaling of frequency with carrier density (i.e., ωn1/4). Indeed, this is observed where two orders of magnitude change in carrier density shift the center frequency by approximately a factor of three. This will prevent tuning of a single structure into different frequency ranges—say from the THz to the mid-IR. However, this tuning is more than sufficient to tune a graphene structure across the leaky-wave bandwidth (i.e., hundreds of GHz) for a beam-scanning or phase shifter application. For example, a phase shift of π radians can be obtained across a single unit cell by changing the carrier density by a factor of ten. Fig. 4 shows that the balanced condition is approximately preserved even as the carrier density is changed. Indeed, this is predicted by the analytic transmission-line treatment given in the  Appendix for which only geometric parameters determine the balanced condition. However, it is seen that slight unbalancing occurs as the carrier density is tuned, whether this is within acceptable limits will depend upon the application.

FIG. 4.

Dispersion relationships for Design A for various graphene carrier densities labeled in units cm−2.

FIG. 4.

Dispersion relationships for Design A for various graphene carrier densities labeled in units cm−2.

Close modal

Table I presents the dimensions of several CRLH waveguides designed for center frequencies of 2.8 THz, 5 THz, and 12 THz—labeled designs A, B, and C, respectively. Dispersion relations for Design A, B, and C are plotted in Fig. 5, along with transmission line model fits obtained using the model and parameters given in the  Appendix. It is seen that excellent agreement is obtained although the use of several adjustable parameters is required to properly fit the numerical results. Unfortunately, the q1/2 nature of graphene dispersion means that the width w must be reduced quadratically to achieve linear increases in frequency. This is seen in our designs, where frequency is scaled up by a factor of 4.2 from 2.8 THz (Design A) to 12 THz (Design C), requiring the width w to be downscaled by a factor of 17, and period p by a factor of 12.5. The gap size a shrinks only by a factor of 5 from 100 nm to 19 nm; the ability to fabricate increasingly small gaps will no doubt be a limiting factor in any attempt to realize mid-IR CRLH designs at 20–100 THz. We do not present such higher frequency designs here since the small dimensions result in increasingly large corrections on the graphene conductivity model due to spatial dispersion (i.e., non-local conductivity).30,31 Such effects are not readily incorporated into the finite-element eigenfrequency simulations used here.

TABLE I.

Dimensions of frequency scaled designs.

DesignCenter frequency (THz)Periodicity p (μm)Width w (μm)Gap a (nm)
2.8 2.5 4.55 100 
0.8 1.42 40 
12 0.2 0.260 19 
DesignCenter frequency (THz)Periodicity p (μm)Width w (μm)Gap a (nm)
2.8 2.5 4.55 100 
0.8 1.42 40 
12 0.2 0.260 19 
FIG. 5.

Simulated dispersion relations for the frequency scaled designs in Table I (points), along with fits to transmission line model presented in the  Appendix. Insets are zoomed in views of dispersion near β = 0.

FIG. 5.

Simulated dispersion relations for the frequency scaled designs in Table I (points), along with fits to transmission line model presented in the  Appendix. Insets are zoomed in views of dispersion near β = 0.

Close modal

In principle, graphene holds the promise of high mobility and low plasmonic dissipation—however, in practice, the plasmon damping rates depend sensitively upon the graphene quality and substrate interactions. Furthermore, for modes within the light line, radiative losses must also be considered. In this section, we consider both sources of loss and the consequences for CRLH graphene waveguides.

Radiative loss is calculated using the same finite element eigenfrequency solver described previously in Sec. III A. However, as previously discussed, the eigenfrequency solver overestimates Drude losses. Therefore, we implemented a hybrid method, where Drude losses are deliberately suppressed within the simulation by setting τ>1 ns. As a result, the numerically simulated quality factor contains only radiative losses, which we label Qr. The effect of free-carrier loss is separately estimated assuming τ = 1 ps using the analytic expression for the Drude quality factor that is derived in the  Appendix (see (A11)). The total quality factor is then given by Q1=Qr1+QD1. Both the quality factor and power loss coefficient α are plotted for the three designs in Fig. 6.

FIG. 6.

Drude, radiative, and total quality factor assuming τ = 1 ps for (a) Design A centered at 2.8 THz, (b) Design B centered at 5 THz, and (c) Design C centered at 12 THz. The shaded areas correspond to the air and silicon light cones. The insets plot the total and radiative power attenuation coefficient α. The cartoon in (a) illustrates the equivalent electric current sources for shunt and series resonances.

FIG. 6.

Drude, radiative, and total quality factor assuming τ = 1 ps for (a) Design A centered at 2.8 THz, (b) Design B centered at 5 THz, and (c) Design C centered at 12 THz. The shaded areas correspond to the air and silicon light cones. The insets plot the total and radiative power attenuation coefficient α. The cartoon in (a) illustrates the equivalent electric current sources for shunt and series resonances.

Close modal

Several trends are evident. First, the Drude losses for τ = 1 ps are very large—for Design A at 2.8 THz (αD ∼ 1500 cm−1), which leads to short attenuation lengths (α16μm). The radiative loss is of similar magnitude, which implies a radiative efficiency of 50% or better, but further reduces the propagation length. Second, there is a notable discontinuity of Qr at β = 0. This is due to the fact the CRLH waveguide is not perfectly balanced so that at β = 0 a small 20 GHz gap remains; at this point, the modes retain their distinctive series and shunt resonance character. Only the shunt resonance mode exhibits dipolar character and radiates strongly, whereas the series resonance mode is quadrupolar and is essentially a “dark” mode with a large Qr (see Fig. 6(a)). Away from β = 0, the CRLH propagating mode is a hybrid of both series and shunt resonances, which radiates only via its shunt dipolar component. In fact, the presence of this residual 20 GHz gap is likely a consequence of the fact that only the shunt-mode radiates while the series-mode is dark; closing the gap entirely to achieve a perfect balanced condition would require the series mode to radiate also, as described in Refs. 37 and 38. For designs B and C at higher frequencies, the total quality factor increases for two reasons. First, the Drude loss contribution scales as QDω. Second, the radiative loss contribution also decreases for higher frequencies (Qrω2) since the radiating dipole moment scales with the width w.

How realistic is τ = 1 ps? Within the Drude model, the momentum relaxation time is related to the dc mobility μ according to τ=μEF/evF2, so τ = 1 ps corresponds to a mobility of 15 600 cm2/V s for n=3×1013 cm−2. Dc mobilities greater than 100 000 cm2/V s have been observed at low temperature in suspended graphene,39 and graphene on hexagonal boron nitride (h-BN) underlayers has demonstrated μ60000 cm2/V s (exfoliated sample) as well as 37 000 cm2/V s (CVD grown).40,41 When fully encapsulated in h-BN, graphene was shown to have room temperature mobilities as high as 120 000 cm2/V s at densities of 5×1011 cm−2, and 40 000 cm2/V s at 3×1012 cm−2, values which correspond to effective Drude relaxation times of τ0.81 ps. Direct measurement of plasmon damping at room temperature at λ10μm has shown τ200300 fs in graphene/h-BN structures with no dependence on carrier density.8 In summary graphene with τ1 ps and perhaps even higher is certainly within experimental reach; future improvements in material quality will only help in achieving reduced losses. Values of τ0.1 ps are more typical of graphene on SiO2 substrates; such material quality is not sufficient for plasmonic CRLH waveguide operation in the THz.

A balanced CRLH based graphene structure can be considered for a leaky-wave antenna because its ability to scan its main beam from backwards to broadside to forwards.25,42 Leaky-wave antennas in graphene have been proposed via various schemes that use Bragg scattering from periodic perturbation, such as sinusoidal modulation of the nanoribbon width,21 or backgating to create a periodic reactance.18 This proposed CRLH leaky-wave antenna would differ in two significant ways. First, radiation occurs because the fundamental dispersion of the metamaterial waveguide is within the leaky-wave region—not due to Bragg scattering. As a result, CRLH structures can be balanced so that there is no (in theory) or negligibly small (in practice) stopband at β = 0, which enables propagating waves across the entire leaky-wave band. Second, the radiation of our CRLH structure is mediated through transverse dipolar currents along the graphene microribbon (i.e., through the shunt inductor). Coupling in and out of the metamaterial waveguide thereby occurs via coupling with transverse-electric (TE) plane waves, rather than transverse-magnetic (TM) plane waves, as is typical for conventional surface plasmons.27 

We study beam patterns from Design A using a driven 3D full-wave simulation, where a finite length (300 μm) CRLH waveguide is excited as a LWA. We considered the Drude relaxation time as τ = 1 ps and τ= (radiative losses only).

For the purposes of simulation, a near-field excitation is used to drive the CRLH mode, with a horizontal electrical dipole antenna placed 0.1 μm above the first unit cell as shown in Fig. 7(a). This is designed to excite the transverse current flow associated with the shunt resonance component of the CRLH mode. In practice, a CRLH LWA might be excited by replacing the horizontal dipole with a photoconductive switch for narrowband THz generation via photomixing.43 Alternately, the graphene CRLH structure might be excited using a microstrip-type waveguide—such as from a THz quantum-cascade laser—operating in its m = 1 higher-order mode.26,35 Or, incident TE plane waves can be used to excite single or arrays of multiple CRLH waveguides.27 

FIG. 7.

(a) Schematic of near-field dipole antenna excitation scheme used for driven mode beam simulations. (b) Design A far-field intensity beam pattern in x–z plane for five driving frequencies simulated with τ=. (c) Far-field intensity beam pattern in (c) y–z plane at 2.75 THz and (d) x–z plane simulated with τ = 1 ps. All beams are normalized to unity peak intensity.

FIG. 7.

(a) Schematic of near-field dipole antenna excitation scheme used for driven mode beam simulations. (b) Design A far-field intensity beam pattern in x–z plane for five driving frequencies simulated with τ=. (c) Far-field intensity beam pattern in (c) y–z plane at 2.75 THz and (d) x–z plane simulated with τ = 1 ps. All beams are normalized to unity peak intensity.

Close modal

The CRLH structure primarily radiates into the silicon half-space with a beam that scans from backwards to forwards as the frequency is changed. In practice, radiation from the substrate could be out-coupled and directed using a silicon hemispherical lens as is common for THz photoconductive antennas. If Drude loss is neglected (as seen in Fig. 7(b)), scanning is observed in the x-z plane from 150° to 225°. At 2.75 THz (near β = 0), the LWA radiates broadside with a fan-like beam. A typical beam pattern in the y-z plane is shown in Fig. 7(c). In the far-field, the electric field is predominately linearly polarized in the ϕ-direction (i.e., transverse to the LWA axis). This is consistent with our understanding of radiation originating from the transverse electric current dipoles associated with the shunt resonance component of the CRLH mode. The cross polarized component is smaller by 10 dB or more. When a Drude relaxation time of τ=1 ps is used, the attenuation length is shorter, so the beam is broader, and scanning is much less pronounced, from 180° to 210° (see Fig. 7(d)). This highlights the importance of high mobility graphene for leaky wave applications.

We now consider the effect of a perfect electric conductor backplane placed a distance h away from the graphene microribbon on the substrate side. There are two possible motivations for this: radiation control and electrostatic gating. We simulate several different backplane distances h using Design A and obtain their corresponding dispersion relationships. First, we consider h = 8 μm which corresponds to approximately λ/4 at 2.8 THz within silicon. At this distance, the ground plane is sufficiently far away such that there is little interaction between it and the fringing fields of the graphene plasmonic mode; the dispersion relation is unaffected (see Fig. 8). The radiative quality factor is also similar in magnitude to the case without the backplane. However, care must be taken when adding a backplane when the dielectric thickness is close to a quarter wavelength, as undesirable coupling with surface waves may result.44 A straightforward solution to suppress these parasitic modes is to etch the dielectric slab into a ridge underneath the graphene. Such a structure will prevent excitation of surface waves, while leaving the near-field and dispersion relation unchanged. Although not shown here, simulations show that this scheme preserves the scanning beam behavior except into the air side rather than the substrate.

FIG. 8.

Dispersion relationship and associated radiative quality factor Qr (inset) for Design A with an infinite groundplane distance h placed below the graphene strips. When the groundplane is placed close at h = 0.2 μm, Design A must be rebalanced; this new design is labeled Design A.

FIG. 8.

Dispersion relationship and associated radiative quality factor Qr (inset) for Design A with an infinite groundplane distance h placed below the graphene strips. When the groundplane is placed close at h = 0.2 μm, Design A must be rebalanced; this new design is labeled Design A.

Close modal

A smaller backplane separation can be used to reduce radiative loss as well as allow for electrostatic gating of the graphene. As seen in Fig. 8, the presence of the backplane significantly changes the dispersion relation for h=1μm as ωsh drops due to the increase in shunt capacitance CR that results from interaction with the ground plane. For the case of h=0.2μm, the change is so large that a redesign is necessary to “re-balance” the CRLH waveguide; the dimensions and dispersion is shown in Fig. 8 for this “Design A.” The width was reduced to 1.9 μm to increase ωsh back to 2.8 THz. However, since the reduction in w also reduces CL, the gap a is reduced to 25 nm to keep ωse constant. The presence of such a close backplane effectively shorts out the radiation of the CRLH structure via the image currents, and the radiative quality factor becomes very large (Qr104). This configuration will also allow tuning of the CRLH dispersion via backgating of the carrier density with modest voltages, perhaps for applications such as phase shifters or tunable resonators. Comparison with recent experimental results with similar gate dielectric thicknesses (in SiO2) suggests that voltages on the order of 50–100 V would be required to tune the Fermi level from the Dirac point to the values of 0.5–1 eV necessary for our designs.12,45 An alternative to reduce the gate voltage and obviate the need for a backgate is to use an ion-gel top-gate, which has been shown to modulate the Fermi level over 1 eV with a voltage swing of <10 V.45 

In this work, we have introduced a new design paradigm for graphene plasmonic CRLH waveguides and presented several designs in the terahertz frequency range. We have also characterized free-carrier and radiative losses, leaky-wave antenna radiation patterns, and the effects of carrier density tuning. A transmission-line model provides qualitative design guidance and quantitative agreement with numerical results with the assistance of fitting parameters. The dimensional requirements are challenging (i.e., the 2.8 THz design has ribbons of width 4.5 μm and gaps of 100 nm, and the 12 THz design has ribbons of width 260 nm and gaps of 19 nm), but within what is achievable using electron-beam lithography. However, the primary challenges for scaling to higher frequencies is the requirement that dimensions shrink quadratically with frequency as dictated by the square root dispersion relation. The gap size a is the smallest critical dimension and will likely be the primary limiting factor for fabrication if these structures were to be scaled up to the mid-IR, since it would need to shrink to just a few nanometers.

The issue of loss is the most critical challenge for practical CRLH graphene devices. For example at 3 THz, for practical values of the graphene mobility (i.e., τ = 1 ps, μ15600 cm2/V s), losses are large and propagation lengths limited. Radiative losses are also large—this benefits radiative efficiency for leaky-wave antennas but results in reduced beam scanning and large beamwidths. Hence, the most critical challenge for practical implementation is the need for large-area high-mobility graphene to minimize plasmon dissipation and achieve times of τ>1 ps. For guided wave operation, radiative losses can be eliminated with the addition of a backplane, which can also act as a backgate for free-carrier tuning. By tuning the carrier concentration by one order of magnitude, 180° of phase shift can be obtained even over the length of one unit cell—an attractive feature for a distributed phase shifter. The situation for losses improves somewhat at higher frequencies (i.e., 5 THz and 12 THz), since the Drude quality factor for CRLH modes scales proportionally to ω, and the radiative quality factor Qr scales as ω2 due to the quadratic reduction in waveguide dimensions. These challenges are not unique to our design—indeed, they are ubiquitous for any device based upon plasmons.46 However, the reduced group velocity of the CRLH waveguide reduces the propagation length by a factor two and the Q by a factor of ∼0.75 compared to an unpatterned graphene plasmon.

In summary, our results show that while demonstrating CRLH graphene waveguides is possible, it will be challenging, with the ultimate performance highly dependent upon availability of high-mobility graphene with long Drude relaxation times. This challenge will need to be weighed against potential benefits of using graphene, such as free-carrier tunability, or large optical nonlinearity.

The authors would like to thank Ziwen Wang for his preliminary simulation work. This research was supported in part by National Science Foundation Grant Nos. ECCS 1150071 and ECCS 1340928.

The theory of CRLH metamaterials is most easily presented in terms of a transmission-line model. In the graphene microribbon waveguide presented here, this is not strictly valid since the CRLH behavior is based upon a higher-order transverse mode, which results in an approximately sinusoidal variation in the electric and magnetic fields across the width of the ribbon. Thus, the waveguide cannot be considered sub-wavelength in the transverse dimension, and the quasi-static limit is not applicable. Nonetheless, a circuit model is still very useful for qualitative understanding and design guidance, provided that it is understood that inductances and capacitances are effective parameters. However, for detailed design, numerical methods must be used. We begin by defining the series inductance LR and shunt capacitance CR for a graphene microribbon with an effective width weff and unit cell length p

LR=pweffωIm{σ1}=pπ2weffe2EF,
(A1)
CR=DRpweff(ε1+ε2)Re{q(ω)},
(A2)
=DRpweffπ2(ε1+ε2)2e2EFω2.
(A3)

The series inductance LR originates from the kinetic inductance of the graphene free carriers, i.e., the imaginary part of the Drude resistivity in (2). We neglect any additional “Faraday inductance” associated with the graphene currents, which is negligible in the non-retarded limit (i.e., qω/c).31 The effective shunt capacitance CR was obtained by equating the dispersion relation of a conventional transmission line β=ωLRCRp to the graphene plasmon dispersion relation given in (4) for an infinite graphene sheet. This requires CR to acquire a frequency/wavenumber dependence. However, an additional factor DR has been added which enhances the capacitance to account for the modification of the dispersion due to the finite width of the microribbon. DR can be obtained either through matching to numerical simulations, or by using the term η calculated in Ref. 34 (DR(2π2η)2 in which η is evaluated at qw=π). Using this expression, we accordingly take DR = 1.8, 1.8, and 1.9 for Designs A, B, and C, respectively. This expression for CR also neglects the effects of quantum-capacitance, which is appropriate here where the non-local corrections to conductivity are negligible.30,31 For the fundamental lateral mode m = 0, the effective width weff is equal to the physical ribbon width w. In this case, these expressions are equivalent to those derived in Ref. 31, taken in the limit that kBTEF, and ignoring the negligible contributions of quantum capacitance, and Faraday inductance.

The CRLH waveguide mode is based upon an m = 1 higher order lateral mode (i.e., an edge mode) and includes series capacitive gaps of size a. We model this according to the double transmission line circuit model shown in Fig. 2(d) operating in its odd mode. The previous expressions for LR and CR still apply, provided that the effective width in this case is weff=Fw/2. This reflects the fact that branch of the TL has physical width of w/2, and the effective width is modified by a factor F to reflect that the current, field, and charge are nonuniform over the width of the ribbon. We can approximate this factor as F=2/πA; ideally A = 1, which is appropriate for a field distribution takes on exactly a half-sinusoidal transverse variation.27 In reality, there is a slight phase shift associated with the reflection of the plasmon from the microribbon edge; we find A = 0.9, 0.96, and 1.0 for designs A, B, and C, respectively. We also define the shunt inductance

LL=weffpωIm{σ1}=weffπ2pe2EF,
(A4)

which is based upon the same kinetic inductance of the graphene used to define LR, except now for a section of length of weff and width of p. This ensures that the effective shunt frequency ωsh=(LLCR)1/2 corresponds the condition where the microribbon width w is nearly equal to one-half of the plasmon wavelength such that q(ωsh)=Aπ/w.

A working expression for the series capacitance CL can be obtained from the capacitance between metal patches of width weff and length p/2 on a dielectric half-space in the limit that pa (derived for a Sievenpiper surface in Ref. 47), which gives

CL=DLweffπ(ε1+ε2)ln(2p/a).
(A5)

DL is a phenomenological factor that is required to reduce ωse=(LRCL)1/2 so that it matches the simulated result. We fit the values of DL = 0.7, 0.66, and 0.6 for the dimensions of designs A, B, and C, respectively. Without this factor, CL is overestimated since the period of the structure is not sufficiently subwavelength to accurately use a lumped element expression.

Despite the presence of the phenomenological parameters, we can still draw important conclusions. For example, by equating these analytic expressions for the series and shunt resonance frequencies, we are able to obtain a transcendental equation for the CRLH balanced condition

DLln(2p/a)=wpDRA.
(A6)

This expression for the balanced condition contains only geometric parameters; there is no dependence on carrier concentration or other material properties. This implies that the balanced dispersion relation will be preserved as the carrier density is changed, which is consistent with the tuning observed in the numerical simulations in Fig. 4.

Drude loss within the graphene can be included by adding resistances

RR=2pπ2weffe2EFτ,RL=2weffπ2pe2EFτ,
(A7)

in series with LR and LL, respectively. We can further use this formalism to derive an analytic expression for loss in a CRLH waveguide. The dispersion for a balanced line is given25 

β=ωp(LR+iRR/ω)CR1pω1(LL+iRL/ω)CL.
(A8)

This gives a power attenuation coefficient

α=2Im{β}=2pτLRCR+2pω2τ1LLCL.
(A9)

For a balanced CRLH structure, at β = 0, the attenuation coefficient is

α(ω=ωsh=ωse)=4pτLRCR=4π2(ε1+ε2)ωτe2EF,
(A10)

which is twice as large as for an unpatterned m = 0 graphene plasmon mode at the same frequency. The quality factor associated with Drude loss can then be obtained using the relation QD=ω/αvg; at β = 0, it has the value of

QD=ω/αvg.
(A11)

This relation is used to estimate the Drude losses in Fig. 6. We note that at β = 0, QD=(3/4)ωτ—three quarters of the value of QD of an unpatterned m = 0 graphene plasmon.

1.
M.
Jablan
,
H.
Buljan
, and
M.
Soljačić
,
Phys. Rev. B
80
,
245435
(
2009
).
2.
E. H.
Hwang
and
S.
Das Sarma
,
Phys. Rev. B
75
,
205418
(
2007
).
3.
F. H. L.
Koppens
,
D. E.
Chang
, and
F. J.
García de Abajo
,
Nano Lett.
11
,
3370
(
2011
).
4.
A. Y.
Nikitin
,
F.
Guinea
,
F. J.
García-Vidal
, and
L.
Martín-Moreno
,
Phys. Rev. B
84
,
161407
(
2011
).
5.
A. N.
Grigorenko
,
M.
Polini
, and
K. S.
Novoselov
,
Nat. Photonics
6
,
749
(
2012
).
6.
Z.
Fei
,
G. O.
Andreev
,
W.
Bao
,
L. M.
Zhang
,
A. S.
McLeod
,
C.
Wang
,
M. K.
Stewart
,
Z.
Zhao
,
G.
Dominguez
,
M.
Thiemens
,
M. M.
Fogler
,
M. J.
Tauber
,
A. H.
Castro-Neto
,
C. N.
Lau
,
F.
Keilmann
, and
D. N.
Basov
,
Nano Lett.
11
,
4701
(
2011
).
7.
Z.
Fei
,
A. S.
Rodin
,
G. O.
Andreev
,
W.
Bao
,
A. S.
McLeod
,
M.
Wagner
,
L. M.
Zhang
,
Z.
Zhao
,
M.
Thiemens
,
G.
Dominguez
,
M. M.
Fogler
,
A. H. C.
Neto
,
C. N.
Lau
,
F.
Keilmann
, and
D. N.
Basov
,
Nature
487
,
82
85
(
2012
).
8.
A.
Woessner
,
M. B.
Lundeberg
,
Y.
Gao
,
A.
Principi
,
P.
Alonso-González
,
M.
Carrega
,
K.
Watanabe
,
T.
Taniguchi
,
G.
Vignale
,
M.
Polini
,
J.
Hone
,
R.
Hillenbrand
, and
F. H. L.
Koppens
,
Nat. Mater.
14
,
421
(
2014
).
9.
J.
Chen
,
M.
Badioli
,
P.
Alonso-González
,
S.
Thongrattanasiri
,
F.
Huth
,
J.
Osmond
,
M.
Spasenović
,
A.
Centeno
,
A.
Pesquera
,
P.
Godignon
,
A. Z.
Elorza
,
N.
Camara
,
F. J. G.
de Abajo
,
R.
Hillenbrand
, and
F. H. L.
Koppens
,
Nature
487
,
77
81
(
2012
).
10.
L.
Ju
,
B.
Geng
,
J.
Horng
,
C.
Girit
,
M.
Martin
,
Z.
Hao
,
H. A.
Bechtel
,
X.
Liang
,
A.
Zettl
,
Y. R.
Shen
, and
F.
Wang
,
Nat. Nanotechnol.
6
,
630
(
2011
).
11.
H.
Yan
,
X.
Li
,
B.
Chandra
,
G.
Tulevski
,
Y.
Wu
,
M.
Freitag
,
W.
Zhu
,
P.
Avouris
, and
F.
Xia
,
Nat. Nanotechnol.
7
,
330
(
2012
).
12.
V. W.
Brar
,
M. S.
Jang
,
M.
Sherrott
,
J. J.
Lopez
, and
H. A.
Atwater
,
Nano Lett.
13
,
2541
(
2013
).
13.
H.
Yan
,
T.
Low
,
W.
Zhu
,
Y.
Wu
,
M.
Freitag
,
X.
Li
,
F.
Guinea
,
P.
Avouris
, and
F.
Xia
,
Nat. Photonics
7
,
394
(
2013
).
14.
X.
Zhu
,
W.
Wang
,
W.
Yan
,
M. B.
Larsen
,
P.
Bøggild
,
T. G.
Pedersen
,
S.
Xiao
,
J.
Zi
, and
N. A.
Mortensen
,
Nano Lett.
14
,
2907
(
2014
).
15.
R.
Yan
,
B.
Sensale-Rodriguez
,
L.
Liu
,
D.
Jena
, and
H. G.
Xing
,
Opt. Express
20
,
28664
(
2012
).
16.
B.
Sensale-Rodriguez
,
R.
Yan
,
M. M.
Kelly
,
T.
Fang
,
K.
Tahy
,
W. S.
Hwang
,
D.
Jena
,
L.
Liu
, and
H. G.
Xing
,
Nat. Commun.
3
,
780
(
2012
).
17.
S. H.
Lee
,
M.
Choi
,
T.-T.
Kim
,
S.
Lee
,
M.
Liu
,
X.
Yin
,
H. K.
Choi
,
S. S.
Lee
,
C.-G.
Choi
,
S.-Y.
Choi
,
X.
Zhang
, and
B.
Min
,
Nat. Mater.
11
,
936
(
2012
).
18.
M.
Esquius-Morote
,
J. S.
Gómez-Díaz
, and
J.
Perruisseau-Carrier
,
IEEE Trans. Terahertz Sci. Technol.
4
,
116
(
2014
).
19.
R.
Filter
,
M.
Farhat
,
M.
Steglich
,
R.
Alaee
,
C.
Rockstuhl
, and
F.
Lederer
,
Opt. Express
21
,
3737
(
2013
).
20.
I.
Llatser
,
C.
Kremers
,
A.
Cabellos-Aparicio
,
E. A.
Josep Miquel Jornet
, and
D. N.
Chigrin
,
Photonics Nanostruct.: Fundam. Appl.
10
,
353
(
2012
).
21.
J. S.
Gómez-Díaz
,
M.
Esquius-Morote
, and
J.
Perruisseau-Carrier
,
Opt. Express
21
,
24856
(
2013
).
22.
C.
Caloz
and
T.
Itoh
, in
Proceedings of the IEEE-AP-S USNC/URSI National Radio Science Meeting
(
2002
), Vol.
2
, p.
412
.
23.
A.
Lai
,
C.
Caloz
, and
T.
Itoh
,
IEEE Microwave Mag.
5
,
34
(
2004
).
24.
G. V.
Eleftheriades
,
A. K.
Iyer
, and
P. C.
Kremer
,
IEEE Trans. Microwave Theory Tech.
50
,
2702
(
2002
).
25.
C.
Caloz
and
T.
Itoh
,
Electromagnetic Metamaterials: Transmission Line Theory and Microwave Applications
(
Wiley
,
2006
).
26.
A. A.
Tavallaee
,
P. W. C.
Hon
,
Q.-S.
Chen
,
T.
Itoh
, and
B. S.
Williams
,
Appl. Phys. Lett.
102
,
021103
(
2013
).
27.
P. W. C.
Hon
,
Z.
Liu
,
T.
Itoh
, and
B. S.
Williams
,
J. Appl. Phys.
113
,
033105
(
2013
).
28.
J.
Horng
,
C.-F.
Chen
,
B.
Geng
,
C.
Girit
,
Y.
Zhang
,
Z.
Hao
,
H.
Bechtel
,
M.
Martin
,
A.
Zettl
,
M.
Crommie
,
Y.
Shen
, and
F.
Wang
,
Phys. Rev. B
83
,
165113
(
2011
).
29.
F. J.
García de Abajo
,
ACS Photonics
1
,
135
(
2014
).
30.
G.
Lovat
,
G. W.
Hanson
,
R.
Araneo
, and
P.
Burghignoli
,
Phys. Rev. B
87
,
115429
(
2013
).
31.
D.
Correas-Serrano
,
J. S.
Gómez-Díaz
,
J.
Perruisseau-Carrier
, and
A.
Álvarez Melcón
,
IEEE Trans. Terahertz Sci. Technol.
61
,
4333
(
2013
).
32.
S.
Thongrattanasiri
,
A.
Manjavacas
, and
F. J.
García de Abajo
,
ACS Nano
6
,
1766
(
2012
).
33.
A.
Fallahi
,
T.
Low
,
M.
Tamagnone
, and
J.
Perruisseau-Carrier
,
Phys. Rev. B
91
,
121405(R)
(
2015
).
34.
J.
Christensen
,
A.
Manjavacas
,
S.
Thongrattanasiri
,
F. H. L.
Koppens
, and
F. J.
García de Abajo
,
ACS Nano
6
,
431
(
2012
).
35.
A. A.
Tavallaee
,
B. S.
Williams
,
P. W. C.
Hon
,
T.
Itoh
, and
Q.-S.
Chen
,
Appl. Phys. Lett.
99
,
141115
(
2011
).
36.
J.
Gomez-Diaz
and
J.
Perruisseau-Carrier
, in
2012 International Symposium on Antennas and Propagation (ISAP)
(
2012
), pp.
239
242
.
37.
J. S.
Gómez-Díaz
,
D.
Canete-Rebenaque
, and
A.
Álvarez Melcón
,
IEEE Antennas Wireless Propag. Lett.
10
,
29
(
2011
).
38.
S.
Otto
,
A.
Rennings
,
K.
Solbach
, and
C.
Caloz
,
IEEE Trans. Antennas Propag.
59
,
3695
(
2011
).
39.
K. I.
Bolotin
,
K. J.
Sikes
,
Z.
Jiang
,
M.
Klima
,
G.
Fudenberg
,
J.
Hone
,
P.
Kim
, and
H. L.
Stormer
,
Solid State Commun.
146
,
351
(
2008
).
40.
W.
Gannett
,
W.
Regan
,
K.
Watanabe
,
T.
Taniguchi
,
M. F.
Crommie
, and
A.
Zettl
,
Appl. Phys. Lett.
98
,
242105
(
2011
).
41.
C. R.
Dean
,
A. F.
Young
,
I.
Meric
,
C.
Lee
,
L.
Wang
,
S.
Sorgenfrei
,
K.
Watanabe
,
T.
Taniguchi
,
P.
Kim
,
K. L.
Shepard
, and
J.
Hone
,
Nat. Nanotechnol.
5
,
722
(
2010
).
42.
S.
Lim
,
C.
Caloz
, and
T.
Itoh
,
IEEE Trans. Microwave Theory Tech.
53
,
161
(
2005
).
43.
E. R.
Brown
,
K. A.
McIntosh
,
K. B.
Nichols
, and
C. L.
Dennis
,
Appl. Phys. Lett.
66
,
285
(
1995
).
44.
A.
Shahvarpour
,
A.
Álvarez Melcón
, and
C.
Caloz
,
IEEE Trans. Antennas Propag.
61
,
4013
(
2013
).
45.
H.
Hu
,
F.
Zhai
,
D.
Hu
,
Z.
Li
,
B.
Bai
,
X.
Yang
, and
Q.
Dai
,
Nanoscale
7
,
19493
(
2015
).
46.
J. B.
Khurgin
,
Nat. Nanotechnol.
10
,
2
(
2015
).
47.
D. F.
Sievenpiper
, “
High-impedance electromagnetic surfaces
,” Ph.D. thesis,
Institution is University of California Los Angeles
,
1999
.