Plasmon Damping Rates in Coulomb-Coupled 2D Layers in a Heterostructure

The Coulomb excitations of charge density oscillation are calculated for a double-layer heterostructure. Specifically, we consider two-dimensional (2D) layers of silicene and graphene on a substrate. From the obtained surface response function, we calculated the plasmon dispersion relations, which demonstrate how the Coulomb interaction renormalizes the plasmon frequencies. Most importantly, we have conducted a thorough investigation of how the decay rates of the plasmons in these heterostructures are affected by the Coulomb coupling between different types of two-dimensional materials whose separations could be varied. A novel effect of nullification of the silicene band gap is noticed when graphene is introduced into the system. To utilize these effects for experimental and industrial purposes, graphical results for the different parameters are presented.


Introduction
A huge number of researchers from various disciplines have been showing their interest in new materials, silicene especially, after the development of its fabrication process in 2012 [1]. Because of its exceptional potential applications in electronic and optoelectronic devices, many industries are making substantial investments to harness its properties. Additionally, before making investments for commercial gain, both theoreticians and experimentalists have been exploring this material for many years. A credit of foremost importance goes to Takeda and Shiraishi [2], who, in 1994, dealt with the atomic and electronic structure of the material for the first time. These authors calculated the band structure of silicon in the corrugated stage having optimized atomic geometry. This work, though very novel, did not receive the attention it deserves until 2004, when single-layer carbon atoms named graphene were fabricated in the laboratory from graphite by Novoselov et al. [3]. Their research not only validated the stability of two-dimensional (2D) material but also opened the door for new research on thin film materials, silicene being one of them.
Both silicene and graphene were studied in parallel. The former has a buckled crystal geometry, whereas the latter has a honeycomb planar geometry. Due to this, differences arise between them. Ab initio calculations showed that the bandgap of silicene is electrically tunable [4][5][6], which is an advantageous property for designing a field effect transistor that works at room temperature. Another distinct difference between these two materials is the strength of the spin-orbital coupling (SOC), which is very weak in graphene. Consequently, the quantum spin Hall effect occurs at extremely low temperatures [7,8]. In contrast to this, silicene displays quantum spin Hall effect at temperature 18 K, far higher than that for graphene.

Theory
In this section, we analyze the plasmonic behavior of the heterostructure consisting of graphene and silicene together for which we employ the low-energy form of the Hamiltonian near the K point in the Brillouin zone. One significant difference between the Hamiltonian of graphene and silicene is that a small band gap, ∆ is present in the silicene energy band structure, which is due to spin-orbit coupling and applied external electric field. This band gap is not seen in intrinsic graphene.

Silicene
We now briefly describe the case pertaining to silicene whose Hamiltonian in the continuum limit is given by: whereτ x,y,z and σ x,y,z are Pauli matrices corresponding to two spin and coordinate subspaces, ξ = ±1 for the K and K valleys, v F (≈ 5 × 10 5 ) m/s is the Fermi velocity for silicene [5,54], k x and k y are the wave vector components measured relative to the K points. The first term represents the low-energy Hamiltonian, whereas the second term denotes the Kane-Mele system [7] for intrinsic spin-orbit coupling with an associated spin-orbit band gap of ∆ so (of order 5 to 30 meV can could reach up to 100 meV) [55]. The last term in the expression describes the sublattice potential difference that arises from the application of a perpendicular electric field. Equation (1) for the Hamiltonian is a block diagonal in 2 × 2 matrices labeled by valley (ξ) and spin σ = ±1 for up and down spin, respectively. These matrices are given by [55]: This gives the low-energy eigenvalues as: where ∆ σξ = |σξ∆ so − ∆ z |.

Graphene
The low-energy model Hamiltonian for monolayer graphene is similar to that in Equation (2) with the diagonal terms replaced by zero and ξ labeling the valley. In this regime, the Hamiltonian for intrinsic graphene is given by [56]: with the linear energy dispersion, E k = ±hv f |k| in either valley.

