{'fraction_non_alphanumeric': 0.03917353397826347, 'fraction_numerical': 0.01695927385644333, 'mean_word_length': 4.319779754341381, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Double-layer graphene and topological insulator thin-film plasmons

We present numerical and analytical results for the optical and acoustic plasmon collective modes of coupled massless-Dirac two-dimensional electron systems. Our results apply to topological insulator (TI) thin films and to two graphene sheets separated by a thin dielectric barrier layer. We find that because of strong bulk dielectric screening TI acoustic modes are locked to the top of the particle-hole continuum and therefore probably unobservable.


I. INTRODUCTION
The physics of closely-spaced but unhybridized twodimensional electron systems (2DESs) has been a subject of theoretical and experimental interest since it was first appreciated 1,2 that electron-electron interactions allow energy and momentum to be transferred between layers, while maintaining separate particle-number conservation. Remote Coulomb coupling has commanded a great deal of attention during the past thirty years or so because it provides a potential alternative to the inductive and capacitive coupling of conventional electronics. Until recently, remote Coulomb coupling research focused on quasi-2D electron systems confined to nearby quantum wells in molecular-beam-epitaxy grown semiconductor heterostructures. The study of Coulomb-coupled 2D systems has now been revitalized by advances which have made it possible to prepare robust and ambipolar 2DESs, based on graphene 3 layers or on the surface states of topological insulators 4 , that are described by an ultrarelativistic wave equation instead of the non-relativistic Schrödinger equation.
Single-and few-layer graphene systems can be produced by mechanical exfoliation of thin graphite or by thermal decomposition of silicon carbide 5 . Isolated graphene layers host massless-Dirac two-dimensional electron systems (MD2DESs) with a four-fold (spin × valley) flavor degeneracy, whereas topologicallyprotected MD2DESs that have no additional spin or valley flavor labels appear automatically 4,6 at the top and bottom surfaces of a three-dimensional (3D) TI thin film. The protected surface states of 3D TIs are associated with spin-orbit interaction driven bulk band inversions. 3D TIs in a slab geometry offer two surface states that can be far enough apart to make single-electron tunneling negligible, but close enough for Coulomb interactions between surfaces to be important. Unhybridized MD2DES pairs can be realized in graphene by separating two layers by a dielectric 7 (such as Al 2 O 3 ) or by a few layers of a one-atom-thick insulator such as BN 8,9 . In both cases inter-layer hybridization is negligible and the nearby graphene layers are, from the point of view of single-particle physics, isolated. Isolated graphene layers can be also found on the surface of bulk graphite 10,11 and in "folded graphene" 12 (a natural byproduct of micromechanical exfoliation), or prepared by chemical vapor deposition 11 . We use the term double-layer graphene (DLG) to refer to a system with two graphene layers that are coupled only by Coulomb interactions, avoiding the term bilayer graphene which typically refers to two adjacent graphene layers in the crystalline Bernal-stacking configuration 13 .
DLG and TI thin films are both described at low energies by a Hamiltonian with two MD2DES 3 coupled only by Coulomb interactions. The importance of electronelectron interactions in MD2DESs has been becoming more obvious as sample quality has improved 14 , motivating investigations of charge and spin or pseudospin dynamics in DLG and thin-film TIs in the regime in which long-range Coulomb forces give rise to robust plasmon collective modes 15,16 . Because of their electrically tunable collective behaviors, DLG and thin-film TIs may have a large impact on plasmonics, a very active subfield of optoelectronics 17-19 whose aim is to exploit plasmon properties in order to compress infrared electromagnetic waves to the nanometer scale of modern electronic devices.
In this Article we use the random phase approximation (RPA) 15,16 to evaluate the optical and acoustic plasmon mode dispersions in DLG and in thin-film TIs. In particular, we obtain an exact analytical formula for the RPA acoustic plasmon group velocity valid for arbitrary substrate and barrier dielectrics that points to a key difference between these two MD2DES's, namely that the velocity in TI thin films is strongly suppressed. The RPA collective modes of DLG have been calculated earlier by Hwang and Das Sarma 20 : below we will comment at length on the relation between our results and theirs. Plasmon collective modes formed from TI surface states have also been considered previously by Raghu et al. 21 in the regime in which coupling between top and bottom surfaces can be neglected. Based on our analysis, we are able to clarify how dielectric screening influences plasma frequencies in this limit.
Plasmons can be observed by a variety of experi- mental tools including inelastic light scattering 22 , which has been widely used to probe plasmons in semiconductor heterostructures 23 , but also by surface-physics techniques like high-resolution electron-energy-loss spectroscopy 24 , and, more indirectly, angle-resolved photoemission spectroscopy 14 . Double-layer field-effect transistors with a grating gate 25 can also be used to detect plasmons. Coupling between far-infrared light and Dirac plasmons in single-layer graphene has recently been achieved by employing an array of graphene nanoribbons 26 and by performing near-field scanning optical microscopy through the tip of an AFM 27 . This manuscript is organized as follows. In Sect. II we present the model we have used to describe a pair of Coulomb-coupled MD2DESs, and introduce the linearresponse functions which describe collective electron dynamics. In Sect. III we present and discuss our main analytical and numerical results for the dispersion of optical and acoustic plasmons in these systems. Finally, in Sect. IV we present a summary of our main conclusions.

II. MODEL HAMILTONIAN AND RANDOM PHASE APPROXIMATION
We consider two unhybridized MD2DESs separated by a finite distance d and embedded in the dielectric environment depicted in Fig. 1. The two systems are assumed to be coupled solely by Coulomb interactions. The Hamiltonian describing this system reads 28 ( = 1) Here v is the bare Dirac velocity, taken to be the same in the = 1, 2 tunnel-decoupled layers, S is the area of each layer, V (q) is the matrix of bare Coulomb potentials, andρ is the density-operator for the -th layer. The Greek letters are honeycomb-sublattice-pseudospin labels and σ = (σ x , σ y ) is a vector of Pauli matrices. A sum over flavor labels is implicit in Eq. (2) in the case of DLG. The relative strength of Coulomb interactions is measured by the dimensionless coupling constant 3 (restoring for a moment) α ee ≡ e 2 /( v) which has a value ≈ 2.2 in DLG and ≈ 4.4 in Bi 2 Te 3 TIs if we use the respective Dirac velocities v G ≈ 10 6 m/s and v TI ≈ 5 × 10 5 m/s. Several important many-body properties of the Hamil-tonianĤ are completely determined by the 2×2 symmetric matrix χ(q, ω) whose elements are the density-density linear-response functions with Â ,B ω the usual Kubo product. 16 Within the RPA these functions satisfy the following matrix equation, where χ 0 (q, ω) is a 2 × 2 diagonal matrix whose elements χ (0) (q, ω) are the well-known 29-31 noninteracting (Lindhard) response functions of each layer at arbitrary doping n . The off-diagonal (diagonal) elements of the ma- The bare intra-and inter-layer Coulomb interactions are influenced by the layered dielectric environment (see Fig. 1). A simple electrostatic calculation 28 implies that the Coulomb interaction in the = 1 (top) layer is given by where The Coulomb interaction in the bottom layer, V 22 (q), can be simply obtained from V 11 (q) by interchanging 3 ↔ 1 . Finally, the inter-layer Coulomb interaction is given by Notice that in the "uniform" 1 = 2 = 3 ≡ limit we recover the familiar expressions . Previous work on TI thin film and DLG collective modes has assumed this limit, which rarely applies experimentally.

III. COLLECTIVE MODES
The collective modes of the system described by the model Hamiltonian (1) can be determined by locating the poles of χ(q, ω) in Eq. (4). A straightforward inversion of Eq. (4) yields the following condition 32,33 : The collective modes occur above the intra-band particlehole continuum where χ (0) is real, positive, and a decreasing function of frequency. Eq. (8) admits two solutions, a higher frequency solution 29,30,34 at ω op (q) which corresponds to in-phase oscillations of the densities in the two layers, and a lower frequency solution at ω ac (q) which corresponds to out-of-phase oscillations.
The plasmon collective modes of MD2DESs are of special interest because of the ease with which they may be altered by changing the carrier densities in either layer using gates. We note in particular that the carrier densities in different layers can easily differ radically. For this reason we present our results in terms of the total 2D carrier density n = n 1 + n 2 , and the density polarization ζ = (n 2 − n 1 )/n ∈ [−1, 1]: ζ = 1 when the carrier density is non-zero only in the bottom layer (n 1 = 0), while ζ = 0 when the two layers have identical carrier densities (n 1 = n 2 ).

A. Analytical results
In this Section we report on exact analytical expressions for the RPA optical and acoustic plasmon dispersions ω op,ac (q) that are valid in the long-wavelength q → 0 limit where ω op (q → 0) ∝ √ q and ω ac (q → 0) ∝ q.
We start by deriving an exact expression for the RPA long-wavelength acoustic-plasmon group velocity, Following Santoro and Giuliani 33 , we first introduce the power expansion for the acoustic-plasmon dispersion relation, and then define a function In the limit q → 0 the function F (q) has the following Laurent-Taylor expansion where the coefficients f i can be extracted from the analytical expression [29][30][31] for the MD2DES Lindhard function χ (0) (q, ω). For Eq. (8) to be valid we have to require that the coefficients f i vanish identically. The coefficient f −1 depends only on c s and by equating its expression to zero we arrive after some tedious but straightforward algebra at the following equation for x = c s /v, the ratio between the plasmon group velocity c s and the Dirac velocity v: where g s (g v ) are real-spin (valley) degeneracy factors.
In the case of DLG, g s = g v = 2, while in the case of thin-film TIs g s = g v = 1. In Eq. (13)d = dk F is a dimensionless inter-layer distance calculated with k F ≡ 4πn/(g s g v ) and n = n 1 + n 2 , and Eq. (13) can be conveniently solved for x by making the change of variables After some straigthforward algebra we find that Eqs. (15)- (16) are the principle results of this Article. We see from this analytic expression that c s is independent of 1 and 3 and depends only on the barrier material dielectric constant, which in the case of TI thin films is simply the TI bulk dielectric constant. This behavior is a consequence of the out-of-phase character of this collective mode in which the double-layer total charge is locally constant but shifts dynamically between layers. Because TIs tend to have narrow gaps they tend to have large dielectric constants ( 2 ∼ 100 in the case 35 of Bi 2 Te 3 ). Thin-film collective modes will therefore tend to have c s /v values that are quite close to 1 unlessd is very large. (For larged the long-wavelength limit formula, which applies when both qd and q/k F are small, will have a limited range of applicability.) It follows from Eq. (15) that the ratio c s /v is larger than unity for any value of the parameters α ee ,d, ζ, and 2 so that the acoustic plasmon always lies outside of the MD2DES particle-hole continuum. This implies than the acoustic plasmon is strictly speaking never Landau damped at small q. (A similar conclusion was reached previously 33 for the case of conventional 2D electron gases, but was limited to the case of identical density and hence identical Fermi velocity.) For moderate values ofd, however, Eq. (15) predicts a TI thin film sound velocity so close to the top of the particle-hole continuum that it will likely be unobservable because of damping effects not captured by the RPA, and because of disorder, which is always present to some degree. For the case of DLG, on the other hand, we expect that acoustic plasmon collective modes will be well defined. This is particularly true in the case of DLG Here kF is the Fermi wave vector evaluated at the total density n = n1 +n2, i.e. kF = 4π(n1 + n2)/(gsgv). The parameters we have used to calculate the curve labeled by "ζ = 0" are: gs = gv = 2, n1 = n2 = 5 × 10 12 cm −2 , αee = 2.2, d = 3.35Å, 1 = 2 = 1, and 3 = 3.9. These parameter values correspond to the case of two graphene layers on SiO2 that are decoupled by rotation. The data labeled by "ζ = 0.8" have been calculated by setting n1 = 1 × 10 12 cm −2 and n2 = 9 × 10 12 cm −2 with the same values for the other parameters. The solid lines plot ω = csq for the ζ = 0 and ζ = 0.8 cases with the plasmon group velocity cs calculated from the analytical result, Eq. (15). These numerical results confirm the validity of our analytic result for cs and the importance of accounting for the delicate dependence of the long-wavelength Lindhard function on ν = ω/(vq). The dashed line plots the upper-bound of the intra-band electron-hole continuum, ω = vq.
with a small number of layers of BN as barrier material. When the BN barrier layer is very thin, the use of macroscopic dielectric parameters to characterize its screening properties is approximate; in that case measurement of the acoustic plasmon group velocity combined with Eqs. (15)-(16) would allow the effective value of 2 to be determined experimentally.
We note that an analytic result for c s was reported previously in Ref. 20 [see their Eq. (5b)] for the special case of DLG embedded in a uniform dielectric, i.e. for 1 = 2 = 3 . In our notation, their result reads This equation is evidently different from Eq. (15) above. We believe that Eq. (15) is the correct RPA result for the acoustic-plasmon group velocity and that Eq. (17) is incorrect. The difference is due to the singular behavior of the Lindhard function χ (0) (q, ω) as a function of wave vector q and frequency ω in the region in which both these quantities are small. (See Sect. 4

.4.3 of Ref. 16.)
In particular, the limit of χ (0) (q, ω) for q → 0 and ω → 0 depends on the ratio ν = ω/(vq), i.e. on the direction along which the origin of the (q, ω) plane is approached: different limits are obtained for different values of ν. In an acoustic plasmon, the ratio ν approaches a constant as q → 0 and thus the limit of χ (0) (q, ω) which matters is the one in which q → 0 while the ratio ω/q is kept constant. This is the limit we have taken 33  A careful comparison between our analytical prediction in Eq. (15) and the result obtained by the brute-force numerical solution of Eq. (8) is shown in Fig. 2. We clearly see that Eq. (15) compares very well with the full numerical result.
The analytical analysis of the long-wavelength optical plasmon mode is simpler since this mode satisfies ω op (q) ∝ √ q for q → 0 and therefore occurs at ν = ω/(vq) → ∞. We obtain an analytic result using the well-known high-frequency (ω vq and ω 2ε F, ) dynamical limit of χ (0) (q, ω): with ε F, = vk F, = v 4πn /(g s g v ). Using Eq. (18) in Eq. (8) we immediately find with¯ = ( 1 + 3 )/2. Note that Eq. (19) does not depend on the inter-layer distance or on the dielectric constant 2 , but only on the average¯ between top and bottom dielectric constants. Notice also that, in the limit n 1 → 0 (ζ = 1), Eq. (19) reduces to the well-known plasmon frequency in a single-layer graphene sheet 29,30,34 with electron density n 2 . This expression applies for qd 1, in which case the entire double-layer MD2DES acts in the optical plasmon mode like a single conducting layer at the interface between dielectric media characterized by constants 1 and 3 .

B. Numerical results
In this Section we briefly report some representative numerical results for the optical and acoustic plasmon dispersion relations obtained by solving Eq. (8), discussing first DLG and then TI thin films.
In Fig. 3 we illustrate the typical properties of DLG plasmon modes for the case with the smallest MD2DES separation, two adjacent layers on a SiO 2 substrate ( 1 = 2 = 1 and 3 = 3.9) that are weakly hybridized e.g. because of a twist between their orientations 36 . Fig. 3a) is for a symmetric system with the same electron concentration on the two layers (ζ = 0), while Fig. 3b) refers to a system with a 50% density imbalance. The characteristic behaviors ω op (q) ∝ √ q of the optical plasmon and ω ac (q) ∝ q of the acoustic plasmon are clearly visible. The collective modes are not Landau damped when they appear in the gap between intra-band and interband particle-hole continua. When the two layers have different densities, their particle-hole continua are different and the gap is smaller for the lower density layer. For adjacent but twisted DLG systemsd is small even when the carrier density is large (d ≈ 0.2 in Fig. 3). It follows that qd is small and the two MD2DESs are strongly coupled over the entire relevant frequency regime. In this smalld example the acoustic plasmon frequency is close to the particle-hole continuum because the capacitive energy associated with charge sloshing between the layers is proportional to the small layer separation.
In Fig. 4 we illustrate the strength of plasmon decay by emission of single electron-hole pairs (Landau damping). Notice that Landau damping occurs when the curves ω op,ac (q) in Fig. 3 hit the inter-band electron-hole continuum of the layer with lower density (layer "1" in our convention). The larger ζ, the sooner this happens. In particular, in the limit in which layer "1" is neutral (ζ = 1), Landau damping is present from vanishingly small wave vectors: damping of the optical plasmon excitation associated with electrons in the high-density layer starts at arbitrarily small wave vectors since decay can easily occur via the emission of inter-band electron-hole pairs in the neutral layer. The many-body properties of two or more decoupled graphene layers can thus be strongly affected by inter-layer Coulomb interactions, even by apparently innocuous geometric features such as the presence of a nearly-neutral layer.
In Fig. 5 we compare optical and acoustic plasmon dispersions for DLG and TI thin-film systems. For the TI thin-film case we have chosen the following parameters: i) 1 = 1, 2 = 100 (this roughly corresponds to the dielectric constant of Bi 2 Te 3 ), and 3 = 4.0; ii) a total electron density on the top and bottom surface states of n = 10 13 cm −2 ; and iii) a thickness of the TI slab of d = 6 nm, corresponding to a six quintile layer MBE-grown Bi 2 Te 3 film (d ≈ 6.7). The DLG example has the same total density and layer separation (d ≈ 3.4; the difference ind in the two cases stems from the g s /g v spin/valley degeneracy factors) and dielectric constants 1 = 1 and 2 = 3 = 4.0, corresponding to two graphene layers separated by approximately 15 BN layers and lying on a BN substrate. In both cases we see that a crossover occurs at intermediate values of q between strong (small q) and weak (large q) coupling of the two collective modes. In the TI case the higher frequency optical plasmon mode deviates much more strongly from simple √ q behavior at this crossover because strong dielectric screening by the TI bulk suppresses the singlesurface plasmon mode. [Note however that the effective dielectric constant for this limit is ( 2 + 1,3 )/2 rather than 2 as used in Ref. 21.] The acoustic plasmon mode of the TI thin film case is, on the other hand, strongly suppressed in the strong-coupling limit, as discusses earlier, and has a velocity much closer to the bare Dirac velocity than in the corresponding DLG case.  Fig. 3a). Note that, within RPA, m χ (0) 1 (q, ωop,ac(q)) is identically zero for wave vectors q up to a critical value q op,ac at which ωop,ac(q) hits the inter-band electron-hole continuum associated with the low-density layer. Since data in this panel correspond to ζ = 0, top-layer and bottom-layer Lindhard functions are identical. Panel b) Same as in panel a) but for ζ = 0.8 (n1 = 1 × 10 12 cm −2 and n2 = 9 × 10 12 cm −2 ). Note that the q op,ac decreases with increasing ζ becoming zero in the limit ζ → 1. Since in this panel ζ = 0 we have plotted both top-layer (low-density) and bottom-layer (high-density) Lindhard functions.

