Quantization of graphene plasmons

TL;DR

This paper presents a comprehensive quantization of graphene plasmons using classical and quantum hydrodynamic models, accounting for dispersion and non-local effects, and analyzes their impact on emitter decay rates.

Contribution

It introduces a novel quantization method for graphene plasmons considering dispersion and non-local effects, with applications to emitter decay rate calculations.

Findings

Mode-functions depend on graphene's dispersive properties

Non-local effects significantly influence decay rates at small Fermi energies

Quantization methods accurately predict emitter decay near graphene

Abstract

In this article we perform the quantization of graphene plasmons using both a macroscopic approach based on the classical average electromagnetic energy and a quantum hydrodynamic model, in which graphene charge carriers are modeled as a charged fluid. Both models allow to take into account the dispersion of graphenes optical response, with the hydrodynamic model also allowing for the inclusion of non-local effects. Using both methods, the electromagnetic field mode-functions, and the respective frequencies, are determined for two different graphene structures. we show how to quantize graphene plasmons, considering that graphene is a dispersive medium, and taking into account both local and nonlocal descriptions. It is found that the dispersion of graphene's optical response leads to a non-trivial normalization condition for the mode-functions. The obtained mode-functions are then used to calculate the decay of an emitter, represented by a dipole, via the excitation of graphene surface plasmon-polaritons. The obtained results are compared with the total spontaneous decay rate of the emitter and a near perfect match is found in the relevant spectral range. It is found that non-local effects in graphene's conductivity, become relevant for the emission rate for small Fermi energies and small distances between the dipole and the graphene sheet.

Equations

$$\mathbf{E}(\mathbf{r},t)$$

$$\mathbf{B}(\mathbf{r},t)$$

$$\begin{gathered} U_{\text{EM}}(\omega )=\int d^3\mathbf{r}({ϵ}_0{\mathbf{E}}_{\omega }^{\ast }(\mathbf{r})\cdot \frac{\partial }{\partial \omega }\left[\omega {\overline{ϵ}}_r(\mathbf{r},\omega )\right]\cdot {\mathbf{E}}_{\omega }(\mathbf{r}) \\ +\frac{1}{{\mu }_0}{\left|{\mathbf{B}}_{\omega }\right|}^2), \end{gathered}$$

$$\mathbf{E}(\mathbf{r},t)$$

$$\mathbf{B}(\mathbf{r},t)$$

$$\mathbf{A}(\mathbf{r},t)=\sum _{\lambda }{\alpha }_{\lambda }{e}^{-i{\omega }_{\lambda }t}{\mathbf{A}}_{\lambda }(\mathbf{r})+\text{c.c.},$$

$$\nabla \times \nabla \times {\mathbf{A}}_{\lambda }(\mathbf{r})=\frac{{\omega }_{\lambda }^2}{c^2}{\overline{ϵ}}_r(\mathbf{r},{\omega }_{\lambda })\cdot {\mathbf{A}}_{\lambda }(\mathbf{r}),$$

$$U_{\text{EM}}=\sum _{\lambda }{U}_{\text{EM}}({\omega }_{\lambda }){\left|{\alpha }_{\lambda }\right|}^2.$$

$${\alpha }_{\lambda }$$

$${\alpha }_{\lambda }^{†}$$

$$\left[{\hat{a}}_{\lambda },{\hat{a}}_{{\lambda }^{\prime }}^{†}\right]={\delta }_{\lambda ,{\lambda }^{\prime }},$$

$$\hat{H}=\sum _{\lambda }\frac{\hslash U({\omega }_{\lambda })}{4L_{\lambda }{ϵ}_0{\omega }_{\lambda }}\left({\hat{a}}_{\lambda }^{†}{\hat{a}}_{\lambda }+{\hat{a}}_{\lambda }{\hat{a}}_{\lambda }^{†}\right),$$

$$\hat{H}=\frac{1}{2}\sum _{\lambda }\hslash {\omega }_{\lambda }\left({\hat{a}}_{\lambda }^{†}{\hat{a}}_{\lambda }+{\hat{a}}_{\lambda }{\hat{a}}_{\lambda }^{†}\right).$$

$$\begin{gathered} \int d^3\mathbf{r}{\left|\nabla \times {\mathbf{A}}_{\lambda }(\mathbf{r})\right|}^2=\int d^3\mathbf{r}{\mathbf{A}}_{\lambda }^{\ast }(\mathbf{r})\cdot \nabla \times \nabla \times {\mathbf{A}}_{\lambda }(\mathbf{r}) \\ =\frac{{\omega }_{\lambda }^2}{c^2}\int d^3\mathbf{r}{\mathbf{A}}_{\lambda }^{\ast }(\mathbf{r})\cdot {\overline{ϵ}}_r(\mathbf{r},{\omega }_{\lambda })\cdot {\mathbf{A}}_{\lambda }(\mathbf{r}), \end{gathered}$$

$$L_{\lambda }=\int d^3\mathbf{r}{\mathbf{A}}_{\lambda }^{\ast }(\mathbf{r})\cdot \left({\overline{ϵ}}_r(\mathbf{r},{\omega }_{\lambda })+\frac{{\omega }_{\lambda }}{2}\frac{\partial }{\partial \omega }{\overline{ϵ}}_r(\mathbf{r},{\omega }_{\lambda })\right)\cdot {\mathbf{A}}_{\lambda }(\mathbf{r}).$$

$$\mathcal{L}={\mathcal{L}}_{\text{EM}}+{\mathcal{L}}_{\text{2D hyd}},$$

$${\mathcal{L}}_{\text{EM}}=\frac{{ϵ}_0}{2}\left({\partial }_t\mathbf{A}\right)\cdot {\overline{ϵ}}_d\cdot \left({\partial }_t\mathbf{A}\right)-\frac{1}{2{\mu }_0}{\left(\nabla \times \mathbf{A}\right)}^2,$$

$$\mathcal{L}_{\text{2D hyd}}=\delta(z)\Biggl{(}\frac{1}{2}n_{0}m\left(\partial_{t}\bm{\upsilon}\right)^{2}+m\beta^{2}n\nabla\cdot\bm{\upsilon}\\ \left.+\frac{m\beta^{2}}{2n_{0}}n^{2}-n_{0}e\partial_{t}\bm{\upsilon}\cdot\mathbf{A}\right),$$

$$\nabla \times \mathbf{B}=\frac{1}{c^2}{\overline{ϵ}}_d{\partial }_t\mathbf{E}-{\mu }_0{n}_{0}e{\partial }_t\mathbf{\upsilon }\delta (z),$$

$${\mathbf{J}}_s^{\text{hyd}}=-e{n}_0{\partial }_t\mathbf{\upsilon }.$$

$$n_0\nabla \cdot \mathbf{\upsilon }+n$$

$$n_{0}m{\partial }_t^2\mathbf{\upsilon }+m{\beta }^2\nabla n$$

$${\rho }_s^{\text{hyd}}=-en.$$

$${\mathcal{L}}^{\prime }={\mathcal{L}}_{\text{EM}}+{\mathcal{L}}_{\text{2D hyd}}^{\prime },$$

$$\mathcal{L}_{\text{2D hyd}}^{\prime}=\delta(z)\Biggl{(}\frac{1}{2}n_{0}m\left(\partial_{t}\bm{\upsilon}\right)^{2}\\ -\frac{1}{2}n_{0}m\beta^{2}\left(\nabla\cdot\bm{\upsilon}\right)^{2}-n_{0}e\partial_{t}\bm{\upsilon}\cdot\mathbf{A}.\Biggr{)}$$

$$-m{\partial }_t^2\mathbf{\upsilon }+m{\beta }^2\nabla \left(\nabla \cdot \mathbf{\upsilon }\right)=e\mathbf{E}(z=0),$$

$$\mathbf{\Pi }$$

$$\mathbf{\pi }$$

$$\begin{gathered} H=\int d^3\mathbf{r}\left(\frac{1}{2{ϵ}_0}\mathbf{\Pi }\cdot {\overline{ϵ}}_d^{-1}\cdot \mathbf{\Pi }+\frac{1}{2{\mu }_0}{\left(\nabla \times \mathbf{A}\right)}^2\right) \\ +\int d^2\mathbf{x}\left(\frac{{\left(\mathbf{\pi }+e{n}_0\mathbf{A}(z=0)\right)}^2}{2n_{0}m}+\frac{1}{2}{n}_{0}m{\beta }^2{\left(\nabla \cdot \mathbf{\upsilon }\right)}^2\right). \end{gathered}$$

$$\mathbf{A}(\mathbf{r},t)$$

$$\mathbf{\upsilon }(\mathbf{x},t)$$

$$\begin{gathered} \left[{\omega }_{\mathbf{q},\lambda }^2{ϵ}_0{\overline{ϵ}}_d(z)-\frac{1}{{\mu }_0}{D}_{\mathbf{q}}\times D_{\mathbf{q}}\times \right]{\mathbf{A}}_{\mathbf{q},\lambda }(z)= \\ =-i{\omega }_{\mathbf{q},\lambda }\delta (z)e{n}_0{\mathbf{\upsilon }}_{\mathbf{q},\lambda }, \end{gathered}$$

$$m{n}_0\left[{\omega }_{\mathbf{q},\lambda }^2-{\beta }^2\mathbf{q}\otimes \mathbf{q}\right]{\mathbf{\upsilon }}_{\mathbf{q},\lambda }=i{\omega }_{\mathbf{q},\lambda }e{n}_0{\mathbf{A}}_{\mathbf{q},\lambda }(0).$$

Full text

I Introduction

In many cases, light-matter interaction can be understood in a semi-classical picture, where matter is quantized and the electromagnetic field (EM) is treated classically. This semi-classical approach holds when the number of photons in the field is large or the light source is coherent. On the other hand, in order to understand the properties of a small number of photons the quantization of the EM field is required. Typical phenomena where the quantization of the EM field is necessary involve entanglement, squeezed light, cavity electrodynamics, interaction of photons with nano-mechanical resonators, and near-field radiative effects (Agarwal, 2013).

In near-field radiative effects, plasmonics emerges as a promising candidate to observe quantum effects of the electromagnetic radiation, an example being the Hong-Ou-Mandel interference of plasmons (Heeres et al., 2013). Many other quantum effects in plasmonics exist, such as the survival of entanglement, particle-wave duality, quantum size effects due to reduced dimensions of metallic nanostructures, quantum tunneling of plasmons (which are simultaneously light and matter), and coupling of quantum emitters to surface plasmons (Jacob, 2012; Tame et al., 2013; Martín-Cano et al., 2014; Fakonas et al., 2015; Fitzgerald et al., 2016; Dheur et al., 2016; Bozhevolnyi and Mortensen, 2017; Xu et al., 2018; Fernández-Domínguez et al., 2018).

The exploration of quantum effects in plasmonics in unusual spectral ranges, such as the THz and the mid-IR, has been deterrent by the poor plasmonic response of noble metals in these regions of the electromagnetic spectrum. However, with the emergence of graphene plasmonics (Ju et al., 2011; Fei et al., 2011) the possibility of exploring quantum effects in these yet unexplored spectral regions is a possible prospect. Despite the fact that, at the time of writing, quantum effects involving graphene plasmons remain illusive, the fact that graphene plasmons are characterized by low losses (Fei et al., 2012; Chen et al., 2012; Woessner et al., 2015) boosts the above hope. In addition, graphene plasmons screened by a nearby metal (also called screened or acoustic plasmons) can be confined down to the atomic limit (Alcaraz Iranzo et al., 2018), which certainly opens the prospects of finding rich grounds for quantum plasmonics and nonlocal effects (Lundeberg et al., 2017; Dias et al., 2018). Indeed, the idea of developing quantum optics with plasmons has already a long history (Chang et al., 2006) and quantization of localized plasmons in metallic nanoparticles was recently performed (Alpeggiani and Andreani, 2014; Bordo, 2019).

The development of quantum theory of the electromagnetic field in the presence of dieletric media has a long history and several approaches have been developed (Hopfield, 1958; Alekseev and Nikitin, 1966; Agarwal, 1975a; Huttner and Barnett, 1992; Milonni, 1995; Gruner and Welsch, 1995, 1996; Matloob et al., 1995; Dung et al., 1998; Philbin, 2010; Hanson et al., 2015; Allameh et al., 2016; Sha et al., 2018). These methods are typically based either on the quantization of the macroscopic classic energy(Milonni, 1995), on the classical dyadic Green’s function for the electric field(Gruner and Welsch, 1995) or on the diagonalization of Hopfield-type Hamiltonians(Huttner and Barnett, 1992). The quantization of evanescent EM waves (Carniglia and Mandel, 1971; Carniglia et al., 1972) and of the EM field in the vicinity of a metal (Grossel et al., 1977) have also been considered in the past.

In this paper, we perform the quantization of graphene plasmons, obtaining the plasmonic electromagnetic field mode-functions and, importantly, their normalization, when losses are neglected. These mode-functions are then used to study the interaction of graphene plasmons with nearby quantum emitters and determine, in a very intuitive way, using Fermi’s golden rule the spontaneous decay rate of the emitter due to plasmon emission. Thee quantization of graphene plasmons is performed in two ways: (i) A macroscopic approach, which starts from the classical time-averaged energy of the electromagnetic field in a dielectric (Landau and Lifshitz, 1984; Alekseev and Nikitin, 1966; Agu, 1979; Milonni, 1995). This method allows for the inclusion of dispersion in the optical response of graphene. (ii) A hydrodynamic approach, where graphene charge carriers are described in terms of an electronic fluid, which couples to the eletromagnetic field (Chaves et al., 2017; Fetter, 1973). The hydrodynamic approach allows not only for the inclusion of dispersion, but also for the inclusion of non-local effects in the optical response of graphene.

The paper is organized as follows: in Section II, we present the general macroscopic approach for the quantization of the eletromagnetic field in dispersive, lossless media and a normalization condition for the mode-functions is determined. In Section III, we present the quantum hydrodynamic model or graphene. We will see that when non-local effects are neglected, the result of the macroscopic approach is recovered. In Section IV, the plasmon dispersion relations, mode-functions and mode-function normalization for a single graphene layer and for a graphene-metal structures are obtained. In Section V, we use the quantized plasmon fields to compute the decay rate of a quantum emitter due to the spontaneous emission of plasmons. The plasmon emission rate is compared with the total decay rate of the emitter, which is obtained from the complete dyadic Green’s function for the electric field. The role of non-local response of graphene is analyzed. Finally, we conclude with Section VI, commenting the obtained results and discussing future research paths. A set of appendices detailing the calculations is also presented.

II Macroscopic quantization of the plasmon electromagnetic field

In this section, we will describe how the plasmon fields can be quantized using the macroscopic classical energy of the electromagnetic field in a dispersive, lossless dielectric medium, a method first used in Refs. (Alekseev and Nikitin, 1966; Agu, 1979; Milonni, 1995). For electric and magnetic fields with a harmonic time dependence,

[TABLE]

close to a central frequency $\omega$, the time-averaged classical electromagnetic energy in the presence of a dispersive, lossless dieletric is given by (Landau and Lifshitz, 1984; Ruppin, 2002)

[TABLE]

where $\bar{\epsilon}_{r}(\mathbf{r},\omega)$ is the relative dielectric tensor of the medium, $\epsilon_{0}$ and $\mu_{0}$ are, respectively, vacuum permittivity and permeability. The idea of this method is to take the above equation as the quantum mechanical energy of a EM field eigenmode with frequency $\omega$. We will work in the Weyl gauge, in which the scalar potential is set to zero $\phi=0$, such that the electric and magnetic fields are obtained only from the vector potential $\mathbf{A}$ as

[TABLE]

The vector potential is then expanded in modes as

[TABLE]

where $\alpha_{\lambda}$ are amplitudes for the mode $\lambda$, with mode-function $\mathbf{A}_{\lambda}(\mathbf{r})$ and corresponding frequency $\omega_{\lambda}$. The mode-functions and frequencies are determined by solving the non-linear eigenvalue problem

[TABLE]

with $c$ vacuum’s speed of light, which is just Ampère’s law in the dielectric medium for the mode-function $\mathbf{A}_{\lambda}(\mathbf{r})$. Next, we assume that the total time-averaged energy for the vector potential Eq. (6) is given by

[TABLE]

The quantization of the theory is achieved by promoting the amplitudes $\alpha_{\lambda}$ to quantum mechanical operators

[TABLE]

where $\hat{a}_{\lambda}$ ($\hat{a}_{\lambda}^{\dagger}$) are bosonic annihilation (creation) operators, which obey the usual equal-time commutation relations

[TABLE]

and $L_{\lambda}$ is a mode-length which is determined by demanding that the quantum mechanical Hamiltonian which is obtained from Eq. (8) by performing the replacement of Eqs. (9) and (10),

[TABLE]

coincides with the Hamiltonian for a collection of quantum harmonic oscillators

[TABLE]

Imposing this condition, we have that the mode-length is given by $L_{\lambda}=U(\omega_{\lambda})/\left(2\epsilon_{0}\omega_{\lambda}^{2}\right)$. Using the mode-function equation (7) to write

[TABLE]

The mode-length can be written as

[TABLE]

Notice that in the absence of dispersion, the second term vanishes and $L_{\lambda}$ reduces to the usual norm in the presence of a position dependent dielectric tensor $\bar{\epsilon}_{r}(\mathbf{r},\omega_{\lambda})$.

Although this phenomenological method appears to be unjustified, it has actually been shown to be correct for the case in when the dielectric is modeled by a Lorentz oscillator (Agu, 1979). We will see in the next section, that Eq. 15 remains valid within a quantum hydrodynamic model of graphene, even when non-local effects are included in the optical response. As a matter of fact Eq. (15) for the mode-length remains valid for any linear optical medium (including effects of dispersion, non-locality, inhomogeneity and anisotropy) as long as losses can be neglected (Amorim and Chaves, 2019).

III plasmon quantization within a quantum hydrodynamic model

In this section, we will perform the quantization of graphene surface plasmons employing a hydrodynamic model. The hydrodynamic model treats both the electron gas and the electromagnetic field using a classical picture and provides a simple and elegant way of including non-local effects (Christensen, 2017). Non-local effects are taken into account by including a pressure term in the Boltzmann equation, that arises due to Pauli’s exclusion principle. A detailed derivation of the hydrodynamic model for graphene can be found in (Chaves et al., 2017; Fetter, 1973). To illustrate the method, we choose the simple case of a single graphene sheet located at the plane $z=0$, embedded by a static dielectric medium with relative dielectric constant $\bar{\epsilon}_{d}(\mathbf{r})$.

III.1 Classical hydrodynamic Lagrangean and Hamiltonian

The classical Lagrangean density for the hydrodynamic model of graphene can be written as:

[TABLE]

where $\mathcal{L}_{\text{EM}}$ is the Lagrangean density of the electromagnetic field and $\mathcal{L}_{\text{2D hyd}}$ describes the electronic fluid of graphene and its coupling to the electromagnetic field. Using once again the Weyl gauge, $\mathcal{L}_{\text{EM}}$ is given by

[TABLE]

where $\bar{\epsilon}_{d}(\mathbf{r})$ is allowed to be position dependent, but is frequency independent. Within the hydrodynamic model, the electronic fluid of the graphene layer is described by the fluctuation, $n$, of the density around the equilibrium density, $n_{0}$, and the displacement vector $\bm{\upsilon}=\left(\upsilon_{x},\upsilon_{y},0\right)$, which should not be confused with the velocity field. In the Weyl gauge, $\mathcal{L}_{\text{2D fluid}}$ is written as(Fetter, 1973; Chaves et al., 2017)

[TABLE]

where the $\delta$-function $\delta(z)$ restricts the electronic fluid to the $z=0$ 2D plane, $m$ is the Drude mass and $\beta$ appears from the relation between the degeneracy pressure and the electronic density and depends on the band dispersion for the carriers (see Ref. (Chaves et al., 2017)). In the approximation of the linear dispersion for graphene electrons, the hydrodynamic parameters are given by (Chaves et al., 2017): $n_{0}=k_{F}^{2}/\pi$, $m=\hbar k_{F}/v_{F}$ and $\beta^{2}=v_{F}^{2}/2$, where $k_{F}$ is the Fermi wavenumber and $v_{F}$ the Fermi velocity of graphene. Equation (18) is the 2D equivalent Lagrangian for the hydrodynamic model presented in (Nakamura, 1983).

Using the Euler-Lagrange equations for Eq. (18) with respect to $\mathbf{A}$, we obtain

[TABLE]

which is nothing more than Ampère’s law in the presence of a surface current given by

[TABLE]

Using the Euler-Lagrange equations for Eq. (18) with respect to the fluid variables $n$ and $\boldsymbol{\upsilon}$, we obtain the continuity and Newton’s second law with a diffusion term, which read respectively

[TABLE]

from which we recognize the fluid electronic surface density

[TABLE]

Equations (19)-(22) correspond to the linearized hydrodynamic model for graphene (Chaves et al., 2017) (see also (Lucas and Fong, 2018)).

Notice that Eq. (21) has no dynamics. Therefore, we can use it to eliminate the field $n$, thus obtaining a new Lagrangean density

[TABLE]

with

[TABLE]

This new Lagrangean is equivalent to the Eq. (16) as both lead to the same dynamics. Applying the Euler-Lagrange equations to Eq. (24) with respect to $\mathbf{A}$ leads to Eq. (19), while the equation obtained for $\bm{\upsilon}$ reads

[TABLE]

which can be obtained by combining Eqs. (21) and (22).

From the Lagrangean density Eq. (24), we define the canonical momenta conjugate to $\mathbf{A}$ and $\bm{\upsilon}$, respectively, as

[TABLE]

In terms of the variables $\mathbf{A}$, $\bm{\Pi}$, $\bm{\upsilon}$ and $\bm{\pi}$, the classical Hamiltonian obtained from Eq. (24) is given by

[TABLE]

III.2 Canonical quantization of hydrodynamic model

In order to quantize the classical Hamiltonian Eq. (29), we start by introducing the eigenmodes of the coupled electronic fluid

• electromagnetic field. Assuming, we have in-plane translational invariance, we write the vector potential and fluid displacement variables as

[TABLE]

where $S$ is the area of the graphene layer, $\alpha_{\mathbf{q},\lambda}(t)=\alpha_{\mathbf{q},\lambda}e^{-i\omega_{\mathbf{q},\lambda}t}$ are mode amplitudes with, $\omega_{\mathbf{q},\lambda}$the mode frequency, and $\mathbf{A}_{\mathbf{q},\lambda}(z)$ and $\bm{\upsilon}_{\mathbf{q},\lambda}$ are mode-functions, which, following from Eqs. (19) and (26), obey the equations

[TABLE]

[TABLE]

where we defined the differential operator $D_{\mathbf{q}}=i\mathbf{q}+\hat{\mathbf{z}}\partial_{z}$.

From Eq. (33), we can write

[TABLE]

where $\bar{\sigma}_{g}^{\text{hyd}}(\mathbf{q},\omega)$ is the conductivity within the hydrodynamic model, which we separate into transverse and longitudinal components as

[TABLE]

respectively given by

[TABLE]

where we identify $D=e^{2}n_{0}/m$ as the Drude weight, which for graphene is given by $D=e^{2}v_{F}k_{F}/\pi$. Notice that in the the limit $\mathbf{q}\rightarrow 0$, $\bar{\sigma}_{g}^{\text{hyd}}(\mathbf{q},\omega)$ reduces to the Drude model. Inserting Eq. (34) into Eq. (32), we obtain

[TABLE]

with the dieletric function, including screening by graphene electrons, being given by

[TABLE]

in agreement with Eq. (7).

Inserting the expansions Eq. (30) and (31) into Eq. (29), and using the orthogonality properties of the mode-functions $\left(\mathbf{A}_{\mathbf{q},\lambda}(z),\,\bm{\upsilon}_{\mathbf{q},\lambda}\right)$ it is possible to write the Hamiltonian for the hydrodynamic model as (see Appendix (B))

[TABLE]

with the mode-length defined here as

[TABLE]

Promoting the mode amplitudes to quantum mechanical operators as

[TABLE]

where \hat{a}_{\mathbf{q},\lambda}^{\dagger}$$\left(\hat{a}_{\mathbf{q},\lambda}\right) are creation (annihilation) operators for the plasmon-polaritons, satisfying the usual equal-time commutation relations Eq. (11), the quantum Hamiltonian for the hydrodynamic model becomes

[TABLE]

We will now see that Eq. (41) can be recast in the same form as Eq. (41). In order to do so, we use Eq. (34) which allows to write Eq. (41) as

[TABLE]

It is easy to see that the above equation can also be written as

[TABLE]

with $\bar{\epsilon}_{r}(\mathbf{q},z,\omega)$ given by Eq. (39).

IV Dispersion relations and mode-functions of graphene plasmon in two

different structures

We now wish to determine the dispersion relation and mode-functions of graphene plasmons in two different structures (see Fig. 1): a single graphene layer and a graphene sheet in the vicinity of a perfect metal. As in the previous section, we make use of the in-plane translational invariance of the structures being considered. Therefore, the non-linear eigenvalue problem for the mode-functions, Eq. 7, can be written as

[TABLE]

where

[TABLE]

where $\bar{\epsilon}_{d}\left(z\right)$ is the dielectric constant of the medium, which we assume to be isotropic and a piecewise homogeneous function of $z$, and $\ell$ labels the graphene layers which are located at the planes $z=z_{\ell}$, with conductivity$\bar{\sigma}_{g\ell}(\mathbf{q},\omega)$. We model the conductivity of each graphene layer with Eq. (35), which when including losses becomes

[TABLE]

where $D=4E_{F}\sigma_{0}/\left(\pi\hbar\right)$ is the Drude weight, with $E_{F}$ graphene’s Fermi energy, $\sigma_{0}=\pi e^{2}/\left(2h\right)$ and $\gamma$ is a relaxation rate. When determining mode-functions we will set $\gamma=0$, but we allow for $\gamma\neq 0$, for the situation when ohmic losses are included in Section (V). The presence of the $\delta$-functions in Eq. (48), implies that boundary conditions at each interface located at $z=z_{\ell}$:

[TABLE]

where $\mathbf{E}_{\mathbf{q},\lambda}\left(z\right)=i\omega_{\mathbf{q},\ell}\mathbf{A}_{\mathbf{q},\lambda}(z)$ and $\mathbf{B}_{\mathbf{q},\lambda}\left(z\right)=D_{\mathbf{q}}\times\mathbf{A}_{\mathbf{q},\lambda}(z)$ are the electric and magnetic fields corresponding to mode $\mathbf{A}_{\mathbf{q},\lambda}(z)$, and $\mathbf{J}_{s,\mathbf{q},\lambda}(z_{\ell})=\bar{\sigma}_{g\ell}(\mathbf{q},\omega)\cdot\mathbf{E}_{\mathbf{q},\lambda}\left(z_{\ell}\right)$ is the surface current in the graphene layer $\ell$. In addition to the boundary conditions Eqs. (51) and (52), to determine of the plasmon modes one must impose that the fields decay for $z\rightarrow\pm\infty$. Having determined the plasmon mode-function, $\mathbf{A}_{\mathbf{q},\text{sp}}(z)$, and dispersion, $\omega_{\mathbf{q},\text{sp}}$, the mode-length can be obtained from Eqs. (46) and (48) as

[TABLE]

We will now analyse the different structures in detail.

IV.1 Single layer graphene

We first discuss the simplest case of a single graphene sheet (see left panel of Fig. 1). The problem of finding the spectrum of the surface plasmons in a graphene sheet was first considered in (Jablan et al., 2009) and a detailed derivation can be found in Refs. (Bludov et al., 2013; Gonçalves and Peres, 2016). We assume that the single layer of graphene is located at $z=0$, with a encapsulating dielectric medium $n=2$ for $z>0$ and a medium $n=1$ for $z<0$, such that

[TABLE]

In order to determine the plasmon mode, we look for $p$-polarized solutions of the electric field (the electric field lies in the plane of incidence) in the form of evanescent waves for $z\rightarrow\pm\infty$. In each piecewise homogeneous region we have that $D_{\mathbf{q}}\cdot\mathbf{E}_{\mathbf{q},\lambda}(z)=0$. The mode-function must then take the form

[TABLE]

where

[TABLE]

with $c_{n}=c/\sqrt{\epsilon_{n}}$ the speed of light in medium $n$, and we introduced the vectors

[TABLE]

Imposing the boundary conditions Eqs. (51) and (52), we obtain the following implicit relation for the surface plasmon dispersion:

[TABLE]

In general, Eq. (58) has no analytical solution, except in the simple case where $\epsilon_{1}=\epsilon_{2}$, in which case its solution reduces to solving a quadratic equation.

The corresponding mode-function can be written as

[TABLE]

and the mode-length is obtained according to Eq. (53) as

[TABLE]

where the last term is due to the dispersion in the graphene layer. We point out, that within the hydrodynamic model used for conductivity of graphene, this contribution is only non-zero if non-local effects are also included, that is if $\beta\neq 0$.

IV.2 Graphene-metal structure

We now move to the more complex case of a graphene sheet near a perfect metal (see right panel of Fig. 1). We assume that the metal interface is located at $z=-d$ and the graphene layer is located at $z=0$. The dieletric constant is given by

[TABLE]

Noticing that the plasmon field should decay for $z\rightarrow\infty$, the mode-function should have the form

[TABLE]

Notice that the presence of a perfect metal at $z=-d$ implies that the tangential component of the electric field should vanish. Imposing the previous condition and the boundary conditions Eqs. (51) and (52) at $z=0$, we obtain a homogeneous system of equations

[TABLE]

where $\xi_{2}=\frac{\epsilon_{2}}{\kappa_{2}}+i\frac{\sigma_{g,L}}{\epsilon_{0}\omega}$. The dispersion relation of the screened plasmons is obtained by looking for the zero the determinant of the previous matrix, which leads to the condition for the dispersion relation

[TABLE]

The boundary conditions imply that the mode-function is given by

[TABLE]

The mode-length can be determined from Eq. 53, and we obtain

[TABLE]

where, as in Eq. (60), the last term is due to the dispersion of graphene.

We note that the dispersion relation for the screened plasmons, Eq. (64), is the same one that would be obtained for the acoustic plasmons in a symmetric graphene double layer structure (center panel of Fig. 1). In this structure, we have two graphene layers located at $z=0$ and $z=-d_{\text{gdl}}$. The dieletric constant of the encapsulating medium is given by

[TABLE]

Since the plasmonic modes should decay for $z\rightarrow\pm\infty$, the plasmon mode-function must have the form

[TABLE]

Imposing the boundary conditions Eqs. (51) and (52) at $z=0$ and $z=d$, we obtain the following homogeneous system of equations

[TABLE]

where $\xi_{2}=\frac{\epsilon_{2}}{\kappa_{2}}+i\frac{\sigma_{g\text{top},L}}{\epsilon_{0}\omega}$ and $\xi_{3}=\frac{\epsilon_{3}}{\kappa_{3}}+i\frac{\sigma_{g\text{bot},L}}{\epsilon_{0}\omega}$. The dispersion relation is obtained from zeroing the determinant of the previous matrix. Since the system is composed of two graphene sheets the double layer structure has two dispersion branches, a low energy one – the acoustic mode – and a high energy one – the optical mode. In the optical mode, the charge oscillations in the two graphene sheets are in-phase; whereas in the acoustic mode, the charge oscillations are out-of-phase. In the particular case where, we have a symmetry structure, with $\epsilon_{2}=\epsilon_{3}$ and the conductivities of the top and bottom graphene layers are the same $\sigma_{g\text{top},L}=\sigma_{g\text{bot},L}=\sigma_{g,L}$, the zeroing of the determinant factorizes into two independent expressions

[TABLE]

for the optical mode, and

[TABLE]

for the acoustic one. Notice that the relation for the acoustic mode dispersion Eq. (70) coincides with the equation for the screened plasmon Eq. (64) provided $d_{\text{dlg}}=2d$. This fact can be understood in terms of image charges as depicted in Fig. 2.

In the bottom panel of Fig. 3 we depict the loss function for the graphene-metal system. Clearly, only one branch is seen, which coincides with the acoustic branch of the double layer graphene upon considering $d$ equal to half that of the double layer system, as explained in Fig. 2.

An alternative way to obtain the plasmon dispersion relation is to look for poles (or resonances in the presence of losses) in the so called loss function(Gonçalves and Peres, 2016) , which is defined as

[TABLE]

where $r_{p}(\mathbf{q},\omega)$ is the reflection coefficient of the structure in consideration for the $p$-polarization, and $\mathbf{q}$ and $\omega$ are the in-plane wavevector and frequency of the impingent radiation, respectively (see Appendix A). For a symmetric graphene double layer ($\epsilon_{1}=\epsilon_{3}$ and $\sigma_{g\text{top},L}=\sigma_{g\text{bot},L}$) and neglecting losses $\gamma\rightarrow 0$, this coefficient has poles at the solutions of Eqs. (70) and (69), as can be seen comparing those equations with Eq. (87). The loss function for the double layer graphene is depicted in the left panel of Fig. 3, as function of $\omega$ and of a dimensionless parameter $s=qc/\omega$ which defines the dispersion relation of the single graphene layer, clad by two different dielectrics of dielectric functions $\epsilon_{1}$ and $\epsilon_{2}$, in the electrostatic limit by

[TABLE]

where $E_{F}$ is the Fermi energy of graphene and $\alpha$ is the fine structure constant of vacuum. In the top panel of Fig. 3 two branches are clearly seen: a high energy one – the optical branch – and the acoustic branch at lower energies. At high energies and high $s$ the two branches merge and converge to the single layer branch. The reader may wonder why the lower branch starts at finite momentum. This happens due to the definition of the $s$ parameter, which involves both the frequency and the real wavenumber $q$. This choice allows to clearly separate the two branches in the $(\omega,s)$ plane.

V Application: quantum emission close to graphene structures

We will now apply the quantization of the plasmon modes in grahene structures to the problem of spontaneous emission by a quantum emitter which is located above the structure. We model the quantum emitter as a two-level system embedded in medium $2$ at position $\mathbf{r}_{0}=\left(0,0,z_{0}\right)$. The quantum emitter couples to the plasmonic field via dipolar coupling:$\hat{H}_{\text{sp-d}}=-\hat{\mathbf{d}}\cdot\hat{\mathbf{E}}_{\text{\text{sp}}}(\mathbf{r}_{0}),$ with

[TABLE]

where $\hat{c}_{g/e}^{\dagger}$($\hat{c}_{g/e}$) is the creation (annihilation) operator for the ground/excited state of the two-level system and $\mathbf{d}_{ge}=\left\langle g\left|\hat{\mathbf{d}}\right|e\right\rangle$ is the dipole matrix element.

The transition rate of an emitter due to emission of surface plasmons in graphene is given by Fermi’s golden rule (Novotny and Hecht, 2006; Archambault et al., 2010):

[TABLE]

where $\omega_{0}$ is the transition frequency, $\left|g;n_{\mathbf{q},\text{sp}}+1\right\rangle$ represents a final state with one more surface plasmon and the emitter in the ground state and $\left|e;n_{\mathbf{q},\text{sp}}\right\rangle$ represents an initial state with $n_{\mathbf{q},\text{sp}}$ plasmons in graphene and the emitter in the excited state. The transition matrix element for spontaneous plasmon emission, when there are no surface plasmons in the initial state,is given by

[TABLE]

With this result the transition rate reads

[TABLE]

Using the mode-function we can write

[TABLE]

$\psi$ is the angle the dipole makes with the axis perpendicular to graphene ($z$-axis), and $\phi$ is the azimuthal angle. The prefactor $\mathcal{N}_{\mathbf{q}}$ is define as $\mathcal{N}_{\mathbf{q}}=1$ for the single-layer graphene case and as $\mathcal{N}_{\mathbf{q}}=\sinh^{2}\left(\kappa_{1,\mathbf{q}}d\right)$ for the graphene+metal structure. Using the in-plane isotropy of the system, the momentum integration in Eq. (77) can be trivially performed, yielding:

[TABLE]

where $q_{0}$ is the momentum of a surface plasmon, with frequency $\omega_{0}$, i.e., $\omega_{0}=\omega_{\mathbf{q}_{0},\text{sp}}$. So far the expression for $\Gamma_{\text{sp}}$ is general. The differences arise from the particular forms of the dispersion $\omega_{\mathbf{q},\text{sp}}$, the mode length $L_{\mathbf{q},\text{sp}}$ and the prefactor $\mathcal{N}_{\mathbf{q}}$.

V.1 Decay rate for local conductivity

We will now study the plasmon emission rate, when non-local effects are neglected, $\beta=0$ in Eq. (50).

We first focus on the case of a quantum emitter close to a single graphene layer and the same dielectric above and bellow the graphene layer, $\epsilon_{1}=\epsilon_{2}=\epsilon$. Using the analytic solution of Eq. (58) for this case, we can write Eq. (79) as

[TABLE]

In Figs. 4 and5, we plot the ratio $\Gamma_{\text{sp}}^{\text{gsl}}/\Gamma_{0}$, where $\Gamma_{0}=d_{ge}^{2}\omega_{0}^{3}/\left(3\pi\epsilon_{0}\hbar c^{3}\right)$ the total decay rate of an emitter in vacuum, for different graphene-quantum emitter distances $z_{0}$ and dipole orientations. For comparison, we will display the ratio $\Gamma_{\text{full}}/\Gamma_{0}$, where $\Gamma_{\text{full}}$ is the total decay rate (including both plasmon and photonic losses) of the quantum emitter. The total decay rate can be computed from the knowledge of the reflection coefficients, which are incorporated into the dyadic EM Green’s function (Agarwal, 1975b; Novotny and Hecht, 2006; Gonçalves and Peres, 2016) (see Ref. (Koppens et al., 2011) for a detailed study of the properties of an emitter near graphene using Dyadic Green’s functions). For a dipole in medium $n=2$, the total decay rate is given by

[TABLE]

where $r_{s/p}$ are reflection coefficients for the $s/p$-polarization and $k_{n}=\omega/c_{n}$ and $k_{z,n}$ are defined in Appendix A.

A number of details are worth mentioning. There are two distinct behaviors of the decay rate $\Gamma_{\text{full}}$. At low frequency the curves shoot up due to Ohmic losses, which are not included in Eq. (80). At intermediate frequencies, the curves develop a clear resonance due to the excitation of surface plasmons in graphene. Also the maximum of the resonances blue-shifts with the decrease of the distance of the dipole to the graphene sheet. This behavior is easily explained remembering that the dispersion of the surface plasmons in single layer graphene is proportional to $\sqrt{q}$. Since the distance $z_{0}$ introduces a momentum scale $q\sim 1/z_{0}$, smaller $z_{0}$ values correspond to higher $q$-values and higher energies of the resonant maximum. Equation (80) for the plasmon emission rate produces exactly the same resonance (same magnitude and same maximum of the resonance position) as $\Gamma_{\text{full}}$, indicating that in this region, the decay rate of the quantum emitter is dominated by plasmon emission. This is further shown in Fig. 5, where losses were arbitrarily reduced in the evaluation of $\Gamma_{\text{full}}$. We can also avoid the superposition of the Ohmic and surface plasmon contributions choosing either a larger Fermi energy or a larger distance from graphene to the metal. In Fig. 5 the ratio $\Gamma/\Gamma_{0}$ is smaller than the ones in Fig. 4 due to the larger distance of the dipole to the graphene sheet.

An analytic expression for the plasmon emission rate can be also obtained for the case of graphene-metal structure assuming that $\kappa_{1}d\ll 1$. In this limit, Eq. 64 can be approximately solved, yielding $\hbar\omega_{\mathbf{q},\text{sp}}=\sqrt{4\alpha dE_{F}\hbar c/\epsilon}q$. Plugging this result in Eq. 79, we obtain

[TABLE]

This plasmon emission rate is shown in Fig. 6 as a function of $\omega$ and the graphene-metal separation $d$. As in the case of a single graphene layer, for each $d$ there is a well defined peak as function of $\omega$.

V.2 Decay rate with non-local effects

We now focus on the role of non-local effects in the graphene conductivity play in the emitter decay rate. We will restrict ourselves to the case of a quantum emitter close to a single layer of graphene. In Fig. 7, we compare the transition rate of a quantum emitter calculated taking into account non-local effects $\beta\neq 0$ in Eq. (37), with the local case $\beta=0$. The left panel shows the result for $E_{F}=0.6$ eV, $z_{0}=7$ nm, while the right panel shows the result for$E_{F}=0.6$ eV, $z_{0}=2$ nm. We can see that the nonlocal effects play an important role at smaller emitter-graphene distances. As the distance between quantum emitter and the graphene layer, $z_{0}$, decreases, non-locality leads to an increase and blueshift of the resonance in the transition rate as a function of the emitter frequency. We have also verified (not shown) that the non-local effects are also more prominent in the case of lower electronic densities.

VI Conclusions

In this paper, we have performed the quantization of graphene plasmons, in the absence of losses, and applied the field quantization to the interaction of an emitter with doped graphene. The quantization was performed using both a macroscopic energy approach and a quantum hydrodynamic model, which allows for the inclusion of non-local effects in the EM response of graphene. The quantization approaches used allow for the determination of the plasmon EM field mode-functions and, importantly, their normalization, which becomes non-trivial when dispersion is included.

When comparing the total decay rate of a quantum emitter (as obtained using the full EM dyadic Green’s function) with the decay rate due to quantum emission, it was shown that plasmon emission completely dominates the decay rate, for typical emitter-graphene separations and emitter frequencies. It was shown that non-local effects in the graphene response, become increasingly important for smaller graphene-emitter separations and smaller Fermi energies.

The advantage of the quantization method developed in this work lies in its simplicity. The mode-functions are obtained from the solution of the Maxwell equations for the vector potential, and the normalization of the mode-functions can be expressed as a simply integral, which only involves the dielectric function of the medium. For situations when only a few modes contribute significantly to the physics, as in the case of quantum emission dominated by plamons, mode-functions allows to a much simpler and physically transparent description than the full EM dyadic Green’s function. Since, the determination of the mode-functions only involves the solution of the classical Maxwell equations, the quantization of the electromagnetic field of more complexes structures can be easily achieved. We also note that the procedure gives both the quantized form of the electric and magnetic fields. Therefore, we could also study the enhancement of the spin relaxation. The only change would be the interaction Hamiltonian, which in this case would be of the form $H_{I}=-\bm{\mu}\cdot\mathbf{B}$, where $\bm{\mu}$ is the magnetic dipole moment of the emitter and $\mathbf{B}$ the quantized magnetic field of the surface plasmon.

In possession of a quantized theory for graphene plasmons, we have set the stage for the future discussion of other quantum effects involving these collective excitation made simultaneously of light and matter.

Acknowledgments

B.A.F. and N.M.R.P. acknowledge support from the European Commission through the project “Graphene-Driven Revolutions in ICT and Beyond” (Ref. No. 785219), and the Portuguese Foundation for Science and Technology (FCT) in the framework of the Strategic Financing UID/FIS/04650/2013. Additionally, N. M. R. P. acknowledges COMPETE2020, PORTUGAL2020, FEDER and the Portuguese Foundation for Science and Technology (FCT) through projects PTDC/FIS-NAN/3668/2013 and POCI-01-0145-FEDER-028114. B. A. acknowledges the hospitality of CeFEMA where he was a visiting researcher during part of the time in which this work was developed.

Appendix A Reflection coefficients of graphene structures

In this appendix we provide the reflection coefficients needed for the evaluation of the full decay rate of a quantum emitter, Eq. (81).

A.1 Reflection coefficients for a single graphene layer

The Fresnel problem for a single graphene sheet has been considered in Ref. (Gonçalves and Peres, 2016). Therefore we provide here the final results only and for simplicity we assume we are dealing with non-magnetic media. For an incoming wave from region $2$, and for the $s-$polarization the reflection and transmission coefficient read

[TABLE]

whereas for the $p-$polarization we have

[TABLE]

In these equations, we have $k_{z,n}=\sqrt{\epsilon_{n}\omega^{2}/c^{2}-q^{2}}=i\kappa_{n}$, and $\sigma_{g,L/T}$ is the optical longitudinal/transverse conductivity of graphene, $\epsilon_{0}$ and $\mu_{0}$ are the vacuum permitivity and permeability, respectively, $q$ is the wavenumber along the graphene sheet, and $\omega$ is the frequency of the electromagnetic radiation.

A.2 Metal-graphene reflectance coefficient

The reflection coefficient for the $p$-polarization for this structure is given by

[TABLE]

where $k_{n}^{2}=\epsilon_{n}\omega^{2}/c^{2}$, $d$ is the graphene-metal distance, and $\kappa_{n}^{2}=k_{n}^{2}-q^{2}$, where $q$ is the wavenumber along the graphene sheet. For the $s-$polarization we have

[TABLE]

Appendix B Diagonalization of Hydrodynamic Hamiltonian

B.1 Orthogonality of mode-functions

We start by showing that the mode-functions $\left(\mathbf{A}_{\mathbf{q},\lambda}(z),\,\bm{\upsilon}_{\mathbf{q},\lambda}\right)$ obey certain orthogonality conditions which will be useful. We start by writing the equations for the mode-functions within the hydrodynamic model, Eqs. (32) and (33), in matrix form as

[TABLE]

Let us now consider another mode $\lambda^{\prime}$, with mode-function $\left(\mathbf{A}_{\mathbf{q},\lambda^{\prime}}(z),\,\bm{\upsilon}_{\mathbf{q},\lambda^{\prime}}\right)$, solution of

[TABLE]

Now we contract Eq. (89) with \left[\begin{array}[]{cc}\mathbf{A}_{\mathbf{q},\lambda^{\prime}}^{\dagger}(z)&\bm{\upsilon}_{\mathbf{q},\lambda^{\prime}}^{\dagger}\end{array}\right], Eq. (89) with \left[\begin{array}[]{cc}\mathbf{A}_{\mathbf{q},\lambda}^{\dagger}(z)&\bm{\upsilon}_{\mathbf{q},\lambda}^{\dagger}\end{array}\right] and integrate both equations over $z$ obtaining

[TABLE]

Taking the conjugate of Eq. (102), we obtain

[TABLE]

Subtracting this last equation from Eq. (96), we obtain

[TABLE]

If $\omega_{\mathbf{q},\lambda}\neq\omega_{\mathbf{q},\lambda^{\prime}}$, we can divide by $\omega_{\mathbf{q},\lambda}-\omega_{\mathbf{q},\lambda^{\prime}}$ obtaining one of the orthogonality conditions:

[TABLE]

We can obtain an additional orthogonality relation. If we take the complex conjugate of Eq. (96), replace $\mathbf{q}\rightarrow-\mathbf{q}$ and $\lambda\rightarrow\lambda^{\prime}$, and using the fact that $\omega_{\mathbf{q},\lambda}=\omega_{-\mathbf{q},\lambda}$, we obtain

[TABLE]

We now contract this equation with \left[\begin{array}[]{cc}\mathbf{A}_{\mathbf{q},\lambda}^{\dagger}(z)&\bm{\upsilon}_{\mathbf{q},\lambda}^{\dagger}\end{array}\right], integrate over $z$ and take its complex conjugate, obtaining

[TABLE]

Contracting Eq. (89) with \left[\begin{array}[]{cc}\mathbf{A}_{-\mathbf{q},\lambda}^{t}(z)&\bm{\upsilon}_{-\mathbf{q},\lambda}^{t}\end{array}\right] and integrating over $z$ we obtain

[TABLE]

Subtracting Eq. (107) from (108), we obtain

[TABLE]

For $\omega_{\mathbf{q},\lambda}+\omega_{\mathbf{q},\lambda^{\prime}}\neq 0$, we obtain a second orthogonality condition:

[TABLE]

B.2 Hamiltonian in term of mode amplitudes

We now insert the expansion of the fields in terms of mode-functions, Eqs. (30) and (31), into the classical Hamiltonian (29). Using Eqs. (27) and (28), we can write

[TABLE]

where

[TABLE]

Using the mode-function equation (89), we can write

[TABLE]

Using the orthogonality condition (110) we conclude that $\tilde{h}_{\mathbf{\lambda,\lambda^{\prime}}}(\mathbf{q})=0$. The orthogonality condition (105) implies that $h_{\mathbf{\lambda,\lambda^{\prime}}}(\mathbf{q})=0$, except when $\lambda=\lambda^{\prime}$ (assuming there are no degeneracies in $\omega_{\mathbf{q},\lambda}$). Therefore, we see that the classical Hamiltonian can be written as in Eq. (40), with the mode-length being biven by

[TABLE]

which can be written as Eq. (41).