We present an experimental and numerical study of a terahertz metamaterial with a nonlinear response that is controllable via the relative structural arrangement of two stacked split ring resonator arrays. The first array is fabricated on an n-doped GaAs substrate, and the second array is fabricated vertically above the first using a polyimide spacer layer. Due to GaAs carrier dynamics, the on-resonance terahertz transmission at 0.4 THz varies in a nonlinear manner with incident terahertz power. The second resonator layer dampens this nonlinear response. In samples where the two layers are aligned, the resonance disappears, and the total nonlinear modulation of the on-resonance transmission decreases. The nonlinear modulation is restored in samples where an alignment offset is imposed between the two resonator arrays. Structurally tunable metamaterials and metasurfaces can therefore act as a design template for tunable nonlinear THz devices by controlling the coupling of confined electric fields to nonlinear phenomena in a complex material substrate or inclusion.

Over the past decade, the field of terahertz (THz) science has grown to include the study and engineering of nonlinear optical materials.1,2 This growth in nonlinear THz science has been driven by improvements in the generation of high power THz pulses in the laboratory. Standard THz generation via photo-conductive antennas or optical rectification (e.g., in ZnTe) produces THz pulses limited to femtojoules of energy.3 By contrast, the newer techniques of tilted pulse front THz generation (TPFG) in LiNbO3 and 4-wave mixing in air generate pulse energies orders of magnitude higher.4,5 In particular, TPFG methods are known to reliably produce THz pulses with microjoules of energy, corresponding to peak electric field strengths on the order of 300–500 kV/cm, and in some instances up to 1 MV/cm.6 

The use of high peak power THz pulses makes time resolved nonlinear experiments possible at THz frequencies. Researchers are able to study the THz nonlinear properties of materials on ultrafast timescales7 and design novel nonlinear optical elements for THz radiation.8,9 Combining THz excitation with x-ray probes allows for time-resolved measurements of nonlinear structural changes in complex materials.10 

Within the field of nonlinear THz science, plasmonic metamaterials and metasurfaces (MMs) play a critical role. MMs are engineered composites with optical properties that are determined by the geometry and layout of the component’s sub-wavelength inclusions.11 MMs have been applied across the electromagnetic spectrum and may be used to precisely control light through manipulation of phase, intensity, and polarization.12 This unprecedented level of control over light has led to notable MM demonstrations of left-handed materials,13 electromagnetic cloaking,14 and perfect absorption.15,16 In addition, MMs are interesting from a materials science standpoint due to their enhancement of light–matter interaction.2,7 This interaction can be seen in a resonant nonlinear MM response, making MMs an ideal platform to study THz nonlinear dynamics and design nonlinear devices, such as saturable absorbers.8 

Nonlinear MM devices have been created at microwave and terahertz frequencies through a variety of methods, including incorporating nonlinear lumped circuit elements into the design,17–19 and by using the inherent nonlinear response of the subwavelength inclusions that make up the MM.20 Another approach used in nonlinear MM devices makes use of resonant inclusions, such as the split ring resonator (SRR), which confine electric fields to localized regions in the unit cell. Resonant field confinement (FC) enhances the local peak electric field intensity and excites a nonlinear response in the MM inclusions or substrate.7,21 Many demonstrations of nonlinear MMs based on FC exist, including devices that use confined fields to enhance charge carrier nonlinearities in VO2, GaAs, InAs, and graphene21,22,23 or to excite phase transitions in both low and high Tc superconductors.24–26 

Dynamic control over the MM response dramatically extends device capability, and integration of control and tuning mechanisms is an active area of study for linear and nonlinear MM devices. Two of the most commonly used methods for electrical control of MM properties are via modulation of substrate conductivity in, for instance, GaAs27 or graphene28 and structural actuation, for instance, via microelectromechanical systems (MEMS) based devices.29,30 Manipulating multi-component unit cells is of large importance for MEMS tuning mechanisms. However, to date, there has been no report on how the enhanced nonlinear charge carrier response of an MM substrate responds to structural arrangement of a multi-component unit cell.