IV. DISCUSSION AND CONCLUSIONS
We have presented an analysis of the electronic collective modes of systems composed of two unhybridized but Coulomb-coupled massless-Dirac two-dimensional electron systems (MD2DESs) separated by a vertical dis- tance d. The primary example we have in mind is topological insulator (TI) thin films, which are always described at low energies by this type of model because topologically protected MD2DESs always appear on both top and bottom surfaces. Also of interest are closely related systems, which we refer to as double-layer graphene (DLG) systems, containing two graphene layers that are weakly hybridized either because they are rotated relative to each other or because they are separated by a dielectric barrier layer. Importantly, we allow for a general dielectric environment in which the material above the top MD2DES layer ( 1 ), between the two layers ( 2 ), and below the bottom MD2DES layer ( 3 ) are all allowed to have different dielectric constants. In the case of TI thin film 2 is the bulk dielectric constant of the TI which is expected to have large values. The carrier collective modes of MD2DESs are expected to be most robust in the gap between intraband and interband particle-hole excitations.
The double-layer systems of interest quite generally have two collective modes which in the limit of small qd involve density-fluctuations in the two-layers that are strongly coupled, and in the limit of large qd weakly coupled single-layer plasmons. One key parameter which controls collective mode properties is the dimensionless product k F d ≡d. Small values of k F d imply that the layer separation is smaller than the typical distance between electrons within a layer and that collective modes at all values of q up to ∼ k F are strongly coupled combinations of the two individual layer density-fluctuation contributions. For large k F d a crossover occurs for q ∈ (0, k F ) between strongly and weakly coupled double-layer collective modes. Both small and large values ofd are achievable in samples where disorder plays an inessential role in both DLG and TI thin film cases.
Our study focuses on the long-wavelength limit in which both qd and q/k F = qd/d are small. We have derived analytic expressions for both frequencies of both the low-energy linearly dispersing acoustic plasmon mode ω ac (q) and for the high-energy optical plasmon mode ω op (q) which has √ q dispersion at long-wavelengths. In this limit we find that ω ac (q) − vq ∝ 1/ 2 whereas ω op (q) ∝ 2/( 1 + 3 ); i.e. the separation of the acoustic plasmon mode from the upper edge of the intra-band particle-hole continuum is very strongly suppressed by a large bulk TI dielectric constant, whereas the coupled double-layer plasmon mode is unaffected. This doublelayer optical plasmon behavior contrasts with that of a large qd single-surface plasmon mode which has a frequency proportional to 2/( 2 + 1,3 ). The long-wavelength limit of ω ac (q) is sensitive not only to the energy associated with inter-layer charge sloshing but also to its microscopic kinetics as captured by the singular sensitivity of the MD2DES Lindhard function to ω/(vq). By carefully accounting for this dependence we are able to correct a previous analytic expression in a way that is quantitatively particularly important in the TI thin film (large 2 ) case. Double-layer collective mode coupling plays an important role in MD2DES correlations whend is small. Even whend is large, strongly-coupled small qd = qd/k F modes will often be experimentally accessible and may play an important role in graphene multi-layer or TI based plasmonics. The analytic results derived in this paper can be used to readily anticipate how these modes depend on system parameters.
From the more theoretical point of view, it will be intriguing to study physical properties of plasmons in Coulomb-coupled MD2DESs beyond the random phase approximation by employing e.g. many-body diagrammatic perturbation theory 37 . While this manuscript was being finalized for publication, we became aware of a study of optical and acoustic plasmons in double-layer graphene 38 . The authors of this work present extensive numerical results for the "uniform medium" limit ( 1 = 2 = 3 ) and discuss the relation between (longitudinal and transverse) plasmons and near-field amplification.