Nanoplasmonics in Metallic Nanostructures and Dirac Systems

Abstract

In this book chapter, we review some of the progress made in nanoplasmonics and related optoelectronics phenomena in the field of two-dimensional (2D) materials and the recent 3D Weyl semimetals. We give a brief overview of plasmonics for three-dimensional (3DEG) and two-dimensional electron gases and draw comparisons with graphene, 3D topological insulators, 3D Weyl semimetals, and nanoplasmonics in nanogeometries. We discuss the decay of plasmons into electron-hole pairs and the subsequent thermalization and cooling of the hot carriers. We present our recent results in the fields of plasmonics in different nanostructures made of noble metals, such as Silver, and plasmonics in Dirac systems such as graphene and 3D topological insulators. We show a possibility of dynamically shifting the plasmon resonances in hybrid metal-semiconductor nanostructures. Plasmonics in 3D topological insulator and 3D Weyl semimetals have been least explored in nanoplasmonics although it can provide a variety of interesting physical phenomena involving spin plasmonics and chirality. Due to the inherent large spin-orbit coupling, locked spin-momentum oscillations can exist under special conditions and in the presence of an external laser field. We explore symmetric and antisymmetric modes in a slab of 3D TIs and present their dependences on the thickness of the slab.

1. Introduction

Electromagnetic properties of metal-dielectric interfaces have attracted a vast amount of research efforts. Ever since, the work of Mie [1] for small particles and Ritchie [2] for flat interfaces, a wide variety of scientists ranging from physicist, chemists, material scientists to biologists have explored plasmonics-based phenomena and their potential applications in practical life. In nanostructures, under the right circumstances, light waves propagating at metal-dielectric interface excite collective modes of electrons at the metal surface, resulting in the generation of charge density waves called surface plasmons (SPs), which can be divided into localized modes called localized surface plasmons (LSPs) and propagating modes called surface plasmon polaritons (SPPs), propagating along the interface like ripples across the surface of water with an effective wavelength much less than that of incident electromagnetic wave. Free electrons respond collectively by oscillating in resonance with the light waves. In optics, scientists have investigated methods to use plasmonics for concentrating, channeling, and changing the phase of light using subwavelength metallic structures. This would lead to miniaturized plasmonic circuits with length scales much smaller than those in current use [3–13]. A creatively engineered metal-dielectric interface can generate surface plasmons with the same frequency as the outside electromagnetic waves but with several times shorter wavelength. This interesting phenomenon can be utilized in a way in which surface plasmons carry information in microprocessors faster than current electronic transistors [14]. Plasmonics holds promise for a higher information density than conventional electronics [15]. While this proposed application needs still to be proven feasible, metallic nanostructures much smaller than the wavelength of light have already been successfully built for amplifying signals in surface enhanced Raman spectroscopy (SERS), providing a powerful method to detect a single molecule [16–27]. Plasmonically enhanced electric fields are already being used for sensing biomolecules [28–36]. The efficient heating property of plasmonic fields can be used for photothermal cancer treatment [37, 38] and also for thermally assisted magnetic recording [39–41]. Plasmonic lasers are able to achieve ultrafast dynamics with sub-wavelength mode confinement [42–50].

Metallic nanoparticles can be made in different sizes and shapes, and the distance among them can be controlled as well. These parameters can be used to tune the plasmon resonance frequency [51–62]. The use of optically excited plasmons as a tunable frequency source that can be mixed with a laser through Raman scattering enables dynamical shifting of the wavelength of light in a controlled manner [63]. Fluegel et al. [64] used a continuous laser beam of a few microwatts power to excite carriers in a carefully engineered narrow GaAs quantum well. These photogenerated carriers are selectively passed through a thick barrier of AlAs into a wide GaAs quantum well in which the two-dimensional electron gas (2DEG) supports high charge density waves with collective electron motion normal to the layer. A signal laser beam operating at a different wavelength undergoes inelastic Raman scattering from the plasmon-phonon (longitudinal optical phonon) modes in the two-dimensional (2D) quantum well. The result is that a signal with 13–15 nm redshifted frequency is generated. Plasmonics can also be exploited in optical tweezers to confine nanoparticles to small dimensions [63]. Grigorenko et al. [65, 66] have made electromagnetically coupled gold pillars. A trapping beam simultaneously excites gap plasmon (GP) modes in the gap between the pillars when they are used as a substrate, resulting in an enhanced field at the trapping site.

In a metallic nanoshell containing a core of dielectric material, due to different dielectric environments in the core and outside of the nanoshell, plasmons are excited in the inner and outer surface of the metallic nanoshell, which can be tuned by varying the ratio of inner to the outer radii of the nanoshell and by altering the dielectric materials [67, 68]. Figure 1 shows that the plasmon resonance energy can be tuned over a wide range by varying the material and/or shaping parameters. Figure 1 compares the plasmon resonance wavelengths of gold, silver, graphene, and topological insulator nanostructures. For a particle diameter much smaller than the wavelength of light, the light-matter interaction leads to an oscillating homogeneous polarization of the particle volume, resulting in an oscillating dipole field. For spherical Au and Ag nanoparticles, the dipole plasmon resonance occurs in the visible part of the spectrum. In particle ensembles, additional shifts are noticed due to electromagnetic coupling between LSP modes. For example, in case of a dimer, plasmons can be viewed as bonding and antibonding combinations, in analogy to molecular orbitals, i.e., hybridization of the individual nanoparticle LSPs occurs, giving rise to resonance shifts following the 1/d³ interaction between two classical dipoles [69]. The observed shifts in LSP resonance and the linewidth narrowing can be understood by means of the Mie scattering theory [1]. Due to the coupling of LSPs, the arrays support one longitudinal and two transverse modes of propagating polarization waves. The transport of energy along such a chain is the main idea for transmitting information using plasmonics [70, 71]. Large heat dissipation has so far slowed down the progress in this field.

Figure 1.

Comparison of the plasmon resonance wavelength spectra of gold, silver, graphene, and topological insulator (TI) nanostructures. The lifetime of plasmon for gold and silver ( $\tau =1-100 fs$) is shorter than that of graphene and TI ( $\tau =10-100 fs$).

Nevertheless, we proposed a method of guiding electromagnetic waves along a chain below the diffraction limit in a controlled manner using a chain of nanoshells [72, 73]. The dynamic control over the plasmon resonances and their coupling gave rise to the idea of an optically controlled plasmonic switch. In this chapter, we present a concept of dynamical control over the plasmon resonances that can be obtained by controlling the dielectric environment of the LSPs using a pump probe technique. In Ref. [74], we show that a shift up to 125 nm can be achieved in an Ag core-TiO₂ coated nanostructure.

In addition to that, we also review plasmonics in Dirac systems. Like in the case of a two-dimensional electron gas (2DEG), the surface plasmon mode in 2D materials such as graphene is tunable by changing the gate voltage through shifting the Fermi energy. Exciting surface plasmons give rise to light absorption enhancement in graphene, which can be utilized for photodetectors based on surface plasmon polaritons [75], optical switching of infrared plasmon polariton [76], and THz plasmonic lasing [77]. Plasmonics in Dirac systems show interesting features due to massless electrons around the Dirac nodes. In particular, we focus on surface plasmons in graphene, in Bismuth-based 3D topological insulators (3DTIs), and in 3D Weyl semimetals. Graphene, 3DTIs, and 3D Weyl semimetals are interesting due to their special electronic and optical properties arising from the linear dispersion relation around the Dirac cones in the Brillouin zone [78]. Around these points, energy dispersion of electrons can be described by a low energy Dirac Hamiltonian: H_G(k) = h(k) ⋅ σ, where h(k) = [h_x(k), h_y(k), 0] and σ = [σ_x, σ_y, σ_z]. For a small wave vector, q = k - K, we have h(k) = ħυ_Fq, showing that the electronic states are helical in nature, therefore allowing for “one-way traffic.” This helical nature of electronic states significantly enhances the charge mobility, resulting in reduced plasmonic loss, and therefore high-quality factor in doped graphene.

Due to the excitation of surface plasmons, it is now possible to engineer the behavior of light on nanometer length scales and to increase the light-matter interaction [79–82]. This interaction is an outcome of the near-field enhancement close to the metal surface, which also leads to plasmon damping through radiative decay and through nonradiative decay inside the material, due to Landau damping, i.e., creation of electron-hole pairs via interband or intraband transitions, electron-phonon interaction, and boundary effects. The intraband transition happens in the conduction band and the interband transition occurs between other bands (such as the d-band) and the conduction band, as shown schematically in Figure 2a [83–85].

Figure 2.

(a) Different types for decay of plasmon. The plasmon is damped by radiative (left) or nonradiative decay through Landau damping (right). (b) Schematic representation of the surface plasmon excitation followed by damping through hot-carrier generation. Schematics for the surface-assisted excitation of hot electron-hole (e-h) pair with energy $ħ{\omega }_p$. (c) An external electric field incident on a metal nanoparticle with characteristic size of L excites surface plasmon. (d) It leads to geometry-assisted hot e-h excitation (e) followed by momentum relaxation through surface scattering in a region with size of the order of $v_F/\omega$.

The radiative decay part of plasmon damping is due to the direct photon emission by coherent electron oscillation. As the size of the nanoparticle increases, the radiative decay of the plasmon is more significant. For larger nanoparticle elements, the radiative decay component is the main reason of plasmon resonance broadening and weakening of the dipole strength. In contrast, decreasing the size of the nanoparticle lets the nonradiative component dominate the plasmon decay. For applications in information technology, a slow dephasing of optical polarization by electron oscillation is essential, which is characterized by the dephasing time (T₂). The dephasing time can be found experimentally by measuring the plasmon linewidth Г = 2ħ/T₂, where T₂ is related to the time constant of the inelastic decay of the plasmons T₁ and elastic scattering time constant T^* by the authors of Refs. [82–84].

$\frac{1}{T_{\text{2}}}=\frac{1}{2T_1}+\frac{1}{T^{\ast }}.$E1

T₁ is a combination of time constants related to the radiative and nonradiative decay processes T₁ = T_1,r + T_1,nr. Nonradiative decay of surface plasmon through Landau damping gives rise to electron-hole pair generation. When the energy of these excited carriers is much larger than the thermally excited electron-hole pairs at ambient temperature, they are called hot carriers. This phenomenon of hot-carrier generation finds lots of applications in energy harvesting, photo detection, and photocatalysis. The hot carriers can be used directly to drive chemical reactions [86] in metals and to drive photocurrents in a semiconductor for the photovoltaic effect or for photoelectrochemical processes [87, 88]. The energy distribution of the generated hot carriers is sensitive to the detailed band structure of the metal. For example, the energy of the excited holes in copper and gold is higher than for electrons; however, in aluminum and silver, the electron and hole energies are almost equal. The advantage of the asymmetry between the produced hot electron and holes by the surface plasmon helps to collect more hot carriers before inelastic scattering occurs. The energy distribution of the carriers generated by decay of plasmons, the subsequent elastic and inelastic scatterings and their transport are important for the design of the hot-carrier-based devices. The inelastic scatterings due to the Coulomb interaction lead to the thermalization of the hot carriers among each other, thereby reaching their carrier equilibrium temperature. The inelastic scatterings due to the electron-phonon interaction result in the cooling of the carriers, thereby reaching the equilibrium temperature of the lattice of the material. For example, in graphene, thermalization takes place on a time scale of 100 fs, while cooling is much slower and happens on a time scale of 10 ps.

This book chapter is organized as follows: in Section 2, we discuss the SP resonances in hybrid metal-semiconductor nanostructures. By altering the dielectric environment of nanostructures dynamically using pump-probe techniques, we show that it is possible to shift the SP resonance wavelength. In Section 3, we discuss in detail the plasmon excitations and their damping pathways in a three-dimensional electron gas (3DEG). We discuss both the radiative and nonradiative damping mechanisms of SPs in 3DEGs. In Section 4, we give a brief overview of plasmons in a two-dimensional electron gas (2DEG). In Section 5, we present the size dependent properties of the SPs in nanostructures. Graphene plasmonics and losses are discussed in Section 6. Section 7 is dedicated to the description of the SPs in Dirac systems. We focus on the SPs in 3DTI materials and 3D Weyl semimetals, and we discuss graphene plasmons as a limiting case of the 3DTI plasmons in the limit where the thickness of the 3DTI slab d → 0. Finally, we wrap up the book chapter with the conclusion in Section 8.

2. Surface plasmon resonances in metal nanostructures