In this paper, we report an experimental, proof of principle result showing how structural changes to a multi-resonator unit cell can be used to control the enhanced nonlinear response of the substrate in a THz MM. Our MM design is based on the common broadside-coupled SRR (BCSRR)31 and is shown schematically in Fig. 1(a). This design is composed of two layered SRR arrays, oriented to maximize electromagnetic interactions between the two arrays. More details on device design are discussed below. We show that by altering the lateral positions of the two component split ring resonators in the unit cell, one can control the FC in the capacitive gaps of the resonators and thus control the coupling of incident THz fields to the nonlinear carrier dynamics of an n-doped GaAs substrate. When the two resonator layers are aligned to overlap as shown in Fig. 1(b), very little nonlinear behavior is seen in the device response at a design frequency of 0.4 THz. When a lateral offset is placed between the upper and lower SRR array [shown in Fig. 1(c)], the FC in the capacitive gap of the lower resonator increases, enabling a nonlinear response from the carrier dynamics of the n-doped epilayer at a relatively low incident THz power. The nonlinear carrier dynamics result in a drop in the epilayer conductivity, turning on the lowest order SRR resonance at 0.4 THz. The resonance is thus highly dependent on incident THz power, and the on-resonance THz transmission is modulated by ∼3 dB as the incident THz power is increased. Below, we present nonlinear THz spectroscopy measurements to characterize device response and numerical simulations to provide insight into the physical mechanisms for the device behavior.

FIG. 1.

(a) Vertically expanded perspective and the side view of the MM unit cell showing dimensions and the direction of THz excitation. (b) and (c) Mid-fabrication photos of two samples showing (b) 0 and (c) 48 μm lateral offset between the stacked resonator arrays. The lower layer of SRRs is shown in gold. The upper layer is shown as outlines.

FIG. 1.

(a) Vertically expanded perspective and the side view of the MM unit cell showing dimensions and the direction of THz excitation. (b) and (c) Mid-fabrication photos of two samples showing (b) 0 and (c) 48 μm lateral offset between the stacked resonator arrays. The lower layer of SRRs is shown in gold. The upper layer is shown as outlines.

Close modal

The MM studied in this work consists of two stacked, planar arrays of gold SRRs as shown in Fig. 1. The lower array is fabricated directly onto a GaAs substrate with a 1 μm thick n-doped GaAs epilayer (n = 2 × 1016 cm−3) using standard photolithographic techniques. The SRRs in this layer are inter-connected with metallic wires for use in applying an electrical bias (not used in this study). A 2 μm thick polyimide spacer layer is then deposited on the epilayer, followed by the second SRR array. The orientation of the second array is rotated by 180° relative to the first array, which maximizes electromagnetic coupling between the two layers.32,33

The dimensions of the component SRRs are chosen to maximize coupling between the two resonator layers and to place the device resonance near the peak of the TPFG THz signal at ∼0.4 THz. The key dimensions of each resonator are marked on the upper SRR in Fig. 1(a). These are the side length, L, the capacitive gap width, g, the linewidth, w, and the square array periodicity, P. Below, the dimensions for the lower SRR array are written with a subscript 1, while dimensions for the upper array are written with a subscript 2. In the lower array, L1 = 28 μm, linewidth w1 = 6 μm, and g1 = 2 μm. In the upper array, L2 = 48 μm, w2 = 6 μm, and g2 = 2 μm. The unit cell periodicity is the same for both arrays and is P = 96 μm. Two samples are fabricated, each with a varying lateral shift, Lshift, between the centers of the two SRRs, as shown in the side view of Fig. 1(a). Photographs of the two fabricated samples are shown in Figs. 1(b) and 1(c). In sample 1 [Fig. 1(b)], the two SRR arrays are directly aligned with Lshift = 0 μm. In sample 2 [Fig. 1(c)], the two arrays are laterally offset by half of the unit cell periodicity (Lshift = 48 μm).

To characterize the nonlinear response of the two samples, THz pulses with field strengths between 50 and 400 kV/cm were generated using TPFG in LiNO36 and focused onto the MM at normal incidence, with the THz electric field polarized perpendicular to the SRR capacitive gaps, as shown in Fig. 1(a). The transmitted pulses are then measured using electro-optic sampling in ZnTe in a standard THz time-domain spectroscopy (THz-TDS) configuration. The transmission spectra are obtained through Fourier transform and normalized to a THz-TDS reference measurement of a bare GaAs substrate with a 1 μm n-doped epilayer.

