Plasmons in layered structures including graphene

We investigate the optical properties of layered structures with graphene at the interface for arbitrary linear polarization at finite temperature including full retardation by working in the Weyl gauge. As a special case, we obtain the full response and the related dielectric function of a layered structure with two interfaces. We apply our results to discuss the longitudinal plasmon spectrum of several single and double layer devices such as systems with finite and zero electronic densities. We further show that a nonhomogeneous dielectric background can shift the relative weight of the in-phase and out-of-phase mode and discuss how the plasmonic mode of the upper layer can be tuned into an acoustic mode with specific sound velocity.


Introduction
Graphene, the two-dimensional allotrope of carbon, has become one of the most active fields in today's both experimental and theoretical condensed matter physics [1,2,3,4]. Among many extraordinary phenomena and proposals, the optical properties of graphene have generated particular interest due to potential applications, [5,7,8,6] but also because they are intimately related to the discovery of exfoliated graphene. [9] Recently, plasmonics based on graphene has become a new emerging subfield, trying to take advantage of the strong electronic confinement and long propagation lengths of its carriers and, most importantly, the possibility of applying an electrostatic gate voltage. [10,11,12,13,14,15] Due to the linear spectrum, different gate voltages can have strong effects on the carrier concentration and thus on the plasmonic spectrum, especially for small nano-islands. [17,16] There is also renewed focus on layered structures due to the experimental advances of exfoliating a number of materials and placing them on top of each other. [18,21,19] Like this, the Coulomb drag of closely separated graphene layers was observed, [20] which has triggered considerable attention by various theoretical groups. [22,23,24,25,26,27,28,29,30] In this paper, we want to combine these two fields and analyze the plasmonic spectrum of layered structures at finite temperature.
To discuss the undamped plasmon dispersion, it suffices to determine the zeros of the dielectric functions [31] or, alternatively, the zeros of the denominator of the transmission or reflection amplitude. [32] But for a finite imaginary part, i.e., especially for finite temperatures, this procedure is ambiguous. We will thus need a different approach and will investigate the energy loss function. For single layer graphene systems, the energy loss function is related to the (negative) imaginary part of the inverse of the dielectric function, but for double layer structures, this function changes its sign. [32] We will thus define the energy loss function as the trace of the (negative) imaginary part of the full response.
Another common approximation to obtain the longitudinal plasmon dispersion consists in using the static Coulomb potential and for small wave numbers q k F in the regime of small layer separation k F d 1, the Coulomb potential can be further simplified, only depending on the average of the outer dielectric media ( 1 + 3 )/2. For the general case, the full electrostatic problem has to be considered, which has been recently done in the context of graphene double layers. [32,33,34] Since for small energies of the order of α F (α being the fine-structure constant) retardation effects lead to strong light-matter interaction, [35] we will here derive the full retarded photon propagator for the longitudinal and transverse channel. Another focus of this work is placed on materials with a large dielectric constant which strongly screen the graphene layers. SrTiO 3 , e.g., has a relative dielectric constant of ∼ 300, which can reach up to ∼ 5000 at liquid helium temperature due to the proximity of a ferroelectric instability. [36] Another example are surface states of a three-dimensional topological insulator which are separated by the width of the sample. For Bi 2 Te 3 , the two electronic Dirac systems are then electrostatically coupled through a dielectric medium with ∼ 100. [33] Particularly, we will show that strong dielectrics shift the relative weight of the in-phase and out-of-phase mode. Used as a substrate, they strongly screen the graphene layers and therefore basically act like metals, leading to a linear plasmon dispersion. [37,38] The paper is organized as follows. In section II, we derive the photon propagator in free space, separating the final result into longitudinal and transverse channels. In section III, we formulate the linear response theory for a layered structure including graphene in the interfaces and give explicit expressions for the double layer geometry. In section IV, we finally discuss the near-field optical properties of layered graphene structures, contrasting between single and double layer, high-and low temperature, large and small dielectric constants. We close with a summary and an outlook. An appendix outlines the analytical discussion how to obtain the linear plasmon dispersion for a double layer system with large-substrate.