A nanoparticle shows tunable optical properties under controlled variation of its geometry. In a pure Ag spherical nanoparticle in vacuum, for example, the plasmon resonance occurs at 320 nm. These plasmon modes are shifted if the nanoparticle is coated with dielectric materials. It has been shown that with increasing shell thickness, the local electric field enhancement factor peak increases and redshifts for ε₂ > ε₃, whereas the local field enhancement factor peak decreases and blueshifts nonlinearly for ε₂ < ε₃, where ε₂ and ε₃ are dielectric functions of the shell and the surrounding materials, respectively [89]. The electric field enhancement factor is defined by E_F = |E/E₀|². It can be easily calculated by means of the finite-difference time domain (FDTD) technique. Figure 3 shows the example of an Au nanoparticle with a radius of 50 nm in air with index of refraction n = 1.

Figure 3.

The resonant excitation of the Au nanoparticle causes large local electric fields close to the particle surface, which is obtained by means of finite-difference time domain (FDTD) calculations. Here, E_F(x,y) is shown in the xz plane at 494 nm.

The quasi-static approximation provides a good estimate for a nanoparticle size of around 1/10 or smaller of the incident light wavelength. For larger nanoparticles, due to the finite speed of light, retardation effects lead to a redshift of the plasmon resonance [90]. In Ref. [91], authors have found an analytical expression for a spheroid that takes into account the depolarization factors and that gives a good approximation for nanogeometries of size up to 150 nm. Figure 4 shows our results for the local field enhancement in the presence of an Ag nanocube. As expected from electrostatics, the largest enhancement occurs at the vertices of the Ag nanocube.

Figure 4.

Local field enhancement E_F for an Ag nanocube of 20 nm size. The surrounding medium has a refractive index of n = 1.65, which is the case for boron nitride (BN). The green arrow points in the propagation direction of the electromagnetic wave, and the red double arrow shows the polarization of the electromagnetic wave.

The optical resonances of a nanoshell exhibit enhanced sensitivity to its local dielectric environment relative to the solid nanoparticle, as shown in Figure 5. For a particle diameter less than the wavelength of light, the light-matter interaction leads to an oscillating homogeneous polarization of the particle volume, resulting in a dipole field. Figure 5a and b shows the dependence of plasmon resonances on the shell thickness and the size of a hybrid metal-semiconductor nanostructure of 15 nm diameter with an Ag core coated by TiO₂ shell obtained in the quasi-static approximation. The expressions are given in Ref. [67].

Figure 5.

Surface plasmons (SPs) field enhancement in a hybrid nanoshell with a Ag core and TiO₂ shell of size of 15 nm at different shell thickness using the quasi-static field approximation (a), and different positions inside the shell (b) [see Ref. [74]. The inset in Figure 5a shows a scheme of a core-shell nanoparticle.

Plasmon resonances in a nanoshell can be tuned dynamically by letting a pump laser pulse of energy equal to the band gap or above generate electron-hole pairs in a semiconducting material surrounding the nanoshell. A probe laser pulse at a plasmon resonance frequency is used to excite plasmons on the metal surface. The generation of free electron-hole pairs alters the dielectric function of the surrounding semiconducting material. Due to the reduced dielectric function caused by the excitation of the electron-hole pairs, the excitation of surface plasmons by a probe pulse requires a higher energy. The frequency of the probe pulse is smaller than that of the pump pulse ensuring that no excitons are excited in the semiconductor during the probing. The change in the dielectric function of the surrounding medium due to the pump pulse can be calculated using Fermi’s golden rule:

$\text{Im}(\epsilon (\omega ))=\frac{4{\pi }^2{e}^2}{m^2{\omega }^2}{\displaystyle \sum _{k,c,v}}{\left|e·P_{c,v}\right|}^{2}f(E_v(k))\left[1-f(E_c(k))\right]\delta (E_c(k)-E_v(k)-ħ\omega ),$E2

where P_c,v is a transition dipole matrix element between conduction and the valence band, e is the direction of the polarization and f(E) is the Fermi-Dirac distribution function. At zero temperature, the factor $f(E_v(k))\left[1-f(E_c(k))\right]=1$ before the optical excitation and becomes less than 1 after optical excitation.

Figure 6 shows shifts in the resonance peak of the surface plasmons occurring at around 620 nm before the generation of excitons in a nanoshell structure (as shown in inset of Figure 5a) with diameter of 15 nm. After the pump pulse, depending on the density of the excitations, the plasmon resonance peaks are excited by the probe pulse shift. The larger the density of excited free electron-hole pairs in the semiconductors, the larger is the blueshift of the plasmon resonance peak. For a density of excitation of 5 × 10²¹ cm⁻³, a resonance shift of up to 125 nm can be achieved.

Figure 6.

Local field enhancement factor E_F on the surface of a nanoshell of 15 nm diameter. The core is made of Ag with a diameter of 7 nm. At higher density of excitation, the resonance peak is gradually blueshifted.

3. Plasmon excitation and damping for a three-dimensional electron gas (3DEG)

In this section, we discuss the plasmon excitation and their damping pathways in 3D materials made of metal. The plasmonic damping pathways in 3D materials include radiative decay, Landau damping, and resistive loss, as depicted in Figure 7. During Landau damping, plasmon quasi-particles lose their energy by exciting hot electron-hole pairs via direct interband or phonon/geometry-assisted intraband transitions. In the case of geometry-assisted intraband transition, the translational symmetry is broken due to electric field confinement or boundaries of the material [92]. In the case of resistive loss, single carriers, electrons or holes that are the building blocks of the plasmon quasi-particle, are kicked out of the phase-coherent collective plasma oscillation through electron-electron or electron-phonon scattering, giving rise to plasmon damping. In 2D materials, the plasmons follow similar damping pathways. Figure 7 shows the stages of the plasmon decay, the initial nonequilibrium configuration after the excitation, the thermalization, and the cooling of the hot carriers. Some of these damping pathways can be used to inject hot carriers into other materials. For example, at a metal-semiconductor interface, hot electron-hole pairs can be separated by means of the Schottky barrier for the purpose of energy harvesting. In metal-graphene or metal-MoS₂ junctions, the surface plasmons can generate hot electron-hole pairs, thereby injecting electrons/holes into n/p-doped 2D materials, giving rise to hot carrier-induced doping [93] or even insulator-to-metal phase transitions [94].

Figure 7.

The stages of energy relaxation of excited plasmon. (a) The incident electric field on a graphene nanostructure excites the surface plasmon. (b) The plasmon quasi-particles lose their energies by different scattering processes to excite hot carriers, which have nonthermal distribution on the order of tens to hundreds of femtoseconds. (c) Electron-electron and other scattering processes redistribute the hot carriers and thermalized them in the order of one to hundreds of picoseconds to generate ling-lived hot carrier distribution. (d) Cooling stages and energy transferring to the substrate happen by phonon and geometry-assisted scattering processes over a longer time scale (hundreds of picoseconds to tens of nanoseconds) for relaxation of the hot carrier distribution to the equilibrium.

The initial distribution of the hot carriers can be estimated using the jellium model for metal nanoparticles and nanoshells [95, 96], but this approach cannot explain the material dependence of this process because the specific band structure of the metal is completely neglected. In order to capture the material properties, it is necessary to combine FDTD calculations for obtaining the plasmon modes with ab initio density functional theory (DFT) calculations to determine the initial energy distribution of generated hot carriers by decay of surface plasmon. Figure 7b shows the schematic of plasmonic hot-carrier generation [97].

The plasmon resonance frequency and dispersion can be obtained by evaluating the dynamic polarizability in the presence of the carrier-carrier Coulomb interaction. The dynamic polarizability in the random phase approximation (RPA) is given by

${\chi }^{\text{RPA}}(q,\omega )=\frac{{\chi }^0(q,\omega )}{{\epsilon }^{\text{RPA}}(q,\omega )},$E3

where χ⁰(q,ω) is the noninteracting (zeroth-order) polarizability (single pair bubble) and ${\epsilon }^{RPA}(q,\omega )={\epsilon }_m-v_c(q){\chi }^0(q,\omega )$, with ε_m being the permittivity of the environment, and $v_c(q)=e^2/2{\epsilon }_{0}q$ the Coulomb interaction between the carriers. The RPA method corresponds to the expansion of $1/{\epsilon }^{\text{RPA}}(q,\omega )$, leading to an infinite power series over the bubble diagrams.

If optical phonons are also considered, the effective dielectric function in the RPA expansion takes the form [98, 99]

${\epsilon }^{\text{RPA}}(q,\omega )={\epsilon }_m-v_c(q){\chi }^0(q,\omega )-{\epsilon }_m{\displaystyle {\displaystyle \sum _l}{v}_{\text{sph},l}(q,\omega )}{\chi }^0(q,\omega )-{\epsilon }_m{v}_{\text{oph}}(q,\omega ){\chi }_{j,j}^0(q,\omega ),$E4

where ε_m is the average of dielectric constants of the environment. The second term represents the effective Coulomb interaction of electrons in graphene, and $v_c(q)=e^2/2q{\epsilon }_0$ is the direct Coulomb interaction. The third term is the effective dielectric function for different phonon modes (l) coming from electron-electron interaction mediated by substrate optical phonons, which couple to the electrons by means of the Fröhlich interaction, i.e.,

$v_{\text{sph},l}(q,\omega )={\left|M_{\text{sph}}\right|}^2{G}_l^0(\omega ),$E5

where |M_sph|² is the scattering and $G_l^0$ is the free phonon Green’s function. The last term of Eq. (4) corresponds to the optical phonon mediated electron-electron interaction.

$v_{\text{oph}}(q,\omega )={\left|M_{\text{oph}}\right|}^2{G}^0(\omega ).$E6

Here |M_oph|² defines the scattering matrix element and G^o(ω) is the free phonon Green’s function. In Eq. (4), ${\chi }_{j,j}^0(q,\omega )$ is the current-current correlation function. This description is very general and can be applied to any metallic system.

For the simplest case, when only the carrier-carrier Coulomb interaction is present, we can derive the dynamical plasmon dispersion relation following standard textbooks [100]. The first step is to calculate noninteracting dynamical polarizability.

${\chi }^0(q,\omega )=\frac{1}{\Omega }{\displaystyle {\displaystyle \sum _{k,s}}\frac{n_{0k+q,s}-n_{0k,s}}{{\epsilon }_{k+q}-{\epsilon }_k-ħ\omega -iħ\eta }},$E7

where Ω is the volume of the sample, $n_{0k,s}$ is the equilibrium electron density, and ε_k is the single-particle energy dispersion of the electrons. Assuming |q| ≪ k_F, we can expand χ⁰ (q, ω) in q by taking advantage of the expansions ${\epsilon }_{k+q}={\epsilon }_k+q\cdot {\nabla }_k{\epsilon }_k+\cdots$ and $n_{0k+q}=n_{0k}+\frac{\partial n_0}{\partial \epsilon }q\cdot {\nabla }_k{\epsilon }_k+\cdots$. At zero temperature, we have $\frac{\partial n_0}{\partial {\epsilon }_k}=-\delta \left({\epsilon }_k-{\epsilon }_F\right)$ and the velocity can be approximated by the Fermi velocity, i.e., ${\nabla }_k{\epsilon }_k=ħv_k\approx ħv_F$, leading to the approximation

$\begin{array}{c}{\chi }^0(q,\omega )\approx -\frac{2}{{\left(2\pi \right)}^3}{\displaystyle \int d^{3}k\frac{q\cdot v_F\delta \left({\epsilon }_k-{\epsilon }_F\right)}{q\cdot v_F-\omega -i\eta }}\hfill \\ =\frac{2}{{\left(2\pi \right)}^2}{\displaystyle {\displaystyle \underset{-1}{\overset{1}{\int }}}d\left(\cos \theta \right){\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}}{k}^{2}dk}\frac{q{v}_F\cos \theta }{q{v}_F\cos \theta -\omega -i\eta }\delta \left({\epsilon }_k-{\epsilon }_F\right)}\hfill \\ =\frac{k_F^2}{2{\pi }^2ħv_F}\left[\frac{2}{3}{\left(\frac{q{v}_F}{\omega +i\eta }\right)}^2+\frac{2}{5}{\left(\frac{q{v}_F}{\omega +i\eta }\right)}^4+\cdots \right]\hfill \\ \approx \frac{k_F^3{q}^2}{3{\pi }^{2}m{\left(\omega +i\eta \right)}^2}\left[1+\frac{3}{5}\frac{q^2{v}_F^2}{{\left(\omega +i\eta \right)}^2}\right]=\frac{n_0{q}^2}{m{\left(\omega +i\eta \right)}^2}\left[1+\frac{3}{5}\frac{q^2{v}_F^2}{{\left(\omega +i\eta \right)}^2}\right].\hfill \end{array}$E8