Experimental transmission spectra as a function of field strength for both samples are shown in Fig. 2. For the sample with Lshift = 0 μm [Fig. 2(a)], only a small nonlinear response can be seen in the data. As the THz field strength is increased from 50 to 400 kV/cm, the overall transmission through the sample slightly increases, but no clear resonance is seen, regardless of incident field strength. Figure 2(b) shows noticeably different behavior exhibited by the structure in which the SRRs are offset. The MM now exhibits a strong resonance near ∼0.4 THz with a larger nonlinear modulation. As the THz field strength is increased from 50 to 400 kV/cm, the on-resonance transmission is modulated by ∼3 dB.

FIG. 2.

THz-TDS data for samples with (a) Lshift = 0 and (b) Lshift = 48 μm. Insets show photographs of the relative positioning for the two SRR layers for the data shown in the plot. On each plot, the dashed vertical lines mark the position of the lowest order resonance of the BC-SRR system.

FIG. 2.

THz-TDS data for samples with (a) Lshift = 0 and (b) Lshift = 48 μm. Insets show photographs of the relative positioning for the two SRR layers for the data shown in the plot. On each plot, the dashed vertical lines mark the position of the lowest order resonance of the BC-SRR system.

Close modal

Thus, the presence of a strong nonlinear modulation in the MM resonance can be controlled solely by structural positioning of the two MM layers. As another comparison of the stark difference in the modulation range, Fig. 3 plots the transmission at 0.4 THz for both samples as a function of incident THz field strength. Not only is the modulation range noticeably increased in the case for Lshift = 48 μm, the direction of modulation is also reversed compared to the unshifted structure.

FIG. 3.

Modulation of the transmission minimum at 0.4 THz vs applied THz field strength for both shifted and unshifted BC-SRR samples. Cross markers on lines signify measured data points.

FIG. 3.

Modulation of the transmission minimum at 0.4 THz vs applied THz field strength for both shifted and unshifted BC-SRR samples. Cross markers on lines signify measured data points.

Close modal

Using numerical simulations, we investigate the physical origins of the nonlinear response of the BC-SRR MM discussed above. In order to show the mode structure of the MM, we simulate the THz transmission spectra of the BC-SRR structure for Lshift = 0 and 48 μm using commercial solvers based on finite difference time domain techniques.34 The SRR gold patterning is modeled as a lossy metal, while the GaAs substrate is modeled with a 1 µm lossy semiconductor epilayer on a loss-free semi-insulating GaAs substrate. These simulated spectra are shown together in Fig. 4(a). Due to the limited frequency resolution in the tilted pulse front THz-TDS system, the resonances in simulation are narrower than those in the experimental results.

FIG. 4.

(a) Simulated transmission spectra for the Lshift = 0 and 48 μm samples assuming no loss in the GaAs substrate. Letters A and B denote resonances discussed in the main text. (b) Simulated nonlinear response for Lshift = 0 μm. (c) Simulated nonlinear response for Lshift = 48 μm.

FIG. 4.

(a) Simulated transmission spectra for the Lshift = 0 and 48 μm samples assuming no loss in the GaAs substrate. Letters A and B denote resonances discussed in the main text. (b) Simulated nonlinear response for Lshift = 0 μm. (c) Simulated nonlinear response for Lshift = 48 μm.

Close modal

The blue curve in Fig. 4(a) shows the spectra for the sample with Lshift = 0 μm. The incident THz electric field, polarized perpendicular to the capacitive gap of the SRRs, excites two modes. The resonance at 0.73 THz (resonance A) corresponds to the electrically active coupled mode of the BC-SRR.33 The frequency position and oscillator strength of mode A have been shown in a previous work to be highly dependent on the electromagnetic coupling between the two component SRRs. The frequency position of mode A can be approximately described using an LC oscillator model,

(1)

where L is the BC-SRR total inductance and C is the BC-SRR total capacitance.35,36

The red curve in Fig. 4(a) shows the transmission spectrum for the sample with Lshift = 48 μm. The lateral offset of the two SRRs alters the mutual capacitance and inductance of the structure, red-shifting resonance A from 0.73 to 0.47 THz. The mode at 0.66 THz (resonance B) is a surface lattice mode, common in periodic MM structures.37 The frequency position of mode B is largely independent of the mutual interactions between the two SRRs, as expected for a surface lattice resonance.