Photon propagator in a homogeneous medium
In this section and the following section, we will derive the retarded photon propagator in a homogeneous and nonhomogeneous medium, outlining all details. For alternative introductions to confined photon systems, we refer the reader to Ref. [39]. Even though our treatment will be entirely classical based on Maxwell's equations, we will still call the final result, Eq. (15), the photon propagator since this expression coincides with the quantum mechanical photon propagator. [40] The reader only interested in the results presented in Sec. IV may skip this and the following section.
We depart from Ampere's circuital law equation including Maxwell's displacement current assuming a harmonic time evolution with frequency ω: We will first recall the usual representation of the photon propagator in Cartesian coordinates and then introduce the representation in cylindrical coordinates, suitable to discuss layered structures.

Photon propagator in Cartesian coordinates
Introducing the electrostatic potential φ and the vector potential A by we obtain from Eq. (1) the following equation:

Lorentz gauge
In the Lorentz-gauge ∇ · A = iωµµ 0 ε 0 φ, Eq. (4) can be written as four independent inhomogeneous Helmholtz equations with the four-dimensional vector potential A µ = (φ, A) and the four-dimensional current j µ = (ρ/ ε 0 , µµ 0 j). The dispersion relation reads where c denotes the speed of light in vacuum and c 1 = c/ √ µ the (slower) speed of light in the dielectric medium characterized by µ and . The Green's function defined by is thus a scalar and reads where the plus-sign defines propagation out of the source and the minus-sign convergence into the source. Within the Lorentz gauge, the electric field due to a given current density j, as defined in Eq. (2), is thus given by

Weyl gauge
In the following, we will not work in the Lorentz gauge, but will set the electrostatic potential equal to zero, i.e., φ = 0. This gauge condition is often refered to as the "Weyl gauge". The Weyl gauge implies that we only need to consider the propagation of the vector potential and is often used when interactions with non-relativistic particles are involved. Graphene's linear response to the incoming light field is thus entirely defined by its current density. With φ = 0, Eq. (4) becomes In this gauge, the operator acting on A is a vector that connects different spacial directions. For a general solution, we thus need to determine the dyadic Green's function defined by where I is the 3 × 3 unit tensor. At first glance, it seems that the Weyl gauge results in more complicated expressions of the Green's function. This is true in real space, but in Fourier space, the expressions are only slightly more involved. And by having eliminated the electrostatic potential φ, the boundary conditions of the subsequent scattering problem become more compact since they only have to be satisfied by the vector potential A. We can easily obtain the dyadic Green's function by noting that the Maxwell's equations yield the following expression for the electric field: The electric field is thus defined by the same dyadic Green's function as the vector potential in the Weyl gauge. Since the electric field is gauge independent, we can use the expression of Eq. (9) to deduce the dyadic Green's function G from the scalar Green's function G 0 .
In real space, we obtain and in Fourier space With G 0 (k, ω) = (k 2 − k 2 0 ) −1 , we finally obtain the retarded photon propagator D αβ Decomposing it in longitudinal and transverse components, it reads with the functions D L,T (k, ω) given by