Polarization Function: Π(q, ω)
Considerable work has been completed on the dynamical properties involving the use of the dielectric function (q, ω) of various types of free-standing 2D systems [57][58][59] under different conditions. These include temperature effects [60,61], the role of an ambient magnetic field for the 2D electron gas (2DEG), graphene, silicene [62] and the dice lattice [58]. For a single 2D layer, one can extract the plasmon dispersion relation and damping rate by employing the dielectric function. However, the situation is more complicated for a multi-layer heterostructure that relies on knowledge of the surface response function, which we have presented in detail below. However, in either case, we need to calculate the polarization function obtained in the random phase approximation (RPA). For a 2D layer surrounded by a medium with dielectric constant b , the dynamic dielectric function is given by: where V(q) = 2πe 2 4π 0 b q is the Coulomb interaction potential and 0 is the permittivity of free space, q is the wave vector and e is the electron charge. The polarization function is an important quantity in calculations of the transport, collective charge motion and charge screening properties of the material. In the one-loop approximation, the polarization function for gapped graphene is given by [63]: where the angle between k and k + q is θ k,k+q . At zero temperature, the Fermi function f 0 (z) is just a step function. The analytical expression for the polarization function for the silicene and graphene monolayer is given by Tabert et al. [62] and Wunsch et al. [57], respectively. We now turn our attention to a crucial consideration in this paper regarding the structure consisting of a silicene layer, a graphene layer and substrates as depicted in Figure 1. The dielectric constants 1 and 2 are related to the space between the two layers and the semi-infinite region underneath the lower layer. Thus, we assume that there is no material above the upper layer whose susceptibility is χ 1 (q, ω), and it is always assumed to be a vacuum above this top layer. The two 2D layers may be identical or different, possessing different material properties (graphene or silicene in our case), which is reflected in their energy dispersions though the presence or absence of a band gap. The two layers could also have different or equal doping levels (Fermi energies). By employing the boundary condition of continuity of the electrostatic potential and the discontinuity of the electric field across the interface separating two media, we solved for the various coefficients appearing in the potential. Consequently, the result for the surface response function g(q, ω) gives the required conditions for the plasmon dispersion for our case, namely: Here, φ < (z), φ > (z) and φ 1> (z) correspond to the electrostatic potential of regions (I),(II) and (III), respectively, as shown in Figure 1. In order to conduct numerical computation, we make use of linear response theory, for which we have the charge density, i for convenience, we obtain the solution of these equations for different coefficients, leading to: where 1 (ω) is the dielectric function of the substrate between layers "1" and "2", χ 1 and χ 2 correspond to the susceptibilities of these two layers and d is the thickness of the substrate. This substrate thickness is in the order of the wavelength of light considered, for a visible light it could be of the order of a few hundred nanometers. For a thick substrate, the thickness could go up to micrometer in size, and accordingly, the plasmonic mode is modified. The plasmon dispersion equation is obtained from the solutions of D(q, ω) = 0, which we solve below. We note that when we set χ 2 = 0 and take the limit d → ∞, Equation (8) gives the well-established form [64]: which is the surface response function for a 2D layer embedded in a medium whose average background dielectric constant is b = ( 1 + 1)/2. The plasma resonances, which Equation (10) gives from its poles, are in agreement with the zeros of the dielectric function in Equation (5).

Damping Rate
We now turn to a critical issue in this paper, which concerns the rate of damping of the plasmon modes by the single-particle excitations. If this rate of damping for a plasmon mode with frequency Ω p is denoted as γ, then D(Ω p + iγ, q) = 0 is the complex frequency space. Carrying out a Taylor series expansion of both the real and imaginary parts, we have: Therefore, setting the function in Equation (11) equal to zero, we obtain γ to the lowest order as: With these formal results, we now evaluate the plasma spectra for a double layer heterostructure. The expression shows the dependence of γ on the imaginary part of D(Ω p ) and the Real part of D(ω), which in turn, are dependent on the type of layer and the substrate considered. Eventually, we can infer that the viability of plasmon modes can be tuned by the dielectric substrate thickness and by the choice of 2D layer. In addition, the rate of decay also helps us in maintaining the intensity and the frequency of the obtained plasmon mode. This could have great impact in the development of the quantum information sharing technology and the data storing devices.

Numerical Results and Discussion
In our numerical calculations, energy is scaled in units of E F is equivalent to ∼60 meV. From the preceding discussion, in Section 2, it is clear that the plasmon modes for any system are given by the zeros of the dielectric function obtained from Equation (9). Thus, making use of this dispersion equation, we computed the plasmon mode dispersion relation for a heterostructure based on graphene and silicene with/without a substrate. For this, we first obtained a graphical result for graphene, as shown in Figure 2. One can clearly see that a single branch plasmon mode originates from the origin in q − ω space, which increases monotonically and is subject to Landau damping when the plasmon mode reaches the interband particle-hole excitation region. Figure 2 is plotted for monolayer graphene embedded in material with background dielectric constant 1 = 2 , which we set equal to 1 for simplicity. The situation is similar ( 1 = 2 = 1) for Figures 3-7. In Figure 2, the plasmon branches for two values of the Fermi energy are shown in panels (a,b). The damping rates of these plasmon modes are demonstrated in panel (c,d), where it is distinctly shown by an arrow pointing at the boundary of the region where Landau damping takes place. The rate of decay for both types of graphene are monotonically increasing, signifying that the deeper into the single particle excitation region plasmon mode enters, the larger the rate of decay becomes. This implies that the lifetime of the plasmon mode is decreased in the same manner.  panels (b,d)). Two top panels (a,b) demonstrate the plasmon dispersion (either damped or undamped obtained as Re( (q, ω)) = 0, the lower plots (c,d) describe the corresponding damping rate along the plasmon branches, also calculated and shown as insets (i1) and (i2).    Going next to the case when we have a structure with two graphene layers together separated by various distances, a set of plots (as shown in Figure 3) is obtained with two branches of plasmon modes originating from the origin in the q − ω space. One can see that when the graphene sheets are brought closer, the plasmon modes move further apart, signifying weak interactions between the modes. In Figure 3a Figure 4, which shows that the plasmon decay for the lower plasmon branch always starts at a larger wave vector value in comparison to the upper plasmon branch. As the distance of separation is increased, the two plasmon branches come closer and their decay also starts from the same value of the wave vector and the rate of decay is almost the same in value.

Silicene-Silicene plasmon
Now, in addition, we carried out an investigation of the plasmon modes and their decay rate for the structure with two silicene layers for various separations. The graphical results for these calculations are shown in Figure 5 where we again have two plasmon modes originating from the origin of the q − ω plane. As in the case of a two-graphene-layer structure, we again notice a similar effect on two plasmon branches coming closer to each other when their separation increases. This is demonstrated in Figure 5a Figure 6 shows that the upper plasmon mode does not decay at all and the lower plasmon branch decays after reaching a critical wave vector. This behavior is due to the presence of a band gap for silicene, resulting in an opening in the single particle excitation region, which provides a larger area in q − ω space for the plasmon mode to survive. The upper plasmon branch in this case has a larger space and is more likely to self-sustain for a longer period without damping. On the other hand, the lower plasmon branch enters the intraband single-particle excitation region where it decays. The rate of decay starts from a critical value of the wave vector and the magnitude of this decay rate monotonically increases.
A comparison of plasmon modes and their decay for graphene-graphene and silicenesilicene structures is shown along with the single-particle excitation regions in Figure 7. The figure in panel (a) of Figure 7 shows that two plasmon modes that stem from the origin of the frequency-momentum space increase steadily and decay when it reaches the boundary of the single-particle excitation region. Corresponding red and blue lines are drawn to further clarify the point where the actual decay begins. The dark triangular region is the area where the plasmon mode survives without Landau damping and mathematically, in this region, the imaginary part of the polarization function of graphene is zero. This means that the plasmon mode has self-sustaining oscillations. The green region where the imaginary part of the polarization function is nonzero is the single particle excitation region where the plasmon mode decays into particle-hole mode. The corresponding decay rate figure, below this panel, shows that the rate of decay for the upper plasmon branch is greater and its critical wave vector is smaller compared to the lower plasmon branch.
Similar plots for silicene-silicene structures were demonstrated in Figure 7b where one can see the opening of a gap in the single-particle excitation region yielding two parts, which is a significant effect arising from the band gap. The imaginary part of the polarization function in this gap region is zero where the plasmon mode can sustain its oscillation for a long time. The upper breakaway region is a single-particle excitation region due to interband transitions of electrons from the valence to the conduction band and the lower breakaway region is the intraband single-particle excitations region, which is due to transitions within the same band from below to above the Fermi level. In Figure 7b, two plasmon modes originate from the origin as demonstrated in the figure. The upper plasmon mode survives without damping over a wide range of wave vectors and the plasmon branch enters the gap created by the opening within the single-particle excitation region. The corresponding decay rate appearing below the plasmon dispersion shows that the upper plasmon branch does not decay at all, whereas the lower plasmon branch with the part in closer contact with the intraband single-particle excitation region has a small plasmon decay rate as illustrated in the corresponding figure in the panel of Figure 7d below. As the plasmon mode rises, it is separated from the single-particle excitation region where the decay rate is zero and as it moves further away from the origin, the plasmon branch becomes closer to the single-particle excitation region where we notice the Landau damping again. Correspondingly, the decay rate increases monotonically and reaches a maximum before dropping down, indicating the reappearance of an undamped plasmon branch at a larger value of the wave vector. Another noticeable effect observed here is the closeness of the plasmon branches and the plasmon decay rate, which can be altered by altering the layer separation; this effect may be used as another plasmon mode tuning parameter.
These two branches appear as a result of the Coulomb interaction between the two layers, which couple the plasmon excitations arising on each layer. Therefore, the resulting plasmons are physically similar to two coupled oscillators. The obtained plasmon branches are defined as acoustic (lower frequencies) and optical (higher frequencies), or as in-phase and out-of-phase [65]. The number of branches is equal to the number of layers [66], or, more precisely, the number of separate plasmon excitations in them. However, one cannot attribute one branch to the first layer and the other to the second one.
We should also mention that the frequency of the optical plasmon branch depends on the distance between the layers. One can easily verify this by analyzing how far the branch moves from the main diagonal ω = v F q when the distance between the layers is increased. However, this dependence is not as dramatic as for the acoustic plasmons, which also change their shape when the distance between the layers is increased, and the Coulomb coupling is faded away. Finally, we should say that the main subject of the present paper is the plasma branches for a system of different layers and their damping rates, which are affected by the distance between the layers and the type of substrate in between.
To extract more information about the plasmonic behavior, in Figure 8a, we have presented the figure to show the result highlighting the changes in the plasmonic nature for a structure with different types of layers and substrates. In Figure 8a, we demonstrate the plasmon mode for a structure with silicene and graphene separated by a distance of F ) −1 with the vacuum in between. One can observe two plasmon modes originating from the origin in q − ω space. p /ω 2 , for the substrate. Panels (c,d) represent the plasmon damping rate for the structure corresponding to the plasmon modes in (a) and (b), respectively. The bottom substrate is given by 2 = 1 for both cases.
A special effect of overcoming the single-particle excitation region of silicene by the single-particle excitation region of graphene is observed, which causes the shortening of the plasmon branches that used to be there for the silicene-silicene structure. As soon as the plasmon branch reaches the particle-hole mode region, the plasmon mode decays by replacing one silicene layer with a graphene layer in the silicene-silicene structure. The effect, due to the band gap in silicene, is just nullified. In other words, the plasmon modes in the regime, where they used to have plasmon modes before, now do not have them, because the plasmon mode decays into the particle-hole mode of graphene. Furthermore, the result of adding a substrate between the silicene and graphene layer is illustrated in panel (b) of Figure 8. In this case, we could see a new plasmon branch originating from the bulk plasma frequency and one plasmon branch originating from the origin. Here, due to the presence of a substrate, the lower plasmon branch bends sharply toward the intraband single-particle excitation region where it decays causing complete disappearance. The upper plasmon branch and the plasmon from the bulk plasmon frequency become closer and move toward the interband single-particle excitation region where they become damped.
The results in Figure 8 correspond to a graphene layer located on the top and the silicene layer located below it in the heterostructure. The order in which the layers are placed affects the plasma dispersions only if asymmetric substrates are involved. In contrast, the plasmon dispersion relation for a pair of Coulomb-coupled layers embedded in a uniform background is determined by a 2 × 2 determinant equation [67], which is symmetric to switching the layers. This is what our Equation (9) is reduced to when 1 = 2 = 1 and the background dielectric constants are independent of the frequency ω.
These new effects on the plasmon branches in this type of structure were not reported previously. Results of this type help develop electronic and quantum computing devices where knowledge of the plasmonic behavior of materials and their damping nature is very essential.

Concluding Remarks
In summary, we have investigated the key properties of the plasmonic mode and damping for different combinations of graphene and silicene layers. The effect of the addition of substrate in between the two layers is further analyzed. This resulted in the development of a novel technique of tuning the plasmon excitation mode associated with the two dimensional systems. In our system, a complete new plasmon branch emerges from the bulk plasmon frequency, this would be very helpful in engineering modern computing devices.
Another discovery we have from this study is the disappearance of the lower plasmon branch and the suppression of the silicene band gap effect. Along with the study of interesting features in the plasmon modes, we have also developed an approach for calculating the decay rates for the plasmons due to Landau damping by the particle-hole modes.
The principal goal of our investigation and the main results of our work are concerned with the understanding of plasmonic nature and analyzing their Landau damping rates in a multi-layer structure. In contrast to Ref. [67], which includes only the long-wavelength limit for graphene without any discussion of gapped materials (such as silicene), our work is concerned with a thorough and detailed investigation of plasmon and damping rates for finite-value wave vectors and energies.
Additionally, our results infer that the number of plasmon branches emerging from the origin can be varied by choosing the number of the 2D layer. In brief, we can say that our study gives an important idea about the plasmonic behavior of a graphene-silicene-based heterostructure, which would be very helpful in carrying out further studies of other types of heterostructure, including a variety of low dimensional materials.