Figures 4(b) and 4(c) show simulations of the nonlinear response of the Lshift = 0 and Lshift = 48 μm samples, respectively. Here, the GaAs response is modeled using the Drude model with a carrier mobility that can vary between 100 and 1000 cm2/V s. The nonlinear response arises from terahertz field induced intervalley scattering of charge carriers in the 1 μm n-GaAs epilayer. The resulting decrease in mobility reduces the conductivity of the epilayer, resulting in a terahertz field dependent increase in resonance depth for resonances A and B for both Lshift = 0 and Lshift = 48 μm.

In the experiment, the overall loss in the n-doped epilayer broadens all resonances. The resolution limit of the THz-TDS experiment also artificially broadens the observed resonances. For Lshift = 0 μm, neither resonance is seen in the experiment [Fig. 2(a)] for any value of incident THz field strength due to the above-mentioned broadening effects. However, the net decrease in carrier mobility is still seen in the data as a broadband increase in transmission.

For Lshift = 48 μm [Fig. 2(b)], resonance A is clearly visible in the experimental results since the overall oscillator strength, and thus resonance depth, of resonance A is now much greater. Here, the decrease in carrier mobility decreases the on-resonance loss of the BC-SRR structure, leading to a stronger resonance and explaining the difference in the modulation direction discussed in Fig. 3.

The lateral positioning of the two SRR layers allows for tuning of the oscillator strength by controlling the local FC within the MM unit cell. For mode A, the FC within the capacitive gap region is directly proportional to the oscillator strength of mode A. As explained below, the oscillator strength of the resonance depends on the lateral positioning of the two resonators.

The lowest order resonant current distribution of the component SRRs has been calculated in the literature.36 The resonant current is zero at the ends of the capacitive gap and maximum on the side of the SRR opposite the capacitive gap. This is similar to the current distribution of a folded half-wave dipole antenna, sometimes called a halo antenna. In this analogy, the length, d, of the dipole is the circumference of the SRR. Thus, the current in an SRR can be expressed approximately using the cosine distribution for the current of a half-wave dipole antenna,

(2)

where Io is the peak current on resonance, s is the position along the SRR circumference, and ω is the angular frequency of incident THz excitation. Although the connection between the folded dipole antenna and the SRR assumes a circular SRR, the above-mentioned analysis is still a very good approximation for the square SRRs used in this study. As the resonant current oscillates, charge conservation forces a charge build up along the capacitive gap edges of both SRRs. This charge build up oscillates 90° out of phase with the resonant current. Specifically, the charge buildup on one end of the SRR gap can be expressed as

(3)

where Qo is the peak charge at the edge of the SRR capacitive gap. The opposite edge of the capacitive gap builds up an equal amount of charge of the opposite sign. Thus, a dipole charge distribution builds up across the capacitive gap in both resonators. In the BC-SRR resonance studied in this paper, currents oscillate around both rings in opposite directions.35 Thus, the dipole charge distributions in the two resonators have opposing orientations.

With no lateral offset, the electric fields from the two electric dipoles in the upper and lower SRR superimpose and destructively interfere, leading to a low net electric field strength in the unit cell of the MM. As the two resonators are laterally offset, the resonant electric fields of the two charge distributions no longer overlap closely, leading to less destructive interference and a higher net electric field strength inside the unit cell. The end result is an increase in the strength of the local electric fields in the lower SRR gap region. This results in a stronger resonance oscillator strength for the sample with Lshift = 48 μm and a larger nonlinear modulation of the resonance as the incident THz field strength is increased.

We can confirm the above-mentioned explanation by performing simulations of the local electric field distributions within the MM unit cell. Figure 5 shows the time domain electric field strength maximum in the plane of the lower SRR for the Lshift = 0 μm [Fig. 5(a)] and Lshift = 48 μm [Fig. 5(b)] samples. The local electric fields are higher by close to a factor of 4 in the Lshift = 48 μm sample. As the local electric field amplification increases, especially in the vicinity of the SRR capacitive gap, so will the overall oscillator strength and resonance depth of the BC-SRR resonance.

FIG. 5.