Using the equation $k_F={\left(3{\pi }^2{n}_0\right)}^{1/3}\iff n_0=\frac{k_F^3}{3{\pi }^2}$ for 3D metals, one obtains

${\chi }^0(q,\omega )\approx \frac{{\rho }_0{q}^2}{m{\left(\omega +i\eta \right)}^2}\left[1+\frac{3}{5}\frac{q^2{v}_F^2}{{\left(\omega +i\eta \right)}^2}\right]=\frac{{\rho }_0{q}^2}{m{\left(\omega +i\eta \right)}^2}R{\left(q,\omega \right)}^2$E9

with

$R{\left(q,\omega \right)}^2=\left[1+\frac{3}{5}\frac{q^2{v}_F^2}{{\left(\omega +i\eta \right)}^2}\right].$E10

Consequently, the dynamical polarizability in the RPA and long wavelength regime is

$\begin{array}{c}{\chi }^{\text{RPA}}(q,\omega )=\frac{{\chi }^0(q,\omega )}{\epsilon (q,\omega )}=\frac{\frac{{\rho }_0{q}^2}{m{\left(\omega +i\eta \right)}^2}R{\left(q,\omega \right)}^2}{1-\frac{4\pi e^2{\rho }_0}{m{\left(\omega +i\eta \right)}^2}R{\left(q,\omega \right)}^2}=\frac{{\rho }_0{q}^{2}R{\left(q,\omega \right)}^2}{m{\left(\omega +i\eta \right)}^2-4\pi e^2{\rho }_{0}R{\left(q,\omega \right)}^2}\\ =\frac{{\rho }_0{q}^{2}R\left(q,\omega \right)}{2m{\omega }_p}\left\{\frac{1}{\omega +i\eta -{\omega }_{p}R\left(q,\omega \right)}-\frac{1}{\omega +i\eta +{\omega }_{p}R\left(q,\omega \right)}\right\},\end{array}$E11

where the quantum plasma frequency of a 3D metal in the free electron approximation is given by the Pines and Bohm equation [101, 102]

${\omega }_p=\sqrt{\frac{4\pi n_0{e}^2}{m}},$E12

with n₀ being the equilibrium electron density and m the electron mass. Using ${\displaystyle \underset{\eta \to 0_{+}}{\lim }}\frac{1}{z-i\eta }=\mathcal{P}\left(\frac{1}{z}\right)+i\pi \delta (z)$, it follows that

$\mathrm{Im}\left[{\chi }^{\text{RPA}}(q,\omega )\right]\approx \frac{\pi {\rho }_0{q}^{2}R\left(q,\omega \right)}{2m{\omega }_p}\left\{\delta \left(\omega -{\omega }_{p}R\left(q,{\omega }_p\right)\right)-\delta \left(\omega +{\omega }_{p}R\left(q,{\omega }_p\right)\right)\right\},$E13

which yields the bulk plasmon dispersion relation for 3D metals, i.e.,

${\omega }_q={\omega }_{p}R\left(q,{\omega }_p\right)={\omega }_p\left[1+\frac{3}{10}\frac{q^2{v}_F^2}{{\left({\omega }_p+i\eta \right)}^2}+\cdots \right].$E14

The same result can be obtained by solving ${\epsilon }^{\text{RPA}}(q,\omega )=0$. Note that the slope of the parabolic dependence at $q=k_F$ is

$\frac{3v_F^2{k}_F}{5{\omega }_p}=\frac{3v_F^2{k}_F}{5\frac{4\pi e^2{n}_0}{m}}=\frac{9\pi m{v}_F^2{k}_F}{20e^2{k}_F^3}=\frac{9\pi m{v}_F^2}{20e^2{\left(\frac{m{v}_F}{ħ}\right)}^2}=\frac{9\pi {ħ}^2}{20e^{2}m}=0.67\frac{\text{m}}{\text{s}}\ll \frac{c}{\sqrt{{\epsilon }_{\text{d}}}}.$E15

This difference in slopes between the 3D bulk plasmon and the photon-like surface plasmon polariton is clearly visible in Figure 8 (see below).

Figure 8.

Dispersion curves of the RPA bulk plasmon resonance (red), the semi-classical surface plasmon polaritons (blue), and the RPA 2DEG plasmons. The shaded area marks the electron-hole continuum. When the bulk plasmon resonance enters the shaded area, the plasmon decays into electron-hole pairs which is called Landau damping.

It is well known that in 3D metals Landau damping occurs when the plasmon resonance energy enters the electron-hole continuum, which is determined by the condition

$\begin{array}{c}ħ\omega ={\epsilon }_{k+q}-{\epsilon }_k=\frac{{ħ}^2{\left(k+q\right)}^2-{ħ}^2{k}^2}{2m}=\frac{{ħ}^2{q}^2+2{ħ}^{2}k\cdot q}{2m}=\frac{{ħ}^2{q}^2+2{ħ}^{2}kq\cos \theta }{2m}\\ =\frac{{ħ}^{2}q\left(q+2k\cos \theta \right)}{2m}\ge 0.\end{array}$E16

Since $-1\le \cos \theta \le 1$ and $0\le k\le k_F$, the electron-hole continuum is given by the gray shaded area in Figure 8. The plasmon can decay into an electron-hole pair when the plasmon dispersion curve enters the electron-hole continuum limit. This decay corresponds to intraband Landau damping.

The presence of a planar boundary for a 3D metal adds a new mode known as surface plasmon, which propagates at the metal-dielectric interface. Since the electron charge density of a metal leaks outside the interface into the dielectric in the order of $1/k_F\approx 1$ Å, a macroscopic description based on Maxwell’s equations is sufficient to understand the surface plasmon qualitatively. Taking the boundary conditions into account, the dispersion relation of the surface plasmon is determined by Ritchie and Eldridge [103]

$q_x=\frac{\omega }{c}\sqrt{\frac{{\epsilon }_{\text{m}}{\epsilon }_{\text{d}}}{{\epsilon }_{\text{m}}+{\epsilon }_{\text{d}}}}=\frac{\omega }{c}\sqrt{\frac{{\omega }^2-{\omega }_p^2}{\left(1+\frac{1}{{\epsilon }_{\text{d}}}\right){\omega }^2-\frac{{\omega }_p^2}{{\epsilon }_{\text{d}}}}}$E17

in the Drude model, where ${\epsilon }_{\text{m}}\left(\omega \right)=1-{\omega }_p^2/{\omega }^2$. This biquadratic equation can be solved analytically, yielding

${\omega }_{\text{sp}}\left(q_x\right)=\pm \frac{1}{\sqrt{2{\epsilon }_{\text{d}}}}\sqrt{c^2{q}_x^2\left(1+{\epsilon }_{\text{d}}\right)+{\epsilon }_{\text{d}}{\omega }_p^2\pm \sqrt{c^4{q}_x^4{\left(1+{\epsilon }_{\text{d}}\right)}^2+2c^2{q}_x^2\left(-1+{\epsilon }_{\text{d}}\right){\epsilon }_{\text{d}}{\omega }_p^2+{\epsilon }_{\text{d}}^2{\omega }_p^4}}.$E18

Since ${\omega }_{\text{sp}}\left(q_x\right)\ge 0$, there are only two physical solutions, which are drawn in Figure 8. While the upper branch corresponds to the photon-like plasmon polariton, the lower branch represents the plasmon-like surface plasmon polariton (SPP) with the following asymptotics:

${\omega }_{\text{sp}}\left(q_x\right)=\{\begin{array}{cc}\frac{c}{\sqrt{{\epsilon }_{\text{d}}}}{q}_x\left[1-\frac{1}{2}{\left(\frac{c}{{\omega }_p}\right)}^2{q}_x^2+\dots \right],& q_x\ll \frac{{\omega }_p}{c}\\ \frac{{\omega }_p}{\sqrt{1+{\epsilon }_{\text{d}}}}\left[1-\frac{1}{2}{\left(\frac{{\omega }_p}{c}\right)}^2{\left(\frac{{\epsilon }_{\text{d}}}{1+{\epsilon }_{\text{d}}}\right)}^2{q}_x^{-2}+\dots \right],& q_x\gg \frac{{\omega }_p}{c}\end{array}$E19

which in the nonretarded regime reduces to Ritchie’s equation [2],

${\omega }_{\text{sp}}=\frac{{\omega }_p}{\sqrt{1+{\epsilon }_d}}=\frac{{\omega }_p}{\sqrt{2}}=\sqrt{\frac{2\pi n_0{e}^2}{m}}={\left(\frac{4e^2\sqrt{2m}}{3\pi {ħ}^3}\right)}^{1/2}{\epsilon }_F{}^{3/4}$E20

for wave vectors q in the range of ${\omega }_s/c\ll q\ll q_F$, where q_F is the Fermi wave vector [104]. Since $k_F={\left(3{\pi }^2{n}_0\right)}^{1/3}$ in 3D metals, the bulk plasmon frequency exhibits a ${\epsilon }_F{}^{3/4}\propto n_0^{1/2}$ dependence.

To estimate the plasmon decay rate, the band structure of the materials should be calculated by means of DFT to find out the exact quasi-particle orbitals and energies. Different electronic structures can be used to calculate the band structure of noble metals. To estimate the decay rate, the electronic states and energies of the metals resulting from PBEsol + U DFT method are used in Eq. (2). To find the plasmon decay distribution, energy conservation is used to determine the regions of the first Brillouin zone, which are responsible for the electron-hole generation [97]. Since aluminum, which is a nearly free electron metal, has a band crossing close to the Fermi level at the W point, interband transitions from the valence band with energies ħω from zero to the plasmon energy are allowed for creating a symmetric amount of both hot electrons and holes. There is an extra possibility for interband transitions at the K point, which mostly contribute to the generation of hot holes with energies >2 eV, giving rise to a modest asymmetry between hot electron and hole distributions. For copper, the interband transitions happen at the X, L and K, points, and for gold they occur at the X and L points [97]. In these two noble metals, since all the transitions originate from the “d” band, the hot hole energies are larger than 2 eV and consequently they have larger kinetic energies than the electrons. In metal-semiconductor interface, this asymmetry between hot hole and electron energy distribution plays a dominant role in the hot-carrier collection efficiency. For gold and copper, the holes have more kinetic energy than the electrons. This means that in a metal-semiconductor Schottky junction such as gold-n type gallium arsenide junctions, the hot electrons do not have enough energy to pass the Schottky barrier and therefore need to tunnel or get more thermal energy to overcome Schottky barrier. However, in gold-p type gallium arsenide junctions, the hot holes have sufficient kinetic energy to overcome the Schottky barrier. This means for copper and gold, the metal-p type semiconductor junction is more efficient for the collection of hot holes. But for aluminum and silver, the efficiency is almost the same for hot holes and electrons. According to the Fowler theory [96], the angular distribution of the hot-carriers should be isotropic. However, experiments show that none of the noble metals exhibit an isotropic distribution [97].

In nanoconfined structures, because of lack of translational symmetry, the crystal momentum q is not a good quantum number and there is no momentum conservation anymore. Therefore the transition creating an electron-hole pair is allowed to be nonvertical. This is the case for geometry-assisted intraband transition, where no phonon is required. A vector potential for one quantum of each symmetric and antisymmetric modes can be used in Fermi’s golden rule to calculate the decay rate of an LSP to an electron-hole pair. By decreasing the thickness, the effect of confinement is to be more pronounced and more hot electrons can be generated by intraband transitions. This effect is more noticeable for the antisymmetric mode, because a smaller part of the light is located inside the gold film, which leads to reduction of interband transitions [97].

In metals, plasmons decay is not only in the ultraviolet and visible spectrum but also in the infrared and microwave regimes [105, 106]. Due to the conservation of momentum in infinite crystal lattices, the direct interband transition induced by plasmon decay is only possible for energies larger than the band gap energy. However, for energies below the visible spectrum, typically phonon-assisted and surface-assisted intraband electron-hole pair generations are able to bypass this selection rule [92].

The plasmon decay rate is related to the imaginary part of the dielectric tensor $\mathrm{Im}{\displaystyle \overline{\epsilon }}\left(\omega \right)$ [92], i.e.,

$\text{Γ}(r)=\frac{1}{2\pi ħ}{E}^{\ast }(r)\cdot \mathrm{Im}{\displaystyle \overline{\epsilon }}\left(\omega \right)\cdot E(r).$E21

Let us consider a 3D semi-infinite metal slab extending in the negative z-direction with a single surface in the xy-plane at position z = 0. Using the electric field profile of a single quantum of the surface plasmon polariton with wave vector k and frequency ω and integration over the space gives the total nonradiative (nr) decay rate.

${\text{Γ}}_{\text{nr}}(r)=\frac{1}{2L(\omega )\left|\gamma (z<0)\right|}{\lambda }^{\ast }\cdot \mathrm{Im}{\displaystyle \overline{\epsilon }}\left(\omega \right)\cdot \lambda$E22

The decay rate of plasmon as a function of frequency can be calculated by substituting the experimental data for the complex dielectric function measured by ellipsometry. Within the random phase approximation (RPA), the nonradiative decay rate induced by direct interband transition is [97]

${\text{Γ}}_{\text{nr(direct)}}=\frac{2{\pi }^2{e}^2}{L(\omega )\left|\gamma (z<0)\right|m_e{}^2{\omega }^2}{\displaystyle {\displaystyle \underset{BZ}{\int }}\frac{dq}{{(2\pi )}^3}{\displaystyle {\sum }_{n^{\prime }n}(f_{qn}-f_{q{n}^{\prime }})}\delta \left({\epsilon }_{q{n}^{\prime }}-{\epsilon }_{qn}-ħ\omega \right)}\times {\left|\lambda \cdot {⟨p⟩}_{n^{\prime }n}^q\right|}^2$E23

where ${⟨p⟩}_{n^{\prime }n}^q=-{\Psi }_{q{n}^{\prime }}^{\sigma \ast }(r)iħ\nabla {\Psi }_{qn}^{\sigma }(r)$ are momentum matrix elements between the quasi-particle orbitals ${\Psi }_{q{n}^{\prime }}^{\sigma \ast }(r)$ and ${\Psi }_{qn}^{\sigma }(r)$. f_qn are the Fermi occupation functions of the quasi-particles, $\lambda \equiv {\displaystyle \widehat{k}}-{\displaystyle \widehat{z}}k/\gamma (z<0)$ is the polarization vector, and γ(z) is the imaginary part of the wave vector in z-direction. Eq. (23) can also be obtained in first quantization and is based on the electric dipole approximation, which means the plasmon mode function slowly changes on the atomic scale. To include the effect of electron-electron and electron-phonon scattering on the decrease in the lifetime of carriers, the δ function is replaced by a Lorentzian with half width $\mathrm{Im}{\Sigma }_{qn}+\mathrm{Im}{\Sigma }_{q{n}^{\prime }}$, where $\mathrm{Im}{\Sigma }_{qn}$ is the total carrier lifetime due to the scatterings.

The contribution of the electron-phonon interaction to the plasmon decay through intraband transitions can be calculated in second-order perturbation theory [107, 108], giving

$\begin{array}{l}{\text{Γ}}_{\text{nr(phonon)}}=\frac{2{\pi }^2{e}^2}{L(\omega )\left|\gamma (z<0)\right|m_e{}^2{\omega }^2}\underset{BZ}{{\displaystyle \int }}\frac{d{q}^{\prime }dq}{{(2\pi )}^6}{\displaystyle \sum }_{n^{\prime }n\alpha \pm }(f_{qn}-f_{q^{\prime }{n}^{\prime }})\left(N_{q^{\prime }-q\text{,}\alpha }+\frac{1}{2}\mp \frac{1}{2}\right)\hfill \\ \times \delta \left({\epsilon }_{q^{\prime }{n}^{\prime }}-{\epsilon }_{qn}-ħ\omega \mp ħ{\omega }_{q^{\prime }-q,\alpha }\right){\left|\lambda \cdot {\displaystyle {\sum }_{n_1}}\left(\frac{g_{q^{\prime }{n}^{\prime },q{n}_1}^{q^{\prime }-q,\alpha }{⟨p⟩}_{n_{1}n}^q}{{\epsilon }_{q{n}_1}-{\epsilon }_{qn}-ħ\omega +i\eta }+\frac{{⟨p⟩}_{n^{\prime }{n}_1}^{q^{\prime }}{g}_{q^{\prime }{n}_1,qn}^{q^{\prime }-q,\alpha }}{{\epsilon }_{q^{\prime }{n}_1}-{\epsilon }_{qn}\mp ħ{\omega }_{q^{\prime }-q,\alpha }+i\eta }\right)\right|}^2,\end{array}$E24

where $ħ{\omega }_{k,\alpha }$ is the energy of a phonon with wave vector k and polarization index α, N_k,α is the Bose occupation factor for the phonons, $g_{q^{\prime }{n}^{\prime },qn}^{k,\alpha }$ is the electron-phonon matrix element with q and n being the wave vector and band index of the electronic states, and the ± signs distinguish between absorption and emission of phonons.

Excitation of surface plasmon generates a strong field confinement on the surface with the exponential decay length $({|\gamma (z<0)|}^{-1})$ inside the metal, which creates a Lorentzian distribution for the momentum of the plasmon in the z-direction. This momentum distribution leads to diagonal intraband transition and introduces a “surface-assisted” plasmon decay [92]

$\mathrm{Im}{\epsilon }_{\text{surface}}\left(\omega \right)=\frac{{\omega }_p^2}{{\omega }^3}\times \left(\frac{3}{4}\left|\gamma (z<0)\right|v_F\right)\frac{2k^2}{k^2+{\left|\gamma (z<0)\right|}^2},$E25

where ${\omega }_p=\sqrt{4\pi n{e}^2/m_e}$ is the bulk plasma frequency of the metal, $v_F=\left(ħ/m_e\right){\left(3{\pi }^{2}n\right)}^{3/2}$ is the Fermi velocity, and n is the bulk carrier density of the metal.

Numerical studies based on the free electron jellium model show that in nanostructures, due to the localization of electronic states and the nonconservation of the crystal momentum, intraband transitions are enhanced [95, 96, 109]. Using Fermi’s golden rule together with the free electron eigen states and the dipole field profile, the nanoconfinement contribution is

$\mathrm{Im}{\epsilon }_{\text{sphere}}\left(\omega \right)=\frac{{\omega }_p^2}{{\omega }^3}\times \frac{6v_F}{{\pi }^{2}R},$E26

where R is the radius of the spherical nanoparticle.

In the Landau damping theory, the lowest-order processes consisting of direct, phonon-assisted, and surface/geometry-assisted electron-hole pair excitation contribute to the decay of plasmons [102, 110]. Higher-order processes leading to the excitation of many electron-hole pairs or many phonons are suppressed due to the phase-space factors at small energies [111]. Only at large energies, the higher-order processes become significant. Being completely different from Landau damping, another source of plasmon damping is the resistive loss in the metal, which can be calculated by means of the linearized Boltzmann equation in the relaxation time approximation [92], giving

$\mathrm{Im}{\epsilon }_{\text{resist}}\left(\omega \right)=\frac{4\pi {\sigma }_0}{\omega \left(1+{\omega }^2{\tau }^2\right)},$E27

where σ₀ is the zero frequency conductivity and τ is the momentum relaxation time, which can be derived from the carrier-carrier Coulomb and electron-phonon matrix elements.

4. Plasmons in a two-dimensional electron gas (2DEG)

According to Ritchie [2] and Stern [112], the dispersion relation for 2DEG plasmons is

${\omega }_{\text{2DEG}}=\sqrt{\frac{n_s{e}^2}{2m^{*}{\epsilon }_0{\displaystyle \overline{\epsilon }}\left(q,\omega \right)}{q}_x},$E28

where n_s is the surface carrier density. A typical dispersion curve is shown in Figure 8. In the case of a typical metal-oxide-semiconductor (MOS) system, the effective dielectric function is given by

${\displaystyle \overline{\epsilon }}\left(q,\omega \right)=\frac{1}{2}\left[{\epsilon }_{\text{SC}}\left(\omega \right)+{\epsilon }_{\text{OX}}\left(\omega \right)\coth \left(d{q}_x\right)\right],$E29

where ${\epsilon }_{\text{SC}}\left(\omega \right)$ and ${\epsilon }_{\text{OX}}\left(\omega \right)$ are the dielectric constants of the semiconductor and oxide layers, respectively. The oxide layer has a thickness of d. For example, Heitmann [113] and Wilkinson [114] were able to observe this type of 2DEG plasmons in AlGaAs-GaAs heterostructures. This result shows clearly that the dispersion curve can be tuned by changing the thickness of the oxide layer and by altering the gate voltage. Similarly, the plasmon dispersion can be tuned in atomically thin 2D materials. Since the dispersion relation in Eq. (28) exhibits a $\sqrt{q_x}$ dependence for a 2DEG, it is plausible to assume that atomically thin metallic 2D materials follow a similar dependence, which is indeed the case, as we discuss below in Sections 6 and 7. Since in a 2DEG

$n_s=\frac{1}{\Omega }{\displaystyle {\displaystyle \sum _{\left|k\right|\le k_F,s}}1}=2{\displaystyle \int \frac{d^{2}k}{{\left(2\pi \right)}^2}1}=2\frac{\pi k_F^2}{{\left(2\pi \right)}^2}=\frac{k_F^2}{2\pi }=\frac{m^{*}{\epsilon }_F}{\pi {ħ}^2},$E30

the 2DEG plasmon frequency exhibits a $\sqrt{n_s}\propto \sqrt{{\epsilon }_F}$ dependence. Note that the electron-hole continuum in a 2DEG is determined by a similar formula to Eq. (16). Therefore the gray shaded area in Figure 8 has similar forms for 3DEGs and 2DEGs.

5. Static geometry of metallic objects and environment

When the size of the nanoparticle is much smaller than the wavelength of the incident light, the particle exhibits a dipolar oscillation mode (Fröhlich mode). As the diameter of the nanoparticle is increased, the electrostatic limit is not a good approximation anymore and the multipolar oscillation modes start to appear. Excitation of these modes gives rise to the broadening of the resonance [115], as seen from Figure 9b. The LSP resonance redshifts with increasing diameter of the sphere, which is due to retardation effects [90]. Decreasing the size of nanoparticle less than mean free path of electrons moves the material band structure and dielectric function away from the bulk properties and increases the surface scattering that gives rise to broadening of the absorption spectrum [115], as shown in Figure 9a. The internal field enhancement of an illuminated spherical nanoparticle is

Figure 9.

Calculated extinction coefficient per unit volume for a spherical Aluminum nanoparticle of various diameters. (a) For the case when the size of the nanoparticle becomes smaller than the mean free path of the electron, the linewidth becomes wider due to increased surface scattering. (b) For the case when the size of the nanoparticle becomes much larger than the mean free path of the electron, excitation of multipole modes leads to broadening of the linewidth and the retardation effect leads to a redshift of the LSP resonance peak.

$\frac{E_i\left(\omega \right)}{E_0}=\frac{3{\epsilon }_m}{\epsilon \left(\omega \right)+{\epsilon }_m},$E31

where ε and ε_m are the dielectric functions of the nanoparticle and the medium, respectively. An ellipsoidal nanoparticle with a shape defined by ${\left(x/a\right)}^2+{\left(y/b\right)}^2+{\left(z/c\right)}^2=1$ supports an internal field enhancement given by

$\frac{E_i(\omega )}{E_0}=\frac{{\epsilon }_m}{{\epsilon }_m+L(\epsilon (\omega )-{\epsilon }_m)},$E32

where L is the shape factor, whose value lies between 0 and 1. For the resonance frequency, the real part of the denominator is zero, which means the resonance occurs at frequencies when

$\epsilon \left(\omega \right)={\epsilon }_m(1-\frac{1}{L})$E33

is satisfied. The shape factor for the axis parallel to the polarization of the incident light (e.g., the axis a) is given by

$L=\frac{abc}{2}{\displaystyle {\int }_0^{\infty }\frac{dq}{(a^2+q)f(q)}},$E34

where

$f(q)={\left[(a^2+q)(b^2+q)(c^2+q)\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}.$E35

For the core-shell spherical nanoparticles the resonance condition is given by

${\epsilon }_c=-2{\epsilon }_s\left[\frac{{\epsilon }_s(1-f)+{\epsilon }_h(2+f)}{{\epsilon }_s(1+2f)+2{\epsilon }_h(1-f)}\right],$E36

where $f={\left(r_c/r_s\right)}^3$ is the core fill fraction and ε_c, ε_s, and ε_h are the core, shell, and host dielectric functions, respectively. The resonance frequency is tunable by modifying the shell thickness and the radius of the core. As shown in Figure 10, due to the hybridization of the plasmon modes, even in the electrostatic limit there are multiple resonance frequencies [116].

Figure 10.

Energy-level diagram of plasmon hybridization in metallic core-shell nanoparticles. The hybridization results from the coupling between the cavity and sphere surface plasmons.

6. Plasmon theory for graphene

Graphene is a two-dimensional (2D) material comprised of a single layer of carbon atoms in a honeycomb lattice. It has unique electrical, optical, and mechanical properties due to its tunable band dispersion relation and atomic thickness. Because of its unique band structure graphene possesses a very high mobility and a fast carrier relaxation time [117–121], making it an attractive candidate for ultrafast electronics and optoelectronics. Exciting surface plasmons on graphene is a distinct technique to increase absorption with low damping rate. The surface plasmon couples the electromagnetic (EM) wave to the conductive medium, giving rise to direct absorption of light by monolayer graphene and providing the opportunity of electrical tunability of the plasmon resonance frequency, high degree of electric field confinement, and low plasmon damping rate [122–125]. The increased light-matter interaction results in an enhanced spontaneous emission rate close to the nanostructure edges [126, 127]. Recent experiments have achieved an absorption of 90% in the mid-IR range by connecting graphene with high carrier mobility to a silicon diffractive grating [128] and designing graphene nanoribbons [128, 129], nanodisks [130], and antidot array [125] theoretically. These high carrier mobilities can be achieved only for mechanically exfoliated graphene. Exciting plasmon on CVD-grown monolayer graphene with lower mobility than the mechanically exfoliated one reduces the absorbance to 19% and 28% for graphene nanoribbons [99, 127, 131] and nanodisks [132, 133], respectively. We show in experiments that the coupling of a patterned CVD-grown graphene sheet to an optical cavity amplifies the excited LSPs and enhances the light absorption to a current world record of 45% [134]. We also show that the theoretically achievable enhancement is 60% for a square lattice of holes [134].

The electric current of graphene in the interaction picture is given by

$j_{xx}(\omega )=\sigma (\omega )E(\omega ),$E37

where E(ω) is the in-plane electric field and σ(ω) is the interband optical conductivity

${\sigma }_{\text{inter}}(\omega )=-\frac{i{e}^2}{\omega A}{\displaystyle {\displaystyle \sum _k}{v}_x{}^2\frac{{\rho }_{vv}-{\rho }_{cc}}{{\epsilon }_{c,k}-{\epsilon }_{v,k}-ħ(\omega +i\eta )}}$E38

with A being the cross-section area, v_x the velocity operator along x-direction, ρ_vv and ρ_cc the digonal elements of the density matrix, ε_c,k /ε_v,k the conduction/valence band energy, and η an infinitesimal value. Using the energy dispersion, ${\epsilon }_{n,k}=n{v}_Fħ\left|k\right|$ with $n=+1/-1$ for conduction/valence band and $v_F={10}^{6}m/s$, and replacing the summation by the integration in Eq. (38) leads to

$\begin{array}{ll}{\sigma }_{\text{inter}}(\omega )\hfill & =\frac{e^2{v}_f{}^2}{2\omega }\frac{g}{4\pi }\underset{0}{\overset{\infty }{{\displaystyle \int }}}2\pi kdk({\rho }_{vv}-{\rho }_{cc})\delta ({\epsilon }_{c,k}-{\epsilon }_{v,k}-ħ\omega )\hfill \\ \hfill & \begin{array}{c}=\frac{e^2{v}_f{}^2}{\omega }\underset{0}{\overset{\infty }{{\displaystyle \int }}}kdk(f({\epsilon }_v)-f({\epsilon }_c))\delta ({\epsilon }_{c,k}-{\epsilon }_{v,k}-ħ\omega )\end{array}\hfill \\ \hfill & \begin{array}{c}=\frac{e^2}{4ħ}\left[f\left(-\frac{ħ\omega }{2}\right)-f\left(\frac{ħ\omega }{2}\right)\right],\end{array}\hfill \end{array}$E39

Where g = 4 is the graphene degeneracy and ρ is replaced by the Fermi distribution function f. By inserting the explicit expression of the Fermi distribution, Eq. (39) is reduced to

${\sigma }_{\text{inter}}\left(\omega \right)=\frac{e^2}{4ħ}\frac{\sinh \left(\frac{ħ\omega }{2k_{B}T}\right)}{\cosh \left(\frac{{\epsilon }_F}{2k_{B}T}\right)+\cosh \left(\frac{ħ\omega }{2k_{B}T}\right)},$E40

where ε_F is the Fermi energy of graphene, k_B is the Boltzmann constant and T is the temperature. Eq. (40) is the graphene optical conductivity due to interband transition, which is valid in the visible light spectrum. This result can be derived from the Kubo formula too. When the energy of the incident light (ħω) and k_BT are of the same order and ε_F is much larger than k_BT, the interband optical conductivity is constant, i.e.,

${\sigma }_{\text{inter}}\left(\omega \right)=\frac{e^2}{4ħ}.$E41

For infrared and THz radiation, the Fermi energy can be tuned to become much larger than the incident photon energy, and therefore due to Pauli blocking there are only intraband transitions. According ot the Boltzmann equation and under the relaxation time approximation, the carrier distribution in the presence of a constant electric field with x-polarization is given by

$f(k)=f^0(k)-q{E}_x{v}_x(k)\tau (k)\frac{\partial f^0(k)}{\partial \epsilon },$E42

where f ⁰(k) is the carrier distribution in the absence of the electric field, q is the carrier charge, and τ(k) is the relaxation time in the presence of electron-phonon interaction, electron-impurity interaction, and other scattering processes. By using Eq. (42), the electric current is given by

$j_x=\frac{q}{A}{\displaystyle {\displaystyle \sum _k}f(k)}{\nu }_x(k)=\frac{-q^2{E}_x}{A}{\displaystyle {\displaystyle \sum _k}{\nu }_x{}^2(k)\tau (k)\frac{\partial f^0(k)}{\partial \epsilon }}.$E43

In the absence of the external electric field, the net electric current is zero, so the summation of the first part of Eq. (43) in the electric current is zero. Since $\frac{\partial f^0(k)}{\partial \epsilon }$ is so narrow around the Fermi surface, only the wave vectors near the Fermi energy contribute to the integration, and τ(k) is approximately constant. By replacing ${\nu }_x{}^2(k)$ with ${\nu }_F{}^2/2$, the intraband optical conductivity becomes

${\sigma }_{\text{intra}}=-\frac{q^2\tau }{A}\frac{{\nu }_F{}^2}{2}{\displaystyle {\displaystyle \sum _k}\frac{\partial f^0(k)}{\partial \epsilon }}=-\frac{q^2\tau v_F{}^2}{2{\pi }^2}{\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}}dk2\pi k\frac{\partial f^0(k)}{\partial \epsilon }}.$E44

For each k there are two energy values in the conduction and the valence bands. Using the linear dispersion relation, the intraband optical conductivity is reduced to

$\begin{array}{c}{\sigma }_{\text{intra}}=-\frac{q^2\tau }{\pi {ħ}^2}\left({\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}}\epsilon d\epsilon \frac{\partial f^0(\epsilon )}{\partial \epsilon }+{\displaystyle {\displaystyle \underset{0}{\overset{-\infty }{\int }}}\epsilon d\epsilon \frac{\partial f^0(\epsilon )}{\partial \epsilon }}}\right)\hfill \\ =-\frac{q^2\tau }{\pi {ħ}^2}{\displaystyle {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}\left|\epsilon \right|d\epsilon \frac{\partial f^0(\epsilon )}{\partial \epsilon }}\hfill \\ =\frac{e^2\tau }{\pi {ħ}^2}2k_{B}T\ln \left(2\cosh \left(\frac{{\epsilon }_F}{2k_{B}T}\right)\right),\hfill \end{array}$E45

where q is replaced by the electron charge (−e). Assuming ${\epsilon }_F\gg k_{B}T$, which is usually the case, the intraband optical conductivty simplifies to

${\sigma }_{\text{intra}}=\frac{e^2}{\pi {ħ}^2}{\epsilon }_F.$E46

In the presense of an oscillating electric field, the relaxation time is a complex function $\left({\tau }^{-1}\to {\tau }^{-1}-i\omega \right)$ and the intraband conductivity is given by

${\sigma }_{\text{intra}}(\omega )=\frac{e^2}{\pi {ħ}^2}\frac{2k_{B}T}{{\tau }^{-1}-i\omega }\ln \left(2\cosh \left(\frac{{\epsilon }_F}{2k_{B}T}\right)\right),$E47

which in the case of ${\epsilon }_F\gg k_{B}T$ is reduced to [118, 120]

${\sigma }_{\text{intra}}(\omega )=\frac{e^2}{\pi {ħ}^2}\frac{{\epsilon }_F}{{\tau }^{-1}-i\omega },$E48

where τ is determined by impurity scattering and electron-phonon interaction ${\tau }^{-1}={\tau }_{\text{imp}}^{-1}+{\tau }_{\text{e-ph}}^{-1}$.

The dielectric function of graphene can be obtained via its AC conductivity by means of [135]

$\epsilon (\omega )=2.5+i\frac{\sigma (\omega )}{{\epsilon }_0\omega d}$E49

where ε_g = 2.5 is the dielectric constant of graphite. Here d is the thickness of graphene.

The bandstructure of graphene is linear in the tight-binding approximation, as shown in Figure 11b. The dispersion relation for the TM mode in the geometry depicted in Figure 11a, which consists of graphene surrounded by dielectrics with constants ε_r1 and ε_r2, is given by [122]

$\frac{{\epsilon }_{r1}}{\sqrt{q^2-\frac{{\epsilon }_{r1}{\omega }^2}{c^2}}}+\frac{{\epsilon }_{r2}}{\sqrt{q^2-\frac{{\epsilon }_{r2}{\omega }^2}{c^2}}}=-\frac{\sigma (q,\omega )i}{\omega {\epsilon }_0},$E50

Figure 11.

(a) Schematic of the patterned graphene layer surrounded by two media of dielectrics ε_r1 and ε_r1. (b) The sketch of intraband (green arrow) and interband (red arrow) transitions after absorbing a photon on the Dirac cone. (c) The graphene plasmon loss function for the structure shown in part (a) without plasmon-optical phonon coupling. This function shows the graphene plasmon dispersion relation. (d) The graphene dispersion relation in presence of substrate/graphene optical phonons.

where σ(ω,q) is the conductance of graphene. The semiclassical conductivity σ(ω) does not depend on q. In the nonretarded regime (q = ω/c), the dispersion relation of the TM mode is reduced to [122]

$q\approx {\epsilon }_0\frac{{\epsilon }_{r1}+{\epsilon }_{r2}}{2}\frac{2i\omega }{\sigma (q,\omega )}.$E51

The larger the q (smaller plasmon wavelength), the higher is the transversal confinement of the TM mode due to excitation of the surface plasmon (SP) on graphene, which is similar to the case of regular metals. The main difference of SPs on metal and 2D materials is due to their band structures. While normal metals have parabolic bands, graphene has Dirac cones, which leads to a novel dispersion relation for the SPs on graphene, as shown in Figure 11c. The interactions between the substrate and graphene optical phonons and the graphene plasmons open a band gap in the graphene plasmon dispersion relation, as shown in Figure 11d, showing the plasmon loss function in the presence of plasmon-optical phonons interaction given by Eqs. (4)–(6) [134]. The electrons in 2D materials are bound in the normal direction, which leads to more confinement of the coupled electromagnetic wave due to plasmon excitation. In graphene, thermalization takes place on a time scale of 100 fs, while cooling is much slower and happens on a time scale of 10 ps, as seen from Figure 11 [136, 137].

There are two different approaches to obtain the dispersion relation of the surface plasmon. In the semiclassical approximation, the Drude-like conductivity is used to obtain the plasmon dispersion relation [122, 138]

$q(\omega )=\frac{\pi {ħ}^2{\epsilon }_0\left({\epsilon }_{r1}+{\epsilon }_{r2}\right)}{e^2{\epsilon }_F}\left(1+\frac{i}{\tau \omega }\right){\omega }^2$E52

and the plasmon loss

$\frac{\mathrm{Re}q}{\mathrm{Im}q}=\omega \tau =\frac{2\pi c\tau }{{\lambda }_{\text{air}}}.$E53

Another approach that works for both intraband and interband regimes is based on the selfconsistent linear response theory, known as random phase approximation (RPA) along with the relaxation time approximation defined by Mermin [139]. The dispersion relation of the plasmon can be obtained by solving

${\epsilon }_{\text{RPA}}(q,\omega )=0,$E54

with the complex wave vector $q=q_1+i{q}_2$. Considering the Coulomb interaction of electrons in graphene and the medium dielectric, first and second terms of Eq. (4), one obtains the plasmon dispersion relation

$\frac{{\epsilon }_{r1}+{\epsilon }_{r1}}{2}+\frac{e^2}{2{\epsilon }_0{q}_1}\mathrm{Re}\chi (q_1,\omega )=0$E55

and the plasmon loss relation [122]

$q_2=\frac{\mathrm{Im}\chi (q_1,\omega )+\frac{1}{\tau }\frac{\partial }{\partial \omega }\mathrm{Re}\chi (q_1,\omega )+\frac{1}{\omega \tau }\mathrm{Re}\chi (q_1,\omega )\left[1-\frac{\chi (q_1,\omega )}{\chi (q_1,0)}\right]}{\frac{1}{q_1}\mathrm{Re}\chi (q_1,\omega )-\frac{\partial }{\partial q_1}\mathrm{Re}\chi (q_1,\omega )}.$E56

The plasmon losses $\left(\text{Req}/\text{Imq}\right)$, wave localization $\left({\lambda }_{\text{air}}/{\lambda }_p\right)$, and the group velocity of the graphene plasmon $\left({\nu }_g=\partial \omega /\partial q\right)$ are calculated by means of the semiclassical approximation and RPA, as shown in Figure 12 for the Fermi level of ${\epsilon }_F=0.135\text{eV}$ and relaxation time of $\tau =1.35\times {10}^{-13}s$ [122]. Below the interband regime, the plasmon loss, the wave localization, and the group velocity calculated by RPA are in very good agreement with the semiclassical approach.

Figure 12.

(a) Plasmon loss, (b) filed confinement, and (c) the group velocity of the graphene surface plasmon for ${\epsilon }_F=0.135\text{eV}$ and the carrier mobility of $\mu =10000{\text{cm}}^{\text{2}}/\text{Vs}$. The solid and dashed line represent the RPA and semiclassical results, respectively, and the rose shaded area is the high-interband loss region.

For 2D materials such as graphene, the amount of plasmon loss can be calculated by using the effective dielectric function of the material. The dynamical polarizability

$\chi (q,i{\omega }_n)=-\frac{1}{A}{\displaystyle {\int }_0^{\beta }d\tau e^{i{\omega }_n\tau }⟨T{\rho }_q(\tau ){\rho }_{-q}(0)⟩}$E57

determines several important quantities, such as the effective electron-electron interaction, plasmon spectra, phonon spectra, and Friedel oscillations. ${\omega }_n=\frac{2\pi n}{\beta }$ are the Matsubara frequencies, ρ_q is the density operator in q-space and A denotes the area. This quantity is calculated in the canonical ensemble for both of the sub-lattice density operators (ρ = ρ_a + ρ_b). The dynamical polarizability in zeroth-order in the electron-electron interaction in the long wavelength limit is [140]

${\chi }^0(q,i{\omega }_n)=\frac{g_s{g}_v}{4{\pi }^2}{\displaystyle \int d^{2}k{\displaystyle {\displaystyle \sum _{s,s^{\prime }=\pm }}{f}^{s{s}^{\prime }}(k,q)}}\frac{f({\epsilon }^s(k))-f({\epsilon }^{s^{\prime }}(\left|k+q\right|))}{{\epsilon }^s(k)-{\epsilon }^{s^{\prime }}(\left|k+q\right|)+iħ{\omega }_n},$E58

where $g_s=g_v=2$ are the spin and valley degeneracy, $f({\epsilon }^s(k))$ is the Fermi distribution function, and ${\epsilon }^s(k)=sħv_{F}k-{\epsilon }_F$ is the graphene energy. The band overlap of the wavefunctions

$f^{s{s}^{\prime }}(k,q)=\frac{1}{2}\left(1+s{s}^{\prime }\frac{k+q\cos \phi }{\left|k+q\right|}\right)$E59

is a specific property of graphene, where φ signifies the angle between k and q. At zero temperature, the Fermi-Dirac distribution functions are step functions. In that case, using ${\omega }_n\to \omega +i\delta$, Eq. (58) is simplified to

${\chi }_{T=0}^{0\pm }(q,\omega )=\frac{g_s{g}_v}{4{\pi }^2ħ}{\displaystyle \int d^{2}k{\displaystyle {\displaystyle \sum _{\alpha =\pm }}\frac{\alpha f^{\pm }(k,q)}{\omega +\alpha v_F(k\mp \left|k+q\right|)+i\delta }}},$E60

where the + and – signs denote the intra and interband transitions, respectively. Integration over φ and k gives the retarded polarizability or charge-charge correlation function

${\chi }^0(q,\omega )={\chi }_{{\epsilon }_F=0}^0(q,\omega )+\text{Δ}{\chi }^0(q,\omega ),$E61

where

${\chi }_{{\epsilon }_F=0}^0(q,\omega )=-i\pi \frac{\mathfrak{I}(q,\omega )}{{ħ}^2{\nu }_F^2},$E62

and

$\text{Δ}{\chi }^0(q,\omega )=-\frac{g{\epsilon }_F}{2\pi {ħ}^2{v}_F^2}+\frac{\mathfrak{I}(q,\omega )}{{ħ}^2{v}_F^2}\left\{\begin{array}{l}\aleph \left(\frac{ħ\omega +2{\epsilon }_F}{ħv_{F}q}\right)-\text{Θ}\left(\frac{2{\epsilon }_F-ħ\omega }{ħv_{F}q}-1\right)\times \left(\aleph \left(\frac{2{\epsilon }_F-ħ\omega }{ħv_{F}q}\right)-i\pi \right)\\ -\text{Θ}\left(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q}+1\right)\times \aleph \left(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q}\right)\end{array}\right\}$E63

The two functions $\mathfrak{I}(q,\omega )$ and $\aleph \left(x\right)$ are defined by

$\mathfrak{I}(q,\omega )=\frac{g}{16\pi }\frac{ħv_F^2{q}^2}{\sqrt{{\omega }^2-{\nu }_F^2{q}^2}}$E64

and

$\aleph \left(x\right)=x\sqrt{x^2-1}-\ln (x+\sqrt{x^2-1}),$E65

Where $g=g_s{g}_v=4$. For $\omega >q{v}_F$ and in the long wavelength limit q → 0, $x=\left|(ħ\omega \pm 2{\epsilon }_F)/ħv_{F}q\right|\gg 1$, thus $x^2-1\approx x^2$ and $\aleph \left(x\right)\approx x^2-2\ln (x)$. We derive here the dynamical polarizability in Eq. (63) in these regimes. The expansion of $\mathfrak{I}(q,\omega )$ gives

$\mathfrak{I}(q,\omega )=\frac{g}{16\pi }\frac{ħv_F^2{q}^2}{\omega }{\left(1-\frac{v_F^2{q}^2}{{\omega }^2}\right)}^{-1/2}\approx \frac{g}{16\pi }\frac{ħv_F^2{q}^2}{\omega }\left(1+\frac{v_F^2{q}^2}{2{\omega }^2}\right).$E66

In this condition and for intraband transitions $ħ\omega <2{\epsilon }_F$ ( $ħ\omega <2\mu$)

$\aleph \left(\frac{ħ\omega +2{\epsilon }_F}{ħv_{F}q}\right)-\aleph \left(\frac{2{\epsilon }_F-ħ\omega }{ħv_{F}q}\right)=\frac{8ħ\omega {\epsilon }_F}{{ħ}^2{v}_F^2{q}^2}+2\ln \left|\frac{2{\epsilon }_F-ħ\omega }{2{\epsilon }_F+ħ\omega }\right|.$E67

As a result, $\text{Δ}{\chi }^0(q,\omega )$ reduces to

$\begin{array}{c}\text{Δ}{\chi }^0(q,\omega )=-\frac{g{\epsilon }_F}{2\pi {ħ}^2{v}_F^2}+\frac{\mathfrak{I}(q,\omega )}{{ħ}^2{v}_F^2}\left\{\frac{8ħ\omega {\epsilon }_F}{{ħ}^2{v}_F^2{q}^2}+2\ln \left|\frac{2{\epsilon }_F-ħ\omega }{2{\epsilon }_F+ħ\omega }\right|+i\pi \right\}\hfill \\ =\frac{g{q}^2}{8\pi ħ\omega }\left\{\frac{2{\epsilon }_F}{ħ\omega }+\frac{1}{2}\ln \left|\frac{2{\epsilon }_F-ħ\omega }{2{\epsilon }_F+ħ\omega }\right|+\frac{i\pi }{2}\right\}.\hfill \end{array}$E68

If $2{\epsilon }_F\gg ħ\omega$, then

$\text{Δ}{\chi }^0(q,\omega )=\frac{g{q}^2}{8\pi ħ\omega }\left\{\frac{2{\epsilon }_F}{ħ\omega }+\frac{i\pi }{2}\right\}.$E69

By taking the decay rate $\omega \to \omega +i{\tau }^{-1}$ into account and inserting Eq. (62) into Eq. (61), the dynamical polarizability reduces to

${\chi }^0(q,\omega )\approx \frac{{\epsilon }_F{q}^2}{\pi {ħ}^2{\left(\omega +i{\tau }^{-1}\right)}^2}.$E70

In the presence of optical phonons, the effective dielectric function in the RPA regime is given by Eq. (4), which we restate here for convenience [98, 99]:

${\epsilon }^{\text{RPA}}(q,\omega )={\epsilon }_m-v_c(q){\chi }^0(q,\omega )-{\epsilon }_m{\displaystyle {\displaystyle \sum _l}{v}_{s\text{ph},l}(q,\omega )}{\chi }^0(q,\omega )-{\epsilon }_m{v}_{\text{oph}}(q,\omega ){\chi }_{j,j}^0(q,\omega ),$E71

where ${\epsilon }_m=({\epsilon }_1+{\epsilon }_2)/2$ is the average of dielectric constants of graphene’s environment. The collective oscillation modes of the electrons can be obtained by solving ${\epsilon }^{\text{RPA}}(q,\omega )=0$. The extinction function is identified by $Z=-\frac{\delta T}{T_0}$, or for the plasmonic structure coupled to an optical cavity $Z=-\frac{\delta R}{R_0}$, where $\delta R=R-R_0$ and R/R₀ is the reflectance with/without plasmon excitation, which corresponds to the enhanced absorbance at resonance frequencies [129, 131]

$Z\propto -\mathrm{Im}\left(\frac{1}{{\epsilon }^{\text{RPA}}}\right).$E72

The loss function represents the amount of energy dissipated by exciting the plasmon coupled to the substrate and graphene optical phonons. The collective oscillation modes of the electrons can be obtained by solving ${\epsilon }^{\text{RPA}}(q,\omega )=0$. Considering the first two terms in Eq. (71) gives the plasmon dispersion relation of graphene

$ħ{\omega }_p(q)=\sqrt{\frac{g{\alpha }_G{\epsilon }_Fħv_{F}q}{2{\epsilon }_0}},$E73

where the graphene fine structure constant is given by

${\alpha }_G\equiv \frac{e^2}{\epsilon ħv_F}=\frac{c}{\epsilon v_F}\alpha ,$E74

with $\alpha =e^2/ħc$ being the fine structure constant. In the case of graphene on a SiO₂ substrate, ${\epsilon }_{r1}=1$ for air and ${\epsilon }_{r2}=3.9$ for SiO₂, which yield ${\alpha }_G=0.9$. For suspended graphene, ${\epsilon }_{r1}={\epsilon }_{r2}=1$ and therefore ${\alpha }_G=2.2$. Thus, it is possible to tune α_G by altering the dielectric materials surrounding the graphene sheet [141].

It is also interesting to note that in contrast to the ${\epsilon }_F{}^{3/4}\propto n_0^{1/2}$ dependence of the bulk plasmon resonance frequency in Eq. (20) and the ${\epsilon }_F{}^{1/2}\propto n_s^{1/2}$ dependence of the 2DEG plasmon, the graphene plasmon exhibits a ${\epsilon }_F^{1/2}\propto n_s^{1/4}$ dependence.

In order to show the results for ${\epsilon }_F>0$, we define two regimes for $\aleph \left(x\right)$, i.e.,

$\begin{array}{cc}{\aleph }_{>}\left(x\right)=x\sqrt{x^2-1}-{\cosh }^{-1}(x),& \text{for }x>1\\ {\aleph }_{<}\left(x\right)=x\sqrt{1-x^2}-{\cosh }^{-1}(x),& \text{for }x<1\end{array}$E75

Using Eq. (60), the imaginary part of dynamical polarization is given by [140]

$\mathrm{Im}\text{Δ}{\chi }^0(q,\omega )=\frac{\mathfrak{I}(q,\omega )}{{ħ}^2{v}_F^2}\times \{\begin{array}{l}{\aleph }_{>}(\frac{2{\epsilon }_F-ħ\omega }{ħv_{F}q})-{\aleph }_{>}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q}),\to 1A\\ \pi \to 1B\\ -{\aleph }_{>}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q}),\to 2A\\ -{\aleph }_{<}(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q}),\to 2B\\ 0,\to 3A\\ 0.\to 3B\end{array}$E76

and the real part is obtained by the Kramers-Kronig relation, yielding

$\mathrm{Re}\text{Δ}{\chi }^0(q,\omega )=-\frac{g{\epsilon }_F}{2\pi {ħ}^2{v}_F^2}+\frac{\mathfrak{I}(q,\omega )}{{ħ}^2{v}_F^2}\times \{\begin{array}{l}\pi ,\to 1A\\ -{\aleph }_{>}(\frac{2{\epsilon }_F-ħ\omega }{ħv_{F}q})+{\aleph }_{>}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q}),\to 1B\\ -{\aleph }_{<}(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q}),\to 2A\\ {\aleph }_{>}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q}),\to 2B\\ -{\aleph }_{<}(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q})+{\aleph }_{<}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q}),\to 3A\\ {\aleph }_{>}(\frac{2{\epsilon }_F+ħ\omega }{ħv_{F}q})-{\aleph }_{>}(\frac{ħ\omega -2{\epsilon }_F}{ħv_{F}q}).\to 3B\end{array}$E77

Different regions are shown in Figure 13. As we mentioned, the plasmon dispersion relation is determined by $\epsilon (q,{\omega }_p-i\gamma )=0$, where γ is the decay rate of the plasmon. For weak damping, the plasmon dispersion relation and damping are given by

Figure 13.

Different regions related to the dynamical polarization. The regions are separated by the straight line $\omega =v_{F}q$ (solid), $ħ\omega =ħv_{F}q-2{\epsilon }_F$ (dashed), and $ħ\omega =2{\epsilon }_F-ħv_{F}q$ (dotted). The red-shaded area depicts the region of interband Landau damping while the green-shaded area marks the region of intraband Landau damping. The radiative interband Landau damping originating from direct (vertical) transitions occurs at q = 0. The red-shaded area with $q\ne 0$ corresponds to the region of nonradiative interband Landau damping.

${\epsilon }_0=v_c\mathrm{Re}\text{Δ}{\chi }^0(q,{\omega }_p),\gamma =\frac{\mathrm{Im}\text{Δ}{\chi }^0(q,{\omega }_p)}{{(\partial /\partial \omega )\mathrm{Re}\text{Δ}{\chi }^0(q,\omega )|}_{{\omega }_p}}.$E78

The solution of the first part of Eq. (78) exists only for $\mathrm{Re}\text{Δ}{\chi }^0(q,{\omega }_p)>0$, which is valid only for finite graphene doping and $\omega >v_{F}q$. Interestingly, a plasmon does not decay if $\mathrm{Im}\text{Δ}{\chi }^0(q,{\omega }_p)=0$, which is the case in region 1B shown in Figure 13. For the finite doping $\left({\epsilon }_F>0\right)$, the acoustic phonon at long wavelength is inside the region 1A in Figure 13.

7. Surface plasmon resonances in 3D topological insulators

So far, we have discussed nanoplasmonics in 3D metals and graphene and their damping mechanisms. Now we focus on the RPA theory of nanoplasmonics in 3D topological insulator (TI) materials. In particular, we are going to identify the graphene plasmons as a limiting case of 3DTI plasmons in the case when the thickness of the 3DTI slab becomes atomically thin. Bismuth selenide (Bi₂Se₃) is the prime example of a 3DTI material that has a rhombohedral crystal structure and consists of five atomic layers arranged along the z-direction, known as quintuple layers (QLs). Two such QLs are coupled by a weak interaction, predominantly of the van der Waals type. At the vicinity of the Γ point in the Brillouin zone of Bi₂Se₃, the low energy of Dirac particles in the xy-plane can be described by the effective Hamiltonian [142, 143], $H_{TI}=ħv_f(\sigma \times k\cdot {\displaystyle \widehat{\mathbf{z}}})-E_F$, where z is the unit vector in a z-direction. In the case of graphene, it is the pseudospin, representing the two sub-lattice indices, which is locked to the momentum of the electron. The chiral nature of the electronic states with the spin locked in perpendicular direction to the momentum results in a unique type of collective excitations which are, indeed, accompanied by transverse spin fluctuations [142, 143].

Let us consider a slab of a 3DTI material of thickness d > 5 nm, which is large enough to suppress any overlap of the single-electron states between the top and the bottom layers. In contrast, long-range Coulomb interaction exists and couples the opposite surfaces as in 2D electron plasmas.

We consider a Hamiltonian that describes the properties of collective oscillations on the top and the bottom surfaces of a 3DTI material as [142]

$H^{l/l}{}^{\text{'}}=ħv_f{\displaystyle {\displaystyle \sum _{\alpha ,\beta ,k,l,l^{\text{'}}}}{a}_{k,l,\alpha }^{+}\left(\tau \times k\right)·{\displaystyle \widehat{z}}{a}_{k,l,\beta }}+\frac{1}{2S}{\displaystyle {\displaystyle \sum _{q,l,l^{\text{'}}}}{\upsilon }_{l,l^{\text{'}}}}\left(q\right)n_{q,l}{n}_{-q,l^{\text{'}}},$E79

where $l\left(l^{\text{'}}\right)$ indicates the top (bottom) surface, $\alpha ,\beta$ are spin indices, a is the annihilation operator, and S is the surface area. $n_{q,l}={\displaystyle {\displaystyle \sum }_{k,\alpha }}{a}_{k-q,l,\alpha }^{+}{a}_{k,l,\alpha }$ is the density operator, and $\tau =\left\{{\tau }_x,\text{}{\tau }_y\right\}$ are the Pauli matrices. With a long-range Coulomb interaction within the RPA, the response function can be written as

${\chi }_{\text{RPA}}^{2x2}=\frac{{\chi }_0^{2x2}\left(q,\omega \right)}{1-\upsilon \left(q\right){\chi }_0^{2x2}\left(q,\omega \right)},$E80

where ${\chi }_0^{2x2}(q,\omega )=\left(\begin{array}{cc}{\chi }^{11}& 0\\ 0& {\chi }^{22}\end{array}\right)$ is a 2 × 2 matrix. $\upsilon \left(q\right)$ is a 2 × 2 matrix whose diagonal and off diagonal elements are the two-dimensional Fourier transforms of Coulomb potentials and can be obtained by solving Poisson’s equation [144]. Similar to the RPA dielectric function for graphene shown in Eq. (71), the equivalent equation for the RPA dielectric function in 3DTIs is

${\epsilon }_{\text{RPA}}^{2x2}\left(q,\omega \right)=1-v\left(q\right){\chi }_0^{2x2}\left(q,\omega \right).$E81

Compared to graphene, the main difference is that $\upsilon \left(q\right)$ is a 2 × 2 matrix accounting for the intrasurface and intersurface interactions. Here we write these expressions as

${\upsilon }^{11/22}\left(q\right)=\frac{4\pi e^2{\epsilon }_2\left(\cosh \left(qd\right)+\left({\epsilon }_{3/1}/{\epsilon }_2\right)\sinh \left(qd\right)\right)}{q\left({\epsilon }_2\left({\epsilon }_1+{\epsilon }_3\right)\cosh \left(qd\right)+\left({\epsilon }_1{\epsilon }_3+{\epsilon }_2\right)\sinh \left(qd\right)\right)}$E82

for the intrasurface interaction, and

${\upsilon }^{12}\left(q\right)=\frac{4\pi e^2{\epsilon }_2}{q\left({\epsilon }_2\left({\epsilon }_1+{\epsilon }_3\right)\cosh \left(qd\right)+\left({\epsilon }_1{\epsilon }_3+{\epsilon }_2\right)\sinh \left(qd\right)\right)}$E83

for the intersurface interaction. ε_i (i = 1, 2, 3) are the dielectric functions for the top, middle, and bottom layers, respectively. In the limit of equal dielectric constant, ${\epsilon }_1={\epsilon }_2={\epsilon }_3=\epsilon$, $V_{11}=V_{22}\to 2{\pi }^2{e}^2/\epsilon q$, and $V_{12}=V_{21}\to \left(2{\pi }^2{e}^2/\epsilon q\right)e^{-qd}$. This shows that for a given wave vector the intrasurface Coulomb interaction depends on the dielectric constant of the materials, the intersurface Coulomb interaction not only depends on the dielectric constant, but is also modulated by the negative distance-dependent exponential factor. For a sufficiently thick slab, i.e., $qd\to \infty$, the top and the bottom surfaces are completely decoupled and therefore we obtain the limit of $V_{12}=V_{21}\to 0$, and thus we only have a nonzero intrasurface potential. This is the case that is equivalent to the Coulomb interaction in a sheet of graphene.

The response function χ provides important information about the collective states that are excited at small transferred momentum. The collective mode frequencies of the system can be obtained by solving $\text{Det}[1-\upsilon \left(q\right){\chi }_0^{2x2}\left(q\right)]=0$. In the region $v_{f}q<\omega <2E_f/ħ-v_{f}q$, collective modes of oscillations are undamped. Beyond that limit, such modes are not observed because the energy of the modes is transferred to the particle-hole excitations. In Section 6, we derived the dynamical polarizability for the case of a graphene sheet. For $v_{f}q<\omega$ with $\mathrm{Im}[\chi ]=0$, the linear response function can be written as [144]

${\chi }_{\text{nn}}^{l/l^{\text{'}}}\left(q,\omega \right)=-\frac{E_F^{l/l}{}^{\text{'}}}{2\pi {ħ}^2{\upsilon }_F{}^2}+\frac{1}{16\pi ħ}\frac{q^2}{\sqrt{{\omega }^2-{\left({\upsilon }_{F}q\right)}^2}}\times \left[G\left(\frac{2E_F^{l/l}{}^{\text{'}}+ħ\omega }{ħ{\upsilon }_{F}q}\right)-G\left(\frac{2E_F^{l/l}{}^{\text{'}}-ħ\omega }{ħ{\upsilon }_{F}q}\right)\right],$E84

where $G\left(x\right)=x\sqrt{x^2-1}-{\cosh }^{-1}\left(x\right)$. Using the relation for the relative permittivity ${\epsilon }_r=1+\chi$, the linear response function given by Eq. (84) is obtained from Eq. (61) in the long wavelength limit. This charge-charge response function describes explicitly the dispersion of the collective modes exist on the two sides of a slab geometry of a 3DTI slab. Using Eq. (80), the pole of $\text{Det}[1-\upsilon \left(q\right){\chi }_0^{2x2}\left(q\right)]=0$ can be solved with the help of Eq. (35). The potential υ(q) is a 2×2 block diagonal (corresponding to spin and charge) for two surfaces. In the limit of an infinitely thick slab, we recover the intrasurface potential for $\upsilon \left(q\right)$ and the corresponding mode for graphene. The two modes in the long wavelength limit are: the symmetric (photon-like) mode and the antisymmetric (plasmon-like) mode. At a thickness below the Thomas-Fermi screening length, the symmetric mode is highly Landau damped. In Figure 14, we show the antisymmetric mode for a slab of thickness 120 nm at different Fermi energies of 0.3, 0.2, and 0.1 eV, taking into account the inter and intrasurface potentials given in Eqs. (82) and (83) in the q → 0 limit of the RPA response function. This result can be compared with the experimental observation of the Dirac plasmon dispersion in Bi₂Se₃ [145].

Figure 14.

Plasmon dispersion in Bi₂Se₃ slab of thickness 120 nm at different Fermi energy with top and bottom dielectric constant of 1 and 9.3, respectively, taking full potential and $q\to 0$ limit. Only the antisymmetric mode is shown here.

In the limit of small thickness of the slab, when the condition qd → 0 is satisfied, the antisymmetric photon-like mode behaves as

${\omega }_{\text{ant}}^2=\frac{v_f^2{k}_f{\alpha }_D}{{\epsilon }_1+{\epsilon }_3}\left[{\left(\frac{1+\xi }{2}\right)}^{1/2}+{\left(\frac{1-\xi }{2}\right)}^{1/2}\right]q$E85

and the symmetric plasmon-like mode as

${\omega }_{\text{sym}}^2=\frac{v_f^2{k}_f{\alpha }_D}{\sqrt{2}{\epsilon }_2}\left[\frac{{\left(1-\xi \right)}^{1/2}{\left(1+\xi \right)}^{1/2}}{{\left(1-\xi \right)}^{1/2}+{\left(1+\xi \right)}^{1/2}}\right]d{q}^2,$E86

where $\xi =\left(n_1-n_2\right)/\left(n_1+n_2\right)$ is the density polarization and $k_f=\sqrt{4\pi (n_1+n_2)}$ is the total Fermi wave vector. ${\alpha }_D=e^2/ħv_f$ is the fine structure constant for the Dirac system. It is to be noted that indices 1 and 2 indicate two different surfaces. For equal Fermi energy at the top and the bottom surfaces, the density polarization reduces to zero, resulting in an equilibrium situation. When the thickness of the slab tends to zero, we recover only the antisymmetric mode, which is the mode obtained for a sheet of graphene. As the slab thickness goes to infinity, the two surfaces interact weakly and the intersurface potential falls rapidly as the thickness increases. In this case, we obtain that the two antisymmetric modes each correspond to a single sheet of graphene. Using a series expansion of the frequency, it is possible to obtain a more accurate solution for the symmetric mode, which is derived elsewhere [144]. The solution is given by

${\omega }_{\text{sym}}=\left(\sqrt{2}{\epsilon }_2+{\alpha }_D{k}_{f}d)v_{f}q/(\sqrt{2{\epsilon }_2\left({\epsilon }_2+\sqrt{2}{\alpha }_D{k}_{f}d\right)}\right)$E87

For a size of the slab smaller than the mode wavelength, i.e., $d<{\lambda }_p$, the antisymmetric mode is independent of thickness, and the symmetric mode has a d² dependence. The plasmon dispersions for the limit of qd → 0 for symmetric and antisymmetric modes with 10 and 20 nm of slab width at different dielectric environments are shown in Figure 15a. The electric field profiles along (E_x) and perpendicular (E_z) to the plasmon polarization direction are shown in Figure 15b for a thin slab of thickness 10 nm with ${\epsilon }_1={\epsilon }_3=3.8$ for an antisymmetric mode. We also show the Landau damping regions with a blue shade given by the line segments $\omega \le {\upsilon }_{F}q$ and $\omega \ge 2E_F-{\upsilon }_{F}q$, in complete analogy to Figure 13. The distribution of excited “hot” carriers is limited to a narrow region around the Fermi level. In quasi-equilibrium, the hot carrier distribution can affect the plasmon lifetime.

Figure 15.

(a) Symmetric (dotted) and antisymmetric (solid) modes of plasmon dispersion in 3DTI in the limit of $qd\to 0$ for the slab thickness of 10 and 20 nm (Eqs. (7) and (9)). The Landau damping region is shown with the blue shade, determined by $\omega \le {\upsilon }_{F}q$ and $\omega \ge 2E_F-{\upsilon }_{F}q$. (b) The electric field profile of surface plasmons along and perpendicular direction of polarization for antisymmetric mode (plasmon-like mode) in a slab of thickness 10 nm.

The topological surface states are extended inside the bulk with a localization length given by $l_o=ħv_f/E_g$ where E_g is the bulk bandgap. l_o is equal to 1.5 nm (~2 quintuple layer (QL)) in Bi₂Se₃. Due to the quantum confinement of the carriers, two different conductivities arise: one due to the topological surface states that extend to 2 QL and the other one due to the 2DEG trapped below the surface, which extends up to 8 QL [146]. Theoretically, it should be possible to separate out the two different effects and write the response functions as $\chi ={\chi }^{TI}+{\chi }^{\text{2DEG}}$ for two independent channels. The contribution of the 2DEG to the plasmon dispersion can be calculated separately [147].

Surface plasmons in 3DTI surfaces do not consist only of charge density waves but are also accompanied by spin density waves due to inherent nature of spin-momentum locking. This can be qualitatively understood by calculating the surface current in terms of the spin and charge quantity by $J=v_f\sigma \times {\displaystyle \widehat{\mathbf{z}}}$. The continuity equation $\partial n_q/\partial t=-\nabla \cdot J$ shows that the charge density n_q oscillates in phase with the transverse spin $s^T={\displaystyle \widehat{\mathbf{z}}}\cdot ({\displaystyle \widehat{\mathbf{q}}}\times s)$, and the amplitude of their ratio is given by $s^T/n_q=1/ħv_f{[E_F{\alpha }_D{v}_f/2q]}^{1/2}$, which can be much larger than unity in the long wavelength limit [142]. To understand the electromagnetic response of spin and charge oscillations in a selfsustained mode, we write the induced potential $V_{\text{ind}}=\upsilon \left(q\right){\chi }_0^{4x4}\left(q,\omega \right)V_{\text{ind}}$ under vanishing external perturbation [147]. One can solve again for the pole $\left(1-\upsilon \left(q\right){\chi }_0^{4x4}\left(q,\omega \right)\right)V_{{}^{\text{ind}}}=0$ in order to obtain the explicit expression for the induced potential. This time, both the spin and charge part of the response functions are included, i.e., the solutions include the explicit information of the spin and charge oscillations. The solutions to the induced potential proportional to

$V^{\pm }{}_{\text{ind}}\propto \left(\begin{array}{c}1\\ 0\\ y^{\pm }\\ 0\end{array}\right),$E88

where $y^{+}=1\text{ and }{y}^{-}=-\frac{E_F^l}{E_F^{l^{\text{'}}}}$. The induced charge density is calculated from the response function $\rho ={\chi }_0^{4x4}\left(q,\omega \right)V_{\text{ind}}$. Spins are inherently attached to the momentum as accounted by the complete response function given by Eq. (80). The symmetric and antisymmetric modes are purely spin-like and purely charge-like for a slab of 3DTI. This can be seen in the following expressions in the limit of $qd\to 0$ as

$\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right)=\rho \left(\begin{array}{c}1\\ -\sqrt{\frac{{\alpha }_D\left(E_F^l+E_F^{l^{\text{'}}}\right)}{ħ{\upsilon }_{F}q\left({\epsilon }^l+{\epsilon }^{l^{\text{'}}}\right)}}\\ \frac{E_F^l}{E_F^{l^{\text{'}}}}\\ \sqrt{\frac{{\alpha }_D\left(E_F^l+E_F^{l^{\text{'}}}\right)}{ħ{\upsilon }_{F}q\left({\epsilon }^l+{\epsilon }^{l^{\text{'}}}\right)}}\frac{E_F^l}{E_F^{l^{\text{'}}}}\end{array}\right)$E89

for the antisymmetric mode, and

$\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right)=\rho \left(\begin{array}{c}1\\ -\sqrt{\frac{{\alpha }_{D}d\left(E_F^l\times E_F^{l^{\text{'}}}\right)}{ħ{\upsilon }_{F}q{\epsilon }_{TI}\left(E_F^l+E_F^{l^{\text{'}}}\right)}}\\ -1\\ -\sqrt{\frac{{\alpha }_{D}d\left(E_F^l\times E_F^{l^{\text{'}}}\right)}{ħ{\upsilon }_{F}q{\epsilon }_{TI}\left(E_F^l+E_F^{l^{\text{'}}}\right)}}\end{array}\right)$E90

for the symmetric mode, where ${\phi }_1=\left(\begin{array}{c}{\rho }^l\\ s^l\end{array}\right)$ and ${\phi }_2=\left(\begin{array}{c}{\rho }^{l^{\text{'}}}\\ s^{l^{\text{'}}}\end{array}\right)$, and the indices l and l′ can be chosen to represent the top and the bottom surfaces. Clearly, charges oscillate in phase and spins oscillate out of phase for antisymmetric mode, whereas spins oscillate in phase and charges oscillate out of phase for the symmetric mode on the top and the bottom surfaces.

It has been observed that the phonon modes in Bi₂Se₃ have frequencies of around 2 THz. Since the phonon modes overlap with the plasmon modes in energy space, there is a possibility of mode hybridization, similar to the case of graphene. In a micro-ring structure of Bi₂Se₃, there are bonding (lower) and antibonding (upper) plasmon modes [148]. The antibonding plasmon mode frequency is much larger than the phonon mode frequency, and the bonding plasmon mode overlaps with the phonon mode, resulting in hybridization in two branches with an interaction frequency of around 0.35 THz.

Recently, in Ref. [149], authors investigated the possibility of obtaining SPPs in Weyl semimetals. In Weyl semimetals such as TaAs, NbAs, YbMnBi₂, and Eu₂Ir₂O₇ [150–152] the valence and the conduction band touch in isolated points of the Brillouin zone close to the chemical potential, and their dispersion is described by an equation similar to the one for Dirac metals. A pair of Weyl nodes appears with opposite chiralities with a distance of b in the reciprocal space. Their topological properties are described by $S_{\theta }=\frac{e^2}{4\pi ħc}{\displaystyle \int dt{\displaystyle \int d^{3}r\theta E\text{.}B}}$, where $\theta =2\left(br-b_{o}t\right)$. $S_{\theta }$ is an important parameter that alters the electrodynamics of the bulk response to an applied electric field. In such a system, the authors found that the SPPs with $b\ne 0$ have properties similar to those in ordinary metal SPPs in the presence of an external magnetic field, despite their unique origin. Therefore, the θ term that contributes to the dielectric tensor originating from the anomalous Hall displacement current induces an anomalous surface magnetoplasmon. Interestingly, for certain orientations of the surface of the Weyl semimetal, the authors predict the nonreciprocal dispersion of SPPs, i.e., the propagation of the SPP is unidirectional.

8. Conclusion

In this book chapter, we present descriptions of the plasmonic properties of metal nanostructures of different geometry, their size dependence, and applications in modern nanotechnology. We show dynamic control over the plasmon resonances where a shift up to 125 nm at a density of 10²² cm⁻³ can be achieved using a pump-probe technique. This provides the opportunity to utilize plasmonics in modern information processing devices. In addition to plasmonics in 3D metal nanostructures, we present a description of graphene and 3DTI plasmonics using classical and quantum perspectives. Using the RPA theory, we obtain symmetric and antisymmetric modes in a slab of 3DTI, which reduces to a graphene plasmonics in the limit of zero thickness of the slab. Surface plasmon (SP) damping mechanisms are interesting due to their potential applications for enhanced current density that comes from SPs nonradiative damping in nanostructures smaller than the skin depth. We present a quantum theory of SPs damping in metals and layered materials like graphene.

There are several potential applications of the nanoplasmonics in graphene and 3DTI. Graphene and 3DTI are potential candidates for nanospasers that utilize Dirac fermions, unlike the massive electrons or holes in the originally proposed spasing scheme by Bergman and Stockman in 2003 [153]. The spaser is a nanoplasmonic counterpart of a laser, where photons are not emitted. In Ref. [154], the authors have proposed a scheme of nanospasing using a sheet of graphene with an electrically pumped cascaded quantum well structure working as a gain medium. For the range of Fermi energy, it was shown that spasing could be potentially obtained for a plasmon relaxation time of several femtoseconds in the mid-IR range. A similar scheme with an optically pumped nanospaser for a slab 10 nm of 3DTI was proposed in Ref. [155]. It is advantageous to use a 3DTI, such as Bi₂S₃, as a nanospaser due to the possibility of using its bulk as a gain medium and the surface as a medium that supports SPPs. This configuration avoids the use of a separate gain medium to provide the feedback for the SPPs. Therefore, a 3DTI nanospaser can be truly nanoscopic and can be used for various applications in physics, chemistry, and biology.

Acknowledgments

We acknowledge the support provided by NSF grant CCF-1514089 and DARPA grant HR0011-16-1-0003.