Photon propagator in cylindrical coordinates
In a layered structure, assumed to be perpendicular to the z axis, the components parallel to the interface, q = (q x , q y ), will be preserved as a good quantum number. It is, therefore, convenient to employ the following representation for the Green's function in a homogeneous medium: with k = (q, k z ). In this representation, the tensor components have a different structure depending on whether α, β = i, j with i, j = x, y or α, β = z. For the in-plane components of the tensor D α,β 0 , we have with q = q 2 − (ω/c 1 ) 2 . Decomposed into longitudinal and transverse contributions, we obtain with the in-plane propagators d 0 l,t (q, ω) given by For the cross terms of D α,β 0 , we get This shows that in-plane longitudinal sources can generate not only in-plane, but also outof-plane fields with a phase shift of π/2. On the other hand, in-plane transverse currents can only generate in-plane transverse fields. Finally, out-of-plane sources generate out-of-plane fields given by the tensor component 3 Linear response of a layered geometry including graphene We now consider the experimentally important situation of layered structures, i.e., we assume the existence of well-defined interfaces at which the material properties are discontinuous. The Maxwell equations written in integral form then yield boundary conditions for the normal and tangential field components, see e.g. Ref. [41]. For the normal components they read with ρ the charge density on the interface. For the tangential components they read with j the current density on the interface. An arbitrarily polarized electromagnetic wave can always be expressed by superposing two linearly polarized waves which are orthogonal to each other. We can thus define p-polarized and s-polarized waves, respectively, according to the plane of incidence. Alternatively, we will use the denomination of longitudinal and transverse polarization.
The boundary conditions for the normal and tangential field components are not independent of each other since they are connected by Maxwell's equations. According to the plane of incidence, we will either use Eqs. (24) in the case of longitudinal polarization or Eqs. (25) in the case of transverse polarization. These conditions then simplify considerably when working in the Weyl gauge.
Notice that the influence and properties of graphene are only accounted for via the charge and current densities on the interface, ρ and j. The properties of graphene thus enter when matching the vector field at the interface of two adjacent dielectric media. Since we have set φ = 0, these properties are entirely contained in the current-current response of graphene, χ 0 ij . The superindex 0 denotes the bare current response, i.e., the response to the total (external plus induced) vector potential.
We will now assume an isotropic system such that current-current response tensor can be split up into a longitudinal and a transverse contribution, i.e., Within the Dirac approximation this decomposition is always possible and only for transitions close to the van Hove singularity the full tensor structure needs to be considered. [45] The full photon propagator in the presence of a single graphene layer at the location z 1 modifies the "vacuum" propagator in the following way: where summation over repeated indices is implied. χ ij represents the total current-current response of graphene to external fields which, decomposed into longitudinal and transverse contributions, is given by Longitudinal and transverse components thus decouple and in the following, we will only explicitly label these different channels when it is necessary. The full photon propagator of a general layered structure is obtained by superposing all possible scattering paths similar to Eq. (27). Let us finally recall that the gauge field at any point r due to a current source at r is given through the photon propagator via Due to this linear correspondence, the photon propagator can be deduced from the scattering problem of the gauge field using the expressions of D αβ 0 in a homogeneous medium. In the following will outline the general scattering problem for one interface, discussing both, the longitudinal and transverse channel. We will then give the explicit expressions for the retarded photon propagator, dielectric function and graphene loss function for a arbitrary double layer structure.

Scattering on one interface
All general properties of the photon propagator of a layered geometry can be deduced from the scattering problem of one single interface. We will, therefore, discuss this simple problem in some detail. The generalizations are then straightforward. In the following, we will set the plane of incidence the xz-plane and the interface at z = 0.

Longitudinal polarization
For longitudinal polarization, the general vector field has a component parallel (x) and normal (z) to the interface: With q i = q 2 − (ω/c i ) 2 and c i the speed of light in the corresponding medium, we make the ansatz (j = , ⊥) The components of A ⊥ are obtained from the components of A via the condition for a transverse field ∇ · A = 0. This gives the following relations: From Eqs. (24), we see that the parallel component of the vector field is continuous at the interface, but the normal component of the displacement field makes a jump if a graphene layer is present leading to a charge density ρ at the interface. This component is related to the vector field via the relation D ⊥ = ε 0 iωA ⊥ . With the continuity equation ωρ − q · j = 0 and the linear response j = −χ 0 l A q , the set of equations closes. Together with Eq. (32), we thus obtain the following two conditions: The transmission and reflection amplitude for longitudinal polarization then read

Transverse polarization
For transverse polarized light and the plane of incidence again in the xz-plane, only the y-component A y of the vector field is non-zero. We thus make the following ansatz: From Eqs. (25) we see that the vector potential is continuous at the interface and that the first derivative makes a jump due to the current generated by the vector field inside the graphene plane. The current is again related to the corresponding transverse currentcurrent susceptibility, χ 0 t , via linear response. [44,45] We thus obtain the following two conditions: The transmission and reflection amplitude for transverse polarization then read

Scattering on two interfaces
We will now discuss in some detail the special case of a double layer graphene structure by applying the basic matching conditions discussed in the previous subsection to two interfaces. We will give explicit expressions for the retarded photon propagator, the dielectric function and the full response of graphene including the definition of the energy loss function.