Local electric field distribution in lower SRR during THz excitation for (a) Lshift = 0 and (b) Lshift = 48 μm. Insets show the unit cell alignment of the upper and lower SRRs for both plots.

FIG. 5.

Local electric field distribution in lower SRR during THz excitation for (a) Lshift = 0 and (b) Lshift = 48 μm. Insets show the unit cell alignment of the upper and lower SRRs for both plots.

Close modal

Figures 6(a) and 6(b) show simulations of the increase in resonance depth and thus oscillator strength, of resonance A, for varying values of Lshift. As Lshift is increased from 0 to 48 μm, the depth of resonance A increases by 80%, corresponding to an 80% increase to oscillator strength. Consequently, resonance A and its nonlinear modulation are visible in the experimental spectra for Lshift = 48 μm but not for Lshift = 0 μm.

FIG. 6.

(a) Simulated transmission spectra for BC-SRR samples with varying values of Lshift. (b) Change in the on-resonance transmission minimum vs Lshift. Crosses mark measured data points.

FIG. 6.

(a) Simulated transmission spectra for BC-SRR samples with varying values of Lshift. (b) Change in the on-resonance transmission minimum vs Lshift. Crosses mark measured data points.

Close modal

In conclusion, we presented a proof-of-principle study showing how a nonlinear metamaterial response can be tuned in magnitude via the relative lateral positioning of the stacked resonator arrays inside a broadside coupled split ring resonator metamaterial. The metamaterial was patterned on an n-doped GaAs substrate. We investigated the behavior of the metamaterial experimentally via terahertz time domain spectroscopy and used numerical simulations to provide physical insight into the device response. In samples where the two resonator arrays were aligned, the device showed only a small nonlinear response in the experiment. In samples where the component split ring resonator arrays were laterally shifted by 48 μm, a resonance appeared at ∼0.41 THz. Due to the charge carrier dynamics of the n-doped GaAs and on-resonance field confinement in the SRR capacitive gaps, the on-resonance THz transmission was strongly nonlinear, decreasing by ∼3 dB as the incident terahertz power increases from 50 to 400 kV/cm. This result is, to the best of our knowledge, the first example illustrating the control of the nonlinear response of a THz MM device solely through the structural positioning of the component inclusions.

This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, and performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multi-mission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA0003525. Work at the UCSD was supported by the DARPA DRINQS program (Grant No. D18AC00014). G.R.K. thanks the Washington College Faculty Enhancement Program for travel and equipment funding throughout the duration of this project.

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
F.
Blanchard
,
L.
Razzari
,
F. H.
Su
,
G.
Sharma
,
R.
Morandotti
,
T.
Ozaki
,
M.
Reid
, and
F. A.
Hegmann
, in
Nonlinear Photonics and Novel Optical Phenomena
, edited by
Z.
Chen
and
R.
Morandotti
(
Springer New York
,
2012
), Vol. 170, pp.
297
323
.
2.
G.
Keiser
and
P.
Klarskov
,
Photonics
6
(
1
),
22
(
2019
).
3.
D. M.
Mittleman
,
J. Appl. Phys.
122
(
23
),
230901
(
2017
).
4.
J. A.
Fülöp
,
L.
Pálfalvi
,
S.
Klingebiel
,
G.
Almási
,
F.
Krausz
,
S.
Karsch
, and
J.
Hebling
,
Opt. Lett.
37
(
4
),
557
559
(
2012
).
5.
D. J.
Cook
and
R. M.
Hochstrasser
,
Opt. Lett.
25
(
16
),
1210
1212
(
2000
).
6.
H.
Hirori
,
A.
Doi
,
F.
Blanchard
, and
K.
Tanaka
,
Appl. Phys. Lett.
98
(
9
),
091106
(
2011
).
7.
M.
Liu
,
H. Y.
Hwang
,
H.
Tao
,
A. C.
Strikwerda
,
K.
Fan
,
G. R.
Keiser
,
A. J.
Sternbach
,
K. G.
West
,
S.
Kittiwatanakul
,
J.
Lu
,
S. A.
Wolf
,
F. G.
Omenetto
,
X.
Zhang
,
K. A.
Nelson
, and
R. D.
Averitt
,
Nature
487
(
7407
),
345
348
(
2012
).
8.
G. R.
Keiser
,
J.
Zhang
,
X.
Zhao
,
X.
Zhang
, and
R. D.
Averitt
,
J. Opt. Soc. Am. B
33
(
12
),
2649
2655
(
2016
).
9.
G. R.
Keiser
,
N.
Karl
,
P. Q.
Liu
,
C.
Tulloss
,
H.-T.
Chen
,
A. J.
Taylor
,
I.
Brener
,
J. L.
Reno
, and
D. M.
Mittleman
,
Appl. Phys. Lett.
111
(
12
),
121101
(
2017
).
10.
A. X.
Gray
,
M. C.
Hoffmann
,
J.
Jeong
,
N. P.
Aetukuri
,
D.
Zhu
,
H. Y.
Hwang
,
N. C.
Brandt
,
H.
Wen
,
A. J.
Sternbach
,
S.
Bonetti
,
A. H.
Reid
,
R.
Kukreja
,
C.
Graves
,
T.
Wang
,
P.
Granitzka
,
Z.
Chen
,
D. J.
Higley
,
T.
Chase
,
E.
Jal
,
E.
Abreu
,
M. K.
Liu
,
T. C.
Weng
,
D.
Sokaras
,
D.
Nordlund
,
M.
Chollet
,
R.
Alonso-Mori
,
H.
Lemke
,
J. M.
Glownia
,
M.
Trigo
,
Y.
Zhu
,
H.
Ohldag
,
J. W.
Freeland
,
M. G.
Samant
,
J.
Berakdar
,
R. D.
Averitt
,
K. A.
Nelson
,
S. S. P.
Parkin
, and
H. A.
Dürr
,
Phys. Rev. B
98
(
4
),
045104
(
2018
).
11.
J. B.
Pendry
,
A. J.
Holden
,
D. J.
Robbins
, and
W. J.
Stewart
,
IEEE Trans. Microwave Theory Tech.
47
(
11
),
2075
2084
(
1999
).
12.
J. Y.
Ou
,
E.
Plum
,
L.
Jiang
, and
N. I.
Zheludev
,
Nano Lett.
11
(
5
),
2142
2144
(
2011
).
13.
R. A.
Shelby
,
D. R.
Smith
, and
S.
Schultz
,
Science
292
(
5514
),
77
79
(
2001
).
14.
D.
Schurig
,
J. J.
Mock
,
B. J.
Justice
,
S. A.
Cummer
,
J. B.
Pendry
,
A. F.
Starr
, and
D. R.
Smith
,
Science
314
(
5801
),
977
980
(
2006
).
15.
N. I.
Landy
,
S.
Sajuyigbe
,
J. J.
Mock
,
D. R.
Smith
, and
W. J.
Padilla
,
Phys. Rev. Lett.
100
(
20
),
207402
(
2008
).
16.
H.
Tao
,
C. M.
Bingham
,
A. C.
Strikwerda
,
D.
Pilon
,
D.
Shrekenhamer
,
N. I.
Landy
,
K.
Fan
,
X.
Zhang
,
W. J.
Padilla
, and
R. D.
Averitt
,
Phys. Rev. B
78
(
24
),
241103
(
2008
).
17.
A.
Rose
,
D. A.
Powell
,
I. V.
Shadrivov
,
D. R.
Smith
, and
Y. S.
Kivshar
,
Phys. Rev. B
88
(
19
),
195148
(
2013
).
18.
A.
Rose
,
D.
Huang
, and
D. R.
Smith
,
Phys. Rev. Lett.
107
(
6
),
063902
(
2011
).
19.
I.
Gil
,
J.
Garcia-Garcia
,
J.
Bonache
,
F.
Martin
,
M.
Sorolla
, and
R.
Marques
,
Electron. Lett.
40
(
21
),
1347
1348
(
2004
).
20.
B.
Wang
,
J.
Zhou
,
T.
Koschny
, and
C. M.
Soukoulis
,
Opt. Express
16
(
20
),
16058
16063
(
2008
).
21.
J.
Zhang
,
X.
Zhao
,
K.
Fan
,
X.
Wang
,
G.-F.
Zhang
,
K.
Geng
,
X.
Zhang
, and
R. D.
Averitt
,
Appl. Phys. Lett.
107
(
23
),
231101
(
2015
).
22.
H. R.
Seren
,
J.
Zhang
,
G. R.
Keiser
,
S. J.
Maddox
,
X.
Zhao
,
K.
Fan
,
S. R.
Bank
,
X.
Zhang
, and
R. D.
Averitt
,
Light: Sci. Appl.
5
,
e16078
(
2016
).
23.
N.
Kakenov
,
M. S.
Ergoktas
,
O.
Balci
, and
C.
Kocabas
,
2D Mater.
5
(
3
),
035018
(
2018
).
24.
D.
Zhang
,
M.
Trepanier
,
O.
Mukhanov
, and
S. M.
Anlage
,
Phys. Rev. X
5
(
4
),
041045
(
2015
).
25.
C.
Zhang
,
B.
Jin
,
J.
Han
,
I.
Kawayama
,
H.
Murakami
,
X.
Jia
,
L.
Liang
,
L.
Kang
,
J.
Chen
,
P.
Wu
, and
M.
Tonouchi
,
New J. Phys.
15
(
5
),
055017
(
2013
).
26.
K. G.
Nathaniel
,
G. P.
Bradford
, Jr.
,
Y. H.
Harold
,
C. B.
Nathaniel
,
T.
Darius
,
S.
Ranjan
,
Y.
Li
,
T.
Daniel
,
A. T.
Stuart
,
Q. X.
Jia
,
J. T.
Antoinette
,
A. N.
Keith
, and
C.
Hou-Tong
,
New J. Phys.
15
(
10
),
105016
(
2013
).
27.
W. L.
Chan
,
H.-T.
Chen
,
A. J.
Taylor
,
I.
Brener
,
M. J.
Cich
, and
D. M.
Mittleman
,
Appl. Phys. Lett.
94
(
21
),
213511
(
2009
).
28.
O.
Balci
,
N.
Kakenov
,
E.
Karademir
,
S.
Balci
,
S.
Cakmakyapan
,
E. O.
Polat
,
H.
Caglayan
,
E.
Özbay
, and
C.
Kocabas
,
Sci. Adv.
4
(
1
),
eaao1749
(
2018
).
29.
Y. H.
Fu
,
A. Q.
Liu
,
W. M.
Zhu
,
X. M.
Zhang
,
D. P.
Tsai
,
J. B.
Zhang
,
T.
Mei
,
J. F.
Tao
,
H. C.
Guo
,
X. H.
Zhang
,
J. H.
Teng
,
N. I.
Zheludev
,
G. Q.
Lo
, and
D. L.
Kwong
,
Adv. Funct. Mater.
21
(
18
),
3589
3594
(
2011
).
30.
X.
Zhao
,
J.
Zhang
,
K.
Fan
,
G.
Duan
,
J.
Schalch
,
G. R.
Keiser
,
R. D.
Averitt
, and
X.
Zhang
,
Phys. Rev. B
99
(
24
),
245111
(
2019
).
31.
R.
Marques
,
F.
Mesa
,
J.
Martel
, and
F.
Medina
,
IEEE Trans. Antennas Propag.
51
(
10
),
2572
2581
(
2003
).
32.
D. A.
Powell
,
K.
Hannam
,
I. V.
Shadrivov
, and
Y. S.
Kivshar
,
Phys. Rev. B
83
(
23
),
235420
(
2011
).
33.
H. A.
Haus
and
W.
Huang
,
Proc. IEEE
79
(
10
),
1505
1518
(
1991
).
34.
K. S.
Kunz
and
R. J.
Luebbers
,
The Finite Difference Time Domain Method for Electromagnetics
(
CRC Press
,
1993
).
35.
E.
Ekmekci
,
A. C.
Strikwerda
,
K.
Fan
,
G.
Keiser
,
X.
Zhang
,
G.
Turhan-Sayan
, and
R. D.
Averitt
,
Phys. Rev. B
83
(
19
),
193103
(
2011
).
36.
D. A.
Powell
,
M.
Lapine
,
M. V.
Gorkunov
,
I. V.
Shadrivov
, and
Y. S.
Kivshar
,
Phys. Rev. B
82
(
15
),
155128
(
2010
).
37.
T. C.
Tan
,
E.
Plum
, and
R.
Singh
,
Photonics
6
(
3
),
75
(
2019
).