Photon propagator
For two and more interfaces it is convenient to write Fourier transformed part of the photon propagator as a n × n-matrix, n denoting the number of interfaces. The full photon propagator then satisfies the usual Dyson-like equations (one for each polarization channel) where for the special case of two graphene layers the matrix χ 0 = diag(χ 0 1 , χ 0 2 ) represents the bare graphene's response in layer 1 (χ 0 1 ) and layer 2 (χ 0 2 ). d 0 is the photon propagator in the absence of graphene (χ 0 i = 0), but with the dielectric geometry of Fig. 1. The entries of the matrix d can be obtained from the standard matching conditions or, equivalently, using multiple scattering formalism. In the latter case, they can be written as where d i = with the obvious expression for index exchange, where r ij (t ij ) are the corresponding coefficients for a single graphene interface from medium i into medium j, previously obtained. Notice that a recursive interpretation of the last formula can be used for calculating the reflection and transmission amplitudes for any layered structure with multiple interfaces. Solving the matching conditions for the scattering problem with two interfaces but without graphene (χ 0 = 0) directly, one can obtain more explicit expressions for the photon propagator in the absence of graphene, d 0 . Using the results of Ref. [32], one obtains for longitudinal polarization the following compact result: with N l,0 = q 2 2 (q 3 1 + q 1 3 ) cosh(q 2 d) + (q 2 2 1 3 + q 1 q 3 2 2 ) sinh(q 2 d). For the transverse part, one obtains the corresponding expression with N t,0 = µ 2 q 2 (µ 3 q 1 + µ 1 q 3 ) cosh(q 2 d) + (µ 2 2 q 1 q 3 + µ 1 µ 3 q 2 2 ) sinh(q 2 d).

Graphene's response
Graphene's response obeys similar equations (one for each polarization channel) which, together with Eq. (42), provide the complete dynamics of the coupled matter-field system. The retarded dielectric function of double layer graphene with nonhomogeneous background is then often defined by [31] (q, Usually −Im −1 is used to discuss the plasmonic spectrum, but this function changes sign and can thus not be interpreted as (positive definite) spectral density. But instead of the determinant, we find it more convenient to discuss the trace of the the full response matrix. Graphene's excitations correspond to the imaginary part of the full response, and to reveal its presence we will discuss the following generalization of the energy loss function S(q, ω) to several layers: Since S(q, ω) is related to the imaginary part of a causal function, it is strictly positive and since it also proportional to the usual definition of the energy loss function for a single layer, Eq. (50) can serve as a straightforward generalization of the energy loss function for arbitrary layered systems.
Let us briefly comment on the physical interpretation of the response matrix and the related energy loss function. The diagonal entries of the response matrix χ are given by the response of a particular layer if the gauge field is only applied to just this particular layer. Diagonalizing the response matrix χ, one can discuss the elementary excitations of the full system, separately. This was done in Ref. [35], where for double-layer graphene the in-phase and out-of-phase excitations were analyzed. Since the trace of the response matrix χ is invariant with respect to unitary transformations, it serves as natural choice for the definition of the generalized energy loss function and for the discussion of the internal excitations of the whole system, i.e., the sum of excitations of all layers. As mentioned above, the imaginary part of the general response is related to graphene's excitations. In the following section, we will use the developed formalism to discuss the spectrum of longitudinal plasmonic excitations for single and double layer structures. For this, retardation can formally be neglected. But we stress that the formalism can also straightforwardly be used for multiple layer structures as well as to discuss the spectrum of transverse plasmonic excitations where retardation effects are crucial.

Plasmons in layered structures at finite temperature
In this section, we will apply our formalism and discuss the longitudinal response of layered structures at finite temperature. To do so, we will use the density-density correlation function P 0 (q, ω) which is related to the longitudinal component of the current-current correlation function χ 0 l (q, ω) via the continuity equation. [45] Within the Dirac-cone approximation, this reads

Polarizability of graphene
Up to now, there is no approximation involved except of decomposing the current response into a longitudinal and transverse contribution which is well justified for any transitions not too close to the van Hove singularity. We will now approximate the full density-density correlation function by the non-interacting polarizability [42]  with E ± (k) = ± v F k the eigenenergies, n F (E) = (e β(E−µ) + 1) −1 the Fermi function, and g s = g v = 2 the spin and valley degeneracy for graphene. A characteristic difference between the polarizability of graphene and that of a two-dimensional electron gas is the appearance of the prefactors f ss (k, q) coming from the band-overlap of the wave function where ϕ denotes the angle between k and q. At zero temperature, we have µ = F with F the Fermi energy and the analytic solution of Ref. [42] can be decomposed in several patches corresponding to inter-and intraband transitions, respectively. In Figs. 2-6, the three most basic patches are separated by black dashed lines.
At finite temperature the chemical potential is determined with respect to the electronic density n by the following relation: with the density-of-states given by ν( ) = g s g v | |/(2πv 2 F ). For an electron-doped system, we have F > µ > 0. At the neutrality point, n = 0 and due to particle-hole symmetry we then have µ = 0. For our numerical calculations, we will make use of the semi-analytical expression of the polarizability presented in Ref. [43].

Loss function at finite and zero doping
We will now consider single and double layer structures at finite and zero doping at low and high temperatures. To clarify the discussion we will choose a homogeneous medium, i.e., 1 = 2 = 3 = 1. Figure 4: (color online): The energy loss function −Imχ(q, ω + i0) in units of T / 2 ( T = v F k T ) for longitudinal polarization for single layer (left) and double layer with d = 2nm (right) at a temperature T = 300K at zero doping n = 0 with the same dielectric medium for all regions = 1 (air). Also shown the zero-temperature plasmon dispersion for finite doping with k T = k B T v F 2 ln 2 (black solid lines).
We first discuss the plasmon dispersion at finite temperature for suspended single and double layer graphene with a layer separation of k F d = 0.35. For an electronic density of n = 10 12 cm −2 , we then have d = 2nm and T F = 1200K, such that T = T F /4 would correspond to approximately room temperature. The energy loss function is shown for T = T F /4 in Fig. 2 and for T = T F for 3.
We compare the energy loss function with the plasmon dispersion at zero temperature by determining det = 0. In the case of a finite imaginary part of P 0 (q, ω), i.e., Landau damping, we set ImP 0 =0. This guarantees the convergence of the two plasmonic modes for large q. [32] We see that at intermediate temperatures the plasmon dispersion is redshifted compared to the zero temperature solution, whereas for high temperatures it is blue-shifted. This follows directly from the behavior of P 0 (q = 0, ω → 0) which is related to the Drude weight D and defines the plasmon dispersion. [42] At finite temperature, this reads which is a non-monotonic function of the temperature due to the temperature dependence of the chemical potential given in Eq. (54). In Fig. 4, we show the plasmon dispersion at the neutrality point n = 0 at T = 300K. It was already discussed earlier that also in this limit (damped) plasmon excitations exist in single layer graphene due to the weak Landau damping. [46] In fact, there is an analytical approximation which uses the zero temperature result of the polarizability [42] with the replacement k F → k T = k B T v F 2 ln 2. [47] This agrees with the physical expectation that at finite temperature there is a finite electronic density due to thermal broadening of the Fermi function.
The parameter in Fig. 4 corresponds to a thermal wave number k T = 0.061nm −1 and electron density n = 1.2 × 10 11 cm −2 . The energy loss function is given in units of T / 2 with T = v F k T ∼ 36meV. As one can see, the analytical approximation agrees well with the maximum of −Imχ(q, ω + i0) even in the case of double layer graphene.

Loss function for nonhomogeneous dielectric media
We will now discuss the effect of a nonhomogeneous dielectric media which can lead to changes of the plasmon dispersion due to the different photon propagator. In Fig. 5 we show the results for a double layer graphene structure for two temperatures T = T F /10 (left) and T = T F (right). We choose SiO 2 with 3 = 3.8 as a substrate and Al 2 O 3 with 2 = 6 as a buffer layer. It is further assumed that the top layer is air with 1 = 1. The two layers have the same electron density and layer separation k F d = 1.77 which for n = 10 12 cm −2 corresponds to d = 10nm.
We use the same parameters as in Ref. [34] and the results of this reference are shown as green solid lines, obtained from the real part of det = 0. Since the authors do not set the imaginary part to zero, the in-phase and out-of-phase modes do not merge for q > ∼ k F in contrary to the physical expectation. The correct way to approximate the plasmon dispersion from the zeros of the dielectric function is thus to set the imaginary part of the current-current correlation function to zero (black solid lines). Even after including temperature broadening and Landau damping, the discussion of Ref. [34] remains incorrect since it is based on the response of single layer graphene.
Apart from the mere (numerical) corrections of the plasmonic dispersion due to a nonhomogeneous dielectric background, it is possible to shift spectral weight from the out-of-phase to the in-phase mode. One can also change the plasmon dispersion from the typical square-root dispersion of charged plasmonic waves to a linear dispersion typical for acoustic sound waves. In Figs. 2-4, where the same dielectric constant was chosen to be the one of air, the out-of-phase mode is more dominant extending towards larger wave numbers. In Fig. 5, we see that the two modes are almost equal in weight merging in the region of finite Landau damping. If we choose now a substrate with a very large dielectric constant, the in-phase mode becomes more dominant whereas the out-of-phase mode basically merges with the Dirac cone. Since the latter behavior is also seen in the bare graphene response, we can thus state that there is only an in-phase plasmonic mode. This is shown on the left hand side of Fig. 6.
The plasmonic mode is also dominant in the case of a topological insulator where the buffer layer is resembled by a strong dielectric medium, i.e., we set 1 = 1, 2 = 100, and 3 = 4. This can be seen on the right hand side of Fig. 6 where we use the same parameters as in Ref. [33]. They are given by g v = g s = 1 with a sample width of k F d = 2.37 and a Fermi velocity of v T I = 5 × 10 5 m/s which scales out.
Let us finally comment on the linear plasmon mode present on the left hand side of Fig. 6. In the appendix, we give details of how to derive and approximate the plasmon dispersion of the upper layer as ω p = v a q. For this we do not include the full expression of the current-current correlation function, but only use the long wavelength limit. This yields the simple result for the group velocity of the linear (acoustic) mode v a = 2α g dk 1 where k 1 F is the Fermi wave number of the upper graphene layer and α g = αc/v F graphene's fine-structure constant. This approximation breaks down for small dk F and large 2 since v a cannot become smaller than the Fermi velocity [32] and a more careful analysis is necessary. [33] Eq. (56) can nevertheless be used to discuss several aspects. First, v a only depends on the buffer substrate; second, only the upper graphene layer enters in the expression since the lower one is perfectly screened. This formula can thus not only be applied to double layer graphene on substrates with large dielectric constants, but also to single layer graphene on top of metals interpreting d as the graphene-metal distance. We can thus calculate the group velocity of the plasmonic mode for a single layer graphene on top of Pt(111) as discussed experimentally in Ref. [38] where the energy loss function was measured. With the parameters F = 0.3eV, d = 0.33nm and 2 = 1 corresponding to the experimental set-up, we have k F d = 0.15 and thus v a = 1.15v F . The screened plasmon dispersion will thus lie close to the Dirac "light-cone" in agreement with experiment.

Summary
We presented results of the energy loss function of layered structures which is the only way to unambiguously discuss the plasmon dispersion in the presence of dissipative terms like Landau damping or finite temperature. We discussed the effect of temperature and nonhomogeneous dielectric medium, including large dielectrics which almost perfectly screen the graphene layers, but our formalism equally applies to systems with unbalanced electronic densities. Our main results are (i) the possibility of shifting relative weight of the several plasmonic branches by changing the nonhomogeneous dielectric background and (ii) a simple formula for the sound velocity of the (linear) plasmonic mode of the upper graphene layer in the case of a substrate with large dielectric constant. These insights might be useful towards the engineering of specific plasmon modes for future plasmonic circuitries based on graphene.

Acknowledgments
This work has been supported by Portugal's Fundação para a Ciências e a Tecnologia

Appendix: Large substrate screening and acoustic plasmons
In this appendix, we show that a large (huge) value of 3 renders medium 3 almost a metal, largely screening the Coulomb interaction in the upper graphene layer and turning the otherwise square-root plasmon into an acoustic mode over a wide wave number range. This is similar to the proposal of Ref. [44] where a perfect metal as substrate is considered. The acoustic mode evolves into a regular two-dimensional plasmon only at much reduced wave numbers, due to the incomplete screening of medium 3.

Linear plasmon mode of the upper graphene layer
We assume the usual geometry where the three dielectric media separate the upper and lower graphene layer, see Fig. 1. We assume a very large value of 3 ≈ 300, much greater than the other regular values of 1,2 , and ignore retardation. Plasmons are solutions of (see Eq. (45)) 1 − r 21 r 23 e −2qd = 0.
The reflection amplitude r 23 is given by (see Eq. (36) in the limit c → ∞)