Controlling plasmon modes and damping in buckled two-dimensional material open systems
Andrii Iurov, Godfrey Gumbs, Danhong Huang, Liubov Zhemchuzhna
CControlling plasmon modes and damping in buckled two-dimensional material opensystems
Andrii Iurov ∗ Center for High Technology Materials, University of New Mexico,1313 Goddard SE, Albuquerque, NM, 87106, USA
Godfrey Gumbs
Department of Physics and Astronomy, Hunter College of the City University of New York,695 Park Avenue, New York, NY 10065, USA andDonostia International Physics Center (DIPC), P de Manuel Lardizabal, 4, 20018 San Sebastian, Basque Country, Spain
Danhong Huang
Air Force Research Laboratory, Space Vehicles Directorate, Kirtland Air Force Base, NM 87117, USA andCenter for High Technology Materials, University of New Mexico,1313 Goddard SE, Albuquerque, NM, 87106, USA
Liubov Zhemchuzhna
Department of Physics and Astronomy, Hunter College of the CityUniversity of New York, 695 Park Avenue, New York, NY 10065, USA (Dated: October 31, 2018)Full ranges of both hybrid plasmon-mode dispersions and their damping are studied systematicallyby our recently developed mean-field theory in open systems involving a conducting substrate and atwo-dimensional (2D) material with a buckled honeycomb lattice, such as silicene, germanene, and agroup IV dichalcogenide as well. In this hybrid system, the single plasmon mode for a free-standing2D layer is split into one acoustic-like and one optical-like mode, leading to a dramatic change in thedamping of plasmon modes. In comparison with gapped graphene, critical features associated withplasmon modes and damping in silicene and molybdenum disulfide are found with various spin-orbitand lattice asymmetry energy bandgaps, doping types and levels, and coupling strengths between2D materials and the conducting substrate. The obtained damping dependence on both spin andvalley degrees of freedom is expected to facilitate measuring the open-system dielectric propertyand the spin-orbit coupling strength of individual 2D materials. The unique linear dispersion of theacoustic-like plasmon mode introduces additional damping from the intraband particle-hole modeswhich is absent for a free-standing 2D material layer, and the use of molybdenum disulfide witha large bandgap simultaneously suppresses the strong damping from the interband particle-holemodes.
PACS numbers: 71.45.Gm, 73.21.-b, 77.22.Ch
I. INTRODUCTION
Plasmons, or self-sustained electron density oscillations, represent a broad interest and have become an importantsubject for both traditional and recently discovered two-dimensional (2D) materials.
Due to the possibility offabricating complex stacking layered structures, low-dimensional materials have become very attractive for novelquantum-electronic devices. Discovery of graphene, successfully fabricated in 2004, has initiated a number of newtransport and optical studies due to its unusual electronic properties originating from the relativistic lineardispersion of its energy bands. In particular, graphene plasmonics has quickly become one of the actively-pursuedresearch focuses. As an example, novel optical devices in a wide-frequency range showed significant improvement inall the crucial device characteristics (see Ref. [12] and the references therein). At high energies, the π -bond plasmonsin graphene demonstrate both anisotropy and splitting. Moreover, a plasmonic nanoarray coupled to a singlegraphene sheet was found to have significant enhancements in resonant Raman scattering, as well as in spectralshifts of diffractively-coupled plasmon resonances. Plasmon resonances can be employed for investigating chemicalproperties of either graphene or another adjacent bulk surface, displaying the maturing capabilities of graphene-basedplasmonics.
A junction between graphene and metallic contacts could also be used to design and fabricate ahigh-performance transistor. Consequently, exact knowledge about the plasmon dispersions and their mode dampingin hybrid devices based on the newly discovered 2D materials seems absolutely necessary for the full development ofthese novel devices and the extension of their applications.
Historically, it is well known that the time evolution of a plasmon excitation in a closed system (e.g., free-standing a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n
2D layers) is determined by two-particle Green’s functions in many-body theory, from which we are able to deriveboth the plasmon dispersion and the plasmon dissipation (damping) rate. However, in an open system, thetime evolution for electronic excitations becomes much more involved since it depends strongly on the Coulombinteraction with the environment (e.g., electron reservoirs). As an example, the classical and quantum dynamicalphenomena in open systems include tunnel-coupling to external electrodes, optical-cavity leakage to free space, and thermal coupling to heat baths. The coupling to an external reservoir is usually accompanied by extra dissipationchannels. Many dynamical energy-dissipation theories for open systems are based on the so-called Lindblad dissipativesuperoperator. In spite of the obvious advantages, such as being fast and tunable, in designing graphene-based devices, creating asizable energy bandgap has become an important issue for the practical use of graphene for transistors. The reasonbehind this issue is electrons in gapless graphene may not be confined well by an electrostatic gate voltage or blockedby an energy barrier.
Scientists suggested many approaches for opening an energy gap around (cid:119) . and graphene nanoribbons with quantized transverse wave vectors, or even byshining an intensive and polarized irradiation to dress electrons in graphene. In this respect, experimental implementation of a 2D lattice with a sizable spin-orbit coupling seems to be animportant advancement. Buckled structures, such as silicene and germanene in which atoms are displaced out ofthe plane due to sp hybridization, exhibit significant asymmetry with respect to the a − and b − sublattices. Thisleads to a new type of bandgap which is tunable by applying a perpendicular electric field. Physically, silicene andother buckled honeycomb lattices have been modeled successfully by introducing the Kane-Mele type Hamiltonian with an intrinsic spin-orbit energy gap (1 . − . This includes recently synthesized germanene with a considerably larger spin-orbit coupling and the bandgap of 24 −
93 meV. The Hamiltonian, energy dispersion,and related electronic properties of germanene are qualitatively similar to those of silicene, although the Fermi velocity,the bandgap induced by spin-orbit coupling and the buckling height in germanene are still different in magnitudes.We, therefore, believe that our previous theory for silicene can equally be applicable to Ge-based hybrid structures.Germanene layers have been fabricated by molecular beam epitaxy on Ag(1 1 1) surfaces through deposition onh-AlN and investigated by x-ray absorption spectroscopy. Here, h-AlN was used to create an insulating buffer layerbetween germanene and its metal substrate. The measured lattice constant is in good agreement with the theoreticalpredictions for a free-standing Ge layer. In addition, a thorough experimental study for the density of states ofgermanene which had been synthesized on Ge / Pt crystals at finite temperatures was performed by using a scanning-tunneling electron microscope. The obtained virtually perfect linearly dependent density of states is clearly a proofof a 2D Dirac system. Furthermore, Friedel oscillations were not observed, implying the possible Klein paradox withinthe considered system.Recently, there has been a number of reported studies on microscopic electronic properties, insulating regimes andtopologically protected edge states within a certain range of applied electric fields, as well as spin- and valley-polarizedquantum Hall effects.
Here, a crucial feature for a buckled structure is the occurrence of the topological-insulator(TI) properties when the external electrostatic field is relatively low and the resulting field-related energy gap ∆ z becomes less than the intrinsic spin-orbit one. If the two gaps are equal, on the other hand, the system turns into aspin-valley polarized metal (VSPM) since the gap for one of the two subbands will close. For a strong electric field,the system behaves just like a regular band insulator (BI). One expects that the TI phase can show some uniqueelectronic properties.
Indeed, the properties of a TI, and the energy bandgap of Si as well, are experimentallyfound tunable by an in-plane biaxial strain.In addition to silicene and germanene, another atomically thin 2D material with a great potential for applica-tions in electronic devices is MoS (monolayer molybdenum disulfide, or ML-MDS) which is a prototype of a metaldichalcogenide. The first-principle studies of the electron structure of this material has predicted a hybridization ofthe d -orbitals of molybdenum atoms with the p -orbitals of sulfur atoms, giving rise to a two-band continuum modelfor MoS monolayer. This model further indicates the existence of two spin and valley degenerate subbands witha very large direct bandgap (1 . . It is important to note that the electronicproperties of a ML-MDS are drastically different from its bulk samples with an indirect bandgap (cid:119) . at the room-temperature exceeds 200 cm V − s − and alsoacquires an ultralow standby power dissipation. Together with its direct bandgap, this makes ML-MDS an excellentcandidate for next generation field-effect transistors. In strong contrast to the bulk, the MoS monolayer can emitlight efficiently, and a possible high-temperature superfluidity, as a two-component Bose gas, has been predicted ina system involving a transition metal dichalcogenide bilayer. Experimentally, germanene has already been successfullysynthesized on MoS substrate. This leads to a huge advantage in comparison with the germanene synthesis on a
2D layer (silicene)spacer layer semi-infnite conductor
FIG. 1: Schematics of a hybrid plasmonic structure, including a 2D layer (silicene, germanene, MoS etc.) interacting with asemi-infinite conducting substrate, separated from the layer by a distance a . metallic substrate in which germanene is strongly hybridized with metals, leading to two parabolic bands, instead oflinearly dispersing bands, at its six K points.Our main focus in this paper is studying electron interactions in (2DMOS), i.e., ananoscale hybrid structure of a 2D layer (graphene with or without energy gap, silicene, germanene, a transition-metal dichalcogenide monolayer), as illustrated by Fig. 1. The most distinguished feature of such an open system is thedynamically-screened Coulomb interaction between electrons in graphene and the conducting substrate. Screeningis included by calculating the nonlocal frequency-dependent inverse dielectric function K ( r , r (cid:48) ; ω ) which is related tothe regular dielectric function (cid:15) ( r , r (cid:48)(cid:48) ; ω ) by (cid:82) d r (cid:48) K ( r , r (cid:48) ; ω ) (cid:15) ( r (cid:48) , r (cid:48)(cid:48) ; ω ) = δ ( r − r (cid:48)(cid:48) ). The resonance occurring in K ( r , r (cid:48) ; ω ) corresponds to the nonlocal plasmon modes. The bare Coulomb potential v ( r ) in this open system wouldbe screened, leading to U eff ( r , ω ) = (cid:82) d r (cid:48) K ( r , r (cid:48) ; ω ) v ( r (cid:48) ). The mean-field formalism for nonlocal plasmon modes ofa 2D layer interacting with a thick conductor was discussed in Refs. [20,69–71]. This theory had given rise to a linearplasmon branch, which was later confirmed by an experiment. Similar plasmon branches were also shown for doublegraphene layers at different temperatures.
Furthermore, our previous work studied the plasmon instability inthe graphene double-layer system, predicting an instability-based terahertz emission, and nonlocal plasmons in ametal-graphene-metal encapsulated structure. As it was reported before the energy gap plays a crucial role on the nonlocal collective modes since it affects bothplasmon branches and Landau damping due to single-particle excitations. Taking this into account, we focus on theeffects of energy gaps and particle-hole modes (PHMs) in all three distinguishable insulating regimes of silicene, i.e.,TI, VSPM and regular BI.
Here, both the energy gaps and PHMs could be independently tuned by applyingan electric field perpendicular to the silicene layer.The rest of the paper is organized as follows. In Sec. II, we first study nonlocal plasmon excitations in the silicenesystem. Our numerical results show plasmon modes and their damping as functions of frequencies and wave numbersfor different parameters including Fermi energies, spin-orbit and sublattice asymmetry bandgaps, and various surface-plasmon frequencies and coupling strengths. In some limiting cases, analytical results for the plasmon modes arepresented and analyzed for their dependencies on sample structure parameters. In Sec. III, we further explore theband-energy dispersions and the electronic states of molybdenum disulfide, indicating relevance to recently discoveredgroup IV dichalcogenides. For these materials, crucial analytical results are obtained for the wave functions, overlapfactor and Fermi energy, which have not yet been addressed adequately in the literature. Additionally, we alsoinvestigate nonlocal plasmon branches in the long-wavelength limit and demonstrate how they are affected by themismatch of n and p doping types and densities. Finally, concluding remarks and discussion of our numerical resultsin this paper are presented in Sec. IV. II. HYBRID PLASMON MODES AND DAMPING IN OPEN SILICENE SYSTEMS
In this section, we discuss hybrid-plasmon dynamics in a silicene open system at low temperatures, and we willaddress the molybdenum-disulfide open system in the next section. For a single silicene layer, the low-energy Hamil-tonian at the corners of the first Brillouin zone isˆ H × = (cid:126) v F ( ξk x ˆ τ x + k y ˆ τ y ) ⊗ ˆ I × − ξ ∆ SO ˆ σ z ⊗ ˆ τ z + ∆ z ˆ τ z ⊗ ˆ I × , (1)where ∆ SO is the intrinsic spin-orbit energy gap, ∆ z ∝ E ⊥ represents the field-dependent sublattice asymmetrybandgap with E ⊥ as an applied electric field perpendicular to the lattice, ˆ I × is the unit matrix, ˆ τ x,y,z and ˆ σ x,y,z are the 2 × ξ = ± K and K (cid:48) valleys, and the silicene Fermi velocity v F is just half of that in graphene.This Hamiltonian in Eq. (1) can be cast into block-diagonal form with two 2 × H ξ,σ = (cid:20) − ξσ ∆ SO + ∆ z (cid:126) v F ( ξk x − ik y ) (cid:126) v F ( ξk x + ik y ) ξσ ∆ SO − ∆ z (cid:21) , (2)where σ = ± σ z and ξ = ± E ξ,σ ( k ), associatedwith Eq. (2) are ± E ξ,σ ( k ) = ± (cid:113) (cid:126) v F k + ∆ ξ,σ , (3)which represent a pair of spin-dependent energy subbands for each valley and have two corresponding non-equivalentbandgaps ∆ ξ,σ ≡ | ∆ SO − ξσ ∆ z | = | ∆ SO ± ∆ z | . The positive (negative) sign in Eq. (3) is for electron (hole) states. Forsimplicity, we, therefore, introduce the notations, ∆ > = ∆ SO + ∆ z and ∆ < = | ∆ SO − ∆ z | , for these two unequivalentbandgaps. The key issues in this section are obtaining spectra of hybrid-plasmon excitations in 2DMOS with newingredients ∆ >,< and calculating screened Coulomb couplings of electrons in silicene layers to an adjacent semi-infinitebulk plasma.It is known that the coupling of electrons in 2D materials to other conduction electrons in 2DMOS will change theplasmon-mode dispersion. Since the damping region is still decided by the PHMs in 2D materials, this will lead toa modification to the damping of plasmon modes in 2DMOS. According to Refs. [20,69–71], the Fourier-transformednonlocal composite inverse dielectric function can be determined by K ( z , z ; q, ω ) = K S ( z , z ; q, ω ) + Π ( q, ω ) K S ( a, z ; q, ω ) S C ( q, ω ) (cid:26)(cid:90) ∞−∞ dz (cid:48) K S ( z , z (cid:48) ; q, ω ) v c ( q, | z (cid:48) − a | ) (cid:27) . (4)In Eq. (4), Π ( q, ω ) is the electron polarizability of silicene (explicitly given below), the interaction of silicene withthe substrate is included in the second term, a is the separation of the silicene layer from the conducting surface,and v c ( q, | z − z (cid:48) | ) = ( e / (cid:15) (cid:15) r ) exp( − q | z − z (cid:48) | ) with (cid:15) r as the average dielectric constant of silicene and spacer layer.Moreover, K S ( z , z ; q, ω ) represents the local inverse dielectric function of the conducting substrate, expressed as K S ( z, z (cid:48) ; q, ω ) = θ ( z ) (cid:26) δ ( | z | − z (cid:48) ) + δ ( z (cid:48) ) e − q | z | (cid:20) − (cid:15) B ( ω )1 + (cid:15) B ( ω ) (cid:21)(cid:27) + θ ( − z ) (cid:26) δ ( | z | + z (cid:48) ) (cid:15) B ( ω ) + δ ( z (cid:48) ) e − q | z | (cid:15) B ( ω ) (cid:20) (cid:15) B ( ω ) − (cid:15) B ( ω ) + 1 (cid:21)(cid:27) , (5)where θ ( z ) is a unit-step function, z > z <
0) corresponds to the spacer layer (conductor), separated by the surfaceat z = 0. We have also assumed a Drude model in Eq. (5) for the substrate dielectric function (cid:15) B ( ω ) = 1 − Ω p /ω ,where Ω p = (cid:112) n e /(cid:15) (cid:15) b m ∗ is the bulk-plasma frequency. Here, Ω p depends on the electron concentration n , substratedielectric constant (cid:15) b , the effective mass m ∗ of electrons, and it can vary in a very large range from ultra-violet (metals)down to infrared or even terahertz (doped semiconductors) frequencies. The use of the Drude model in Eq. (5) canbe justified by a short screening length for high electron concentrations in bulk materials.Finally, we are in a position to calculate the hybrid-plasmon modes in 2DMOS. The plasmon dispersions for asingle silicene layer can be obtained from the dielectric-function equation: ε ( q, ω ) = 1 − (2 πα/q ) Π ( q, ω ) = 0, where α = e / π(cid:15) (cid:15) r . For 2DMOS, on the other hand, ε ( q, ω ) should be replaced by the so-called “dispersion factor” S C ( q, ω ), which appears in Eq. (4) and is calculated as S C ( q, ω ) = 1 − (cid:18) παq (cid:19) Π ( q, ω ) (cid:34) e − qa Ω p ω − Ω p (cid:35) . (6)This verifies that the plasmon dispersions in 2DMOS will indeed be modified by coupling to other conduction electrons(Ω p (cid:54) = 0).Since the Coulomb coupling between electrons in different ( K and K (cid:48) ) valleys involves two uncompensated very largelattice wave numbers, the resulting electron interaction becomes negligible in comparison with those of electrons withinthe same valley. Consequently, by ignoring inter-valley Coulomb coupling and using the one-loop approximation for silicene, we find Π ( q, ω ) = (cid:88) β = >,< Π ( q, ω ; ∆ β ) , (7)and for each subband we getΠ ( q, ω ; ∆ β ) = 14 π (cid:90) d k (cid:88) s,s (cid:48) = ± (cid:34) ss (cid:48) ( (cid:126) v F ) k · ( k + q ) + ∆ β E β ( k ) E β ( | k + q | ) (cid:35) f [ s E β ( k )] − f [ s (cid:48) E β ( | k + q | )] s E β ( k ) − s (cid:48) E β ( | k + q | ) − (cid:126) ω − i + , (8)where s, s (cid:48) = ± f ( E ) = θ ( E − E F ) at zero temperature and E F is theFermi energy of electrons in silicene. A. Approximate analytical results
In the long-wavelength limit q (cid:28) k βF ( k βF is the Fermi wave number for each subband), the bare bubble polarizationfunction for E F > ∆ > is obtained as Π ( q, ω ) = 1 π (cid:88) β = >,< k βF (cid:12)(cid:12)(cid:12)(cid:12) ∂ E β ( k ) ∂k (cid:12)(cid:12)(cid:12)(cid:12) k = k βF q (cid:126) ω = E F π (cid:18) − ∆ < E F − ∆ > E F (cid:19) q (cid:126) ω , (9)where E F = (cid:113) ( (cid:126) v F k βF ) + ∆ β , k βF = (cid:112) πρ β , and ρ β is the electron areal density for each subband. Therefore, from ε ( q, ω ) = 0 for a single silicene layer, we get the following plasmon branch ω p ( q ) = 4 α (cid:126) E F (cid:18) E F − ∆ > + ∆ < (cid:19) q ≡ G q , (10)where, for convenience, we introduce a coefficient G ≡ G ( E F , ∆ β ).It is important to note that the Fermi energy E F for silicene is fixed by the total electron areal density ρ through E F − (cid:0) ∆ > + ∆ < (cid:1) = ( (cid:126) v F ) π ( ρ > + ρ < ) ≡ ( (cid:126) v F ) πρ , (11)and if ρ is small, we can further approximately obtain E F (cid:119) (cid:114) ∆ > + ∆ < (cid:126) v F ) πρ (cid:112) > + ∆ < ) . (12)This leads to G (cid:119) √ αv F πρ / ( (cid:112) ∆ > + ∆ < ). For gapped graphene, we have ∆ < = ∆ > = ∆, and Eq. (12) gives rise to E F − ∆ (cid:119) ( (cid:126) v F ) πρ / (2∆). In this case, we get G (cid:119) αv F πρ / ∆ and from Eq. (10) we find ω p ( q ) ∼ √ ρ q . Actually,such a scaling relation holds true for all 2D materials except for gapless graphene which yields ω p ( q ) ∼ (cid:113) ρ / q .Additionally, in contrast to Eq. (10), if E F < ∆ > with an unoccupied upper subband, we obtain the plasmon mode ω p ( q ) = 2 αE F (cid:126) (cid:18) − ∆ < E F (cid:19) q . (13)In the above discussion, we are only restricted to the plasmon mode for a stand-alone silicene. For the 2DMOS, onthe other hand, when both subbands are occupied, the plasmon modes are determined by Eqs. (6) and (9), yielding1 − G qω (cid:32) e − qa Ω p ω − Ω p (cid:33) = 0 , (14)which gives rise to the following bi-quadratic equation2 (cid:18) ω Ω p (cid:19) − (cid:18) G q Ω p (cid:19) (cid:18) ω Ω p (cid:19) + G q Ω p (cid:0) − e − qa (cid:1) = 0 . (15)Its two solutions are simply given by4 ω p, ± Ω p = (cid:18) G q Ω p (cid:19) ± (cid:34)(cid:18) G q Ω p (cid:19) − G q Ω p (cid:0) − e − qa (cid:1)(cid:35) / , (16)where the sign + ( − ) corresponds to the in-phase (out-of-phase) plasmon mode. For the strong-coupling regime with qa (cid:28) ω p, + ( q ) (cid:119) Ω p √ G q √ p − G + 4 a Ω p √ p G q + O ( q ) ,ω p, − ( q ) (cid:119) q √ a G − √ a G Ω p G q + O ( q ) . (17)In Eq. (17), the linear dispersions and their prefactor scalings are the same as those for graphene. However, thetwo independent bandgaps, ∆ SO and ∆ z , play unique roles in shaping the hybrid-plasmon branches when the dampingfrom different PHMs is considered. It is found that the outer PHM’s boundaries are only determined by ∆ < and thetwo hybrid-plasmon group velocities (slopes) are proportional to G and √ G . These two group velocities drop to zeroas ∆ SO and ∆ z increase since G (cid:119) √ αv F πρ / (cid:112) ∆ > + ∆ < for low doping. If a proper applied electric field is chosen,the VSPM phase can be reached with ∆ < = 0. For ∆ < >
0, on the other hand, we know the plasmon frequencies forboth gapped graphene and silicene are reduced by finite ∆ < . Meanwhile, these plasmon branches will enter intoa gap region between the interband and intraband PHMs. As a result, we find that both the damping-free plasmonregions and the plasmon group velocities can be controlled independently by ∆ SO and ∆ z . This requirement can befulfilled by scanning an external electric field, even for a fixed spin-orbit interaction strength, leading to distinctivebehaviors for TI and BI phases.Moreover, the plasmon group velocity associated with the ω p, − ( q ) mode in 2DMOS depends on √ a in the strong-coupling regime. However, in the weak-coupling regime with qa (cid:29)
1, this plasmon mode becomes proportional to √ q ,as shown in Eq. (10) for a single silicene layer. Meanwhile, the plasmon group velocity for the ω p, + ( q ) mode, which isindependent of a , approaches zero in the weak-coupling regime. B. Full numerical solutions
For our numerical results presented in Figs. 2-6, we use the scale E for the energy and the scale k for the wavenumber q , where E = (cid:126) v F k , k = √ πρ with ρ as the total doped electron areal density. Here, the constant value ρ = 10 cm − is given for these five figures. ( ) a ( ) b ( ) d ( ) c p l a s m o n p l as m on p l as m on p l a s m o n FIG. 2: Numerical results for the hybrid-plasmon branches in 2DMOS with ∆ SO /E = 0 . z /E = 0 . < /E = 0 . a )-( c ) correspond to various cases with k a = 1 .
0, 3 . .
0. Here, the plasma energy (cid:126) Ω p /E = 1 .
0. In all panels, theblue solid curves are obtained from | S c ( q, ω ) | = 0, while the red short-dashed curves are from Re [ S c ( q, ω )] = 0, demonstratingboth undamped and damped plasmon branches. The PHM regions are depicted by partially transparent green areas enclosedby dashed-curve boundaries. Panel ( d ) gives the density plot for the energy loss function Im [1 /(cid:15) ( q, ω )] of free-standing silicene.The populations of two subbands are shown in the inset. ( ) a ( ) b ( ) c ( ) d p l a s m o n p l a s m o n p l a s m o n p l a s m o n FIG. 3: Numerical results for the hybrid-plasmon branches in 2DMOS with ∆ SO /E = 0 . z /E = 0 . < /E = 0 . a )-( c ) correspond to k a = 1 .
0, 3 . .
0, respectively. Here, (cid:126) Ω p /E = 1 .
0. In all panels, the blue solid curves arefor the undamped plasmon modes, while the red short-dashed curves are for the damped plasmon modes. The PHM regionsare depicted by partially transparent green areas enclosed by dashed-curve boundaries. Panel ( d ) gives the density plot forIm [1 /(cid:15) ( q, ω )] of a single silicene layer. The populations of two subbands are shown in the inset. ( ) a ( ) b ( ) d p l a s m o n p l a s m o n p l a s m o n ( ) c p l a s m o n FIG. 4: Numerical results for the hybrid-plasmon branches in 2DMOS with ∆ SO /E = 0 . z /E = 0 . < /E = 0 . a )-( c ) correspond to k F a = 1 .
0, 3 . .
0, respectively. Here, (cid:126) Ω p /E = 1 .
0. In all panels, the blue solid curves arefor the undamped plasmon modes, while the red short-dashed curves for the undamped plasmon modes. The PHM regionsare depicted by partially transparent green areas enclosed by dashed-curve boundaries. Panel ( d ) gives the density plot forIm [1 /(cid:15) ( q, ω )] of a single silicene layer. The populations of two subbands are shown in the inset. The features of the hybrid plasmon modes beyond the long-wavelength limit could be explored numerically forall possible values of the energy bandgaps based on the exact calculation of the polarization function for a silicenelayer and the use of Eq. (6). Here, two subbands can be selectively populated by controlling E F or ρ . Furthermore,the coupling of electrons to the surface of a semi-infinite conductor in 2DMOS can also be tuned by choosing theseparation of the silicene layer from the bulk surface.We first consider a case with a relatively large minimal bandgap ∆ < /E F = 0 . a . The anticrossing of two hybrid plasmon modes canbe seen most clearly in Fig. 2( b ) with k F a = 3 .
0. In comparison with the plasmon damping for free-standing slicenein Fig. 2( d ), the upper optical-like branch is free from damping into the main diagonal ( ω = v F q , intraband PHM)until exceeding a relatively large critical wave number, as shown in Figs. 2( a )-2( c ). For ∆ < /E F = 0 . k F a = 5 . c ), the anticrossing feature becomes almost indistinguishable, and the lower branchapproaches that of free-standing slicene in Fig. 3( d ).The situation with an even smaller bandgap ∆ < /E F = 0 . d ), however, the damping always occurs at one of thePHM boundaries, and the increase of bandgap makes it more favorable for the plasmon-mode damping to occur atthe interband PHM region. In addition, we also find that the anticrossing feature becomes more significant in theweak-coupling regime, as displayed in Fig. 4( c ).It is known that the external parameter in the 2DMOS, i.e., the plasma frequency Ω p , can greatly affect the dampingof the hybrid plasmon modes. Our results for different values of Ω p are shown in Fig. 5. It is very surprising to seefrom Figs. 5( a )-5( e ) that the damping-free range of the lower branch will depend on Ω p but not on the other internalparameters, such as a , ∆ < and ρ . From Fig. 5( f ), on the other hand, we observe that the upper branch could bedoubly damped by both intraband and interband PHM regions, which also exists for a gapped graphene open system.It is reasonable to expect that the open-system damping effect will become more significant if the conductor surfacestays closer to the silicene layer. The numerical results for the plasmon dispersions with much smaller separations FIG. 5: Numerical results for the hybrid-plasmon branches in 2DMOS with ∆ SO /E = 0 . z /E = 0 .
15 (∆ < /E = 0 . a )-( f ) correspond to (cid:126) Ω p /E = 0 .
5, 0 .
8, 1 .
0, 1 .
2, 1 . .
5, respectively. Here, k a = 1 .
0. In all panels, the blue solidcurves are for the undamped plasmon modes, while the red short-dashed curves for the undamped plasmon modes. The PHMregions are depicted by partially transparent green areas enclosed by with dashed-curve boundaries. The populations of twosubbands are shown in the inset. a are presented in Fig. 6, from which we reproduce a recent experimentally confirmed effect in graphene, i.e.,the acoustic-like plasmon branch will be highly damped in the long-wavelength limit as a < . a )-6( f ) with various energy gaps. It is interesting to note fromFig. 6 that the group velocity of the upper branch almost does not depend on the separation a , but strongly dependson the energy gap ∆ < . III. HYBRID PLASMON MODES AND DAMPING IN OPEN MOLYBDENUM-DISULFIDE SYSTEMS
The two-band model Hamiltonian for molybdenum disulfide, as well as for most other transition-metal dichalco-genides, next to the two inequivalent K and K (cid:48) valley points can be written as ˆ H τ,sd = (cid:18) τ s λ + (cid:126) k m e α (cid:19) ˆ I × + (cid:18) ∆2 − τ s λ + (cid:126) k m e β (cid:19) ˆ σ z + t a ˆ Σ τ · k (18)where τ = ± s = ± . λ = 0 .
042 ∆ is thespin-orbit coupling parameter, m e represents the free electron mass, ˆ Σ τ are the Pauli matrices for valley pseudospins, t = 0 .
884 ∆ is the electron hopping parameter, and a = 1 .
843 ˚A which is obtained from the Mo − S atom-atom bondlength 2 .
43 ˚A. Even though λ (cid:28) ∆, the spin-orbit interaction is not negligible, which is reflected in the spin-resolved0 FIG. 6: Numerical results for the hybrid-plasmon branches in 2DMOS. Here, plots ( a ), ( c ) and ( e ) are for k a = 0 .
2, whileplots ( b ), ( d ) and ( f ) for k a = 0 .
7. In addition, we assume ∆ SO /E = 0 . z /E = 0 . < /E = 0 .
5) in ( a ) and ( b ),∆ SO /E = 0 . z / F = 0 . < /E = 0 .
3) in ( c ) and ( d ), and ∆ SO /E = 0 . z /E = 0 . < /E = 0 .
2) in ( e )and ( f ). In all panels, the blue solid curves are for the undamped plasmon modes, while the red short-dashed curves for theundamped plasmon modes. The PHM regions are depicted by partially transparent green areas enclosed by the dashed-curveboundaries. The populations of two subbands are shown in different insets. energy subbands and in the absence of spin degeneracy. In Eq. (18), we use α = 2 .
21 = 5 . β and we find that t a = 4 . × − J · m plays a role of the Fermi velocity and is equal to 0 .
472 of the (cid:126) v F factor for graphene.Moreover, we neglect the trigonal warping term t a ( ˆ Σ τ · k ) ˆ σ x ( ˆ Σ τ · k ), which leads to the slight anisotropy of theenergy for our whole study since t = 0 . .
053 ∆ does not represent a considerable effect on the electronic states.Consequently, the considered hybrid plasmon are also isotropic.It is easy to verify that the Hamiltonian in Eq. (18) is equivalent to that of gapped graphene with a k − dependent“gap” term, ∆ τ,s ( k ) = ∆ / − τ s λ / (cid:126) k β/ (4 m e ), as well as a k − dependent band-shift term, E τ,s ( k ) = τ s λ / (cid:126) k α/ (4 m e ), yielding ε τ,sγ ( k ) = E τ,s ( k ) + γ (cid:113) [∆ τ,s ( k )] + ( t a k ) , (19)where γ = ± O ( k ) for small k values, Eq. (19) turns into ε τ,sγ ( k ) (cid:119) τ sλ + α (cid:126) m e k + γ (cid:113) (∆ − τ s λ ) + [(2 t a ) + (∆ − τ s λ ) β (cid:126) /m e ] k . (20)1Besides the simple plane-wave part, the spinor parts of the wave functions associated with the eigenvalues in Eq. (20)for each valley are given by Ψ τ,sγ ( k ) = 1 (cid:112) δε τ,sγ ( k ) /γ (cid:112) | δε τ,sγ ( k ) + ∆ τ,s ( k ) | γ (cid:112) | δε τ,sγ ( k ) − ∆ τ,s ( k ) | e iθ k , (21)where θ k = tan − ( k y /k x ), δε τ,sγ ( k ) ≡ ε τ,sγ ( k ) − E τ,s ( k ), and the overlap factor is calculated as F τ,sγ,γ (cid:48) ( k, k + q ) ≡ |(cid:104) Ψ τ,sγ ( k ) | Ψ τ,sγ (cid:48) ( | k + q | ) (cid:105)| = 12 (cid:20) γ γ (cid:48) ∆ τ,s ( k )∆ τ,s ( | k + q | ) + k · ( k + q ) | δε τ,sγ ( k ) | | δε τ,sγ ( | k + q | ) | (cid:21) . (22)Here, each part of the expression in Eq. (22) could be calculated explicitly, e.g.,∆ τ,s ( k ) ∆ τ,s ( | k + q | ) = (∆ − τ s λ ) (cid:126) β (∆ − τ s λ )8 m e [2 k · ( k + q ) + q ] + (cid:18) (cid:126) β m e (cid:19) k | k + q | . (23)If we introduce the notation for a composite index µ ≡ τ s = ± α and β terms in Eq. (20),this gives rise to ε µγ ( k ) (cid:119) µλ / γ (cid:112) ( t a ) k + (∆ − µλ ) /
4, and therefore, the wave function in Eq. (21) could besimplified asΨ µγ ( k ) = 1 (cid:112) [2 ε µγ ( k ) − µλ ] /γ (cid:112) | ε µγ ( k ) − µλ + ∆ / | γ (cid:112) | ε µγ ( k ) − ∆ / | e iθ k = 1 √ Q k (cid:112) | Q k + γ (∆ − µλ ) / | γ (cid:112) | Q k − γ (∆ − µλ ) / | e iθ k , (24)where Q k = [ ε µγ ( k ) − µλ / /γ = (cid:112) ( t a ) k + (∆ − µλ ) /
4. It is clear from Eq. (24) that, unlike gapped graphene,two spinor components become inequivalent and their ratio is different for electron and hole states.Furthermore, the density of states ρ d ( E ) can be formally written as ρ d ( E ) = 2 (cid:90) d k (2 π ) (cid:88) γ = ± (cid:88) µ = ± δ (cid:2) E − ε µγ ( k ) (cid:3) . (25)By denoting ˘ (cid:15) µ = µλ /
2, ˘∆ µ = (∆ − µλ ) /
2, ˘ A µ = (∆ − µλ ) (cid:126) β/ (4 m e ) + ( t a ) , and ˘ α = (cid:126) α/ (4 m e ), Eq. (20) isfurther simplified to ε µγ ( k ) (cid:119) ˘ (cid:15) µ + ˘ α k + γ (cid:113) ˘∆ µ + ˘ A µ k . Therefore, if ˘ α (cid:54) = 0, from Eq. (25) we obtain the analyticalresult for ρ d ( E ), given by ρ d ( E ) = 12 π (cid:88) ± (cid:88) γ, µ = ± (cid:12)(cid:12)(cid:12) ˘ α + γ ˘ A µ E − ˘ (cid:15) µ − χ ± µ ( E )] (cid:12)(cid:12)(cid:12) − Θ (cid:20) γ (cid:18) E − µλ (cid:19) −
12 (∆ − µλ ) (cid:21) , (26)where Θ( x ) is a unit-step function, and the energy-dependent function χ ± µ ( E ) is defined as χ ± µ ( E ) = 12˘ α (cid:20) ˘ A µ + 2˘ α ( E − ˘ (cid:15) µ ) ± (cid:113) ˘ A µ + 4˘ α ˘∆ µ + 4 ˘ A µ ˘ α ( E − ˘ (cid:15) µ ) (cid:21) . (27)If ˘ α = 0, on the other hand, we simply find ρ d ( E ) = 2 π (cid:88) γ, µ = ± | E − µλ / || ˘ A µ | Θ (cid:20) γ (cid:18) E − µλ (cid:19) −
12 (∆ − µλ ) (cid:21) . (28)2Using the result in Eq. (26), all previously known cases, including a pair of parabolic bands or Dirac cones, as well asa pair of gapped Dirac cones, could be easily verified.If we neglect the α and β terms in Eq. (20), i.e., setting ˘ A µ = ( t a ) and ˘ α = 0, we obtain from Eq. (28) the resultfor a pair of non-degenerate, spin- and valley-dependent subbands in gapped graphene, given by ρ d ( E ) = 2 π ( t a ) (cid:88) γ, µ = ± (cid:12)(cid:12)(cid:12) E − µλ (cid:12)(cid:12)(cid:12) Θ (cid:20) γ (cid:18) E − µλ (cid:19) −
12 (∆ − µλ ) (cid:21) . (29)It is clear from Eqs. (26), (28) and (29) that the boundaries for non-zero density of states in all three cases are set by E > ∆ / γ = +1) and E < − ∆ / µλ for holes ( γ = − t (cid:28) ∆, using Eq. (20) we arrive at ε µγ ( k ) = 12 [ µλ (1 − γ ) + γ ∆] + (cid:20) (cid:126) m e ( α + γβ ) + γ ( t a ) ∆ − µλ (cid:21) k , (30)where we have used the fact that ∆ (cid:29) λ . This result leads to the density of states given by ρ d ( E ) = 12 π (cid:126) (cid:88) γ, µ = ± (cid:12)(cid:12)(cid:12) α + γβ m e + γ ( t a ) (cid:126) (∆ − µλ ) (cid:12)(cid:12)(cid:12) − Θ (cid:20) γ (cid:18) E − µλ (cid:19) −
12 (∆ − µλ ) (cid:21) , (31)where there exist two energy-independent giant discontinuities for electrons and holes, respectively.For electrons with γ = +1 at E = ∆ /
2, we get the jump in the density of states given by δρ γ =+1 d = 12 π (cid:88) µ = ± (cid:20) ( t a ) ∆ − µλ + ( α + β ) (cid:126) m e (cid:21) − = 0 . t a . (32)Similarly, for holes with γ = − E = − ∆ / µλ , we obtain two discontinuities at different energies, i.e., δρ γ = − d ( − ∆ / µλ ) = 12 π (cid:126) α − β m e − ( t a ) (cid:126) × (∆ − λ ) − × (∆ + λ ) − − = 1 t a × . , for µ = +1 × . , for µ = − . (33)Our above analytical expressions match exactly the numerical results reported in Ref. [8]. We note that the simplifiedparabolic dispersions in Eq. (30) catch the difference in electron and hole effective masses for various spin and valleyindices even for k = ⇒
0. For large k values, this difference becomes more significant. The numerically-calculatedenergy dispersions for electrons and holes and their zero-temperature Fermi energies for fixed doping density arepresented in Fig. 7. It is interesting to note that at a finite energy away from the bandedge, the density of states forMoS is significantly smaller compared to graphene. This agrees with the well-known fact that non-parabolicity ingraphene energy subbands will enhance the density of states.Now, let us turn to studies of nonlocal plasmon dispersions and their damping in a MoS layer interacting with asemi-infinite conductor. For this case, the one-loop electron polarization function for molybdenum disulfide with twopairs of energy subbands can be obtained in a way similar to Eq. (7) by summing over a composite index µ for eachsubband but specifying γ values for n − and p − doping separately. We will further assume that only the lowest holesubband with µ = +1 will be occupied, and two electron subbands become nearly degenerate with each other since∆ (cid:29) λ .In the long-wavelength limit, by using Eq. (9) the polarization function of a MoS monolayer interacting with asemi-infinite conductor separated by a distance D can be expressed as ω γp, − ( q ) = q (cid:112) D L ( γ ) (cid:119) q (cid:112) πα ρ D ( γ + 3) 2( t a ) / (cid:126) (cid:112) ∆ − (1 − γ ) λ ,ω γp, + ( q ) = Ω p √ L ( γ ) q √ p (cid:119) Ω p √ γ + 3 √ p ( t a ) πα ρ (cid:126) [∆ − (1 − γ ) λ ] q , (34)3 ( ) a ( ) i . . -2.0 -1.0 0.0 2.01.88-0.63-1.25-1.88 ( ) b - ( ) i FIG. 7: Electron and hole energy dispersions (solid curves) and their Fermi energies for a monolayer MoS . Plot ( a ) representsthe dispersions of energy subbands near the K -valley ( τ = 1), where the corresponding results in the parabolic approximationare also shown by dashed curves. In addition, the Dirac-cone dispersions ε γ ( k ) = γ t a | k | are included. The inset ( i
1) gives aclose-look view for two very close conduction bands with opposite spins s = ±
1. Plot ( b ) presents the calculated Fermi energies E F /t for electrons ( γ = +1) [holes ( γ = − a πρ in MoS (red solid curvewith neglected k -terms in the Hamiltonian), gapped graphene (blue dashed curve) and Dirac cones (green dash-dotted curve)outside the gap region. The corresponding density-of-states curves for these materials are displayed in the inset ( i where ρ is the areal density for doping, α = e / (4 π(cid:15) (cid:15) r t a ) (cid:119) . (cid:15) r (cid:119) , and L ( γ ) = 2 π ( t a ) α ρ ( γ + 3) (cid:126) [∆ + (1 − γ ) λ ] . (35)Here, the inclusion of the coupling between MoS and the semi-infinite conductor has split plasmons into in-phase(+) and out-of-phase ( − ) modes in Eq. (34). This will certainly lead to a modification of plasmon-mode damping byPHMs.The main advantage for using MoS in a hybrid plasmonic device is its large energy gap (cid:119) ∆, which allows oneto consider clean metals with an extremely high plasma frequency (cid:126) Ω p (cid:119) eV . This arrangement is not possible forgapped graphene or silicene since the plasmon modes at such a frequency would be strongly damped by the interbandPHMs. Another unique feature for MoS is the large difference between the electron and hole doping processes, i.e.,high doping density ρ (cid:119) − cm − only allows the occupation of one hole subband, as assumed in Eq.(34) for γ = −
1. Here, even in the parabolic approximation, the results for n − and p − doping still vary drastically. Althoughthe λ correction to ∆ is very small, the density of states of electrons is almost twice as large as that of holes.In order to determine Landau damping of the plasmon modes, we need to determine the boundaries (cid:126) Ω γc ( q ) forPHMs, defined by (cid:126) ω γp, ± ( q ) ≥ (cid:126) Ω γc ( q ) ≡ ε µγ ( k F + q ) − ε µγ ( k F ) , (36)which corresponds to k (cid:107) q . Here, k F is the electron Fermi wave number. For moderate n − doping ( γ = +1), fromEq. (20) its PHM boundary is found to be (cid:126) Ω γ =+1 c ( q ) = λ (cid:126) α m e ( q + k F ) − E F + 12 (cid:115) ∆ + (cid:18) (cid:126) β ∆ m e + 4 t a (cid:19) ( q + k F ) . (37)In the long-wavelength limit with q (cid:28) k F , we can approximate Eq. (37) by (cid:126) Ω γ =+1 c ( q ) (cid:119) q (cid:115)(cid:20) t a + (cid:126) ∆ m e ( β + α ) (cid:21) ( E F − ∆ / , (38)where we use the facts that ∆ (cid:29) λ , t a k F and (cid:126) βk F /m e and E F is determined by ρ .Alternatively, if the sample is p − doped ( γ = − (cid:126) Ω γ = − c ( q ) with µ = +1 for the occupiedhole subband is found to be4 (cid:126) Ω γ = − c ( q ) (cid:119) q (cid:115)(cid:20) t a + (cid:126) (∆ − λ ) m e ( β − α ) (cid:21) ( E F − ∆ / λ )(∆ − λ ) , (39)These two PHM boundaries, (cid:126) Ω γ = ± c ( q ), determine whether the acoustic-like plasmon branch would be Landaudamped or not. On the other hand, the optical-like plasmon branch originating from (cid:126) Ω p / √ (cid:119) . IV. SUMMARY AND CONCLUDING REMARKS
In conclusion, we have presented in this paper the numerical results for full ranges of hybrid plasmon-mode dis-persions, as well as analytical expressions in the long-wavelength limit, in an open interacting system including a 2Dmaterial and a conducting substrate. Although the plasmon damping is set by the particle-hole modes(PHMs) ofelectrons in the 2D material, the strong coupling between electrons in 2D materials and in the conducting substrategives rise to a splitting of plasmons into one in-phase and one out-of-phase mode. Such dramatic changes in plasmondispersions are expected to have impacts on the damping of these modes. In addition, in comparison with gappedgraphene, the different plasmon modes in silicene or transition-metal dichalcogenides make our damping studies evenmore distinctive, including different energy bandgaps, doping types, occupations of subbands, and coupling between2D materials and the conducting substrate. Here, each plasmon branch and its damping can be independently an-alyzed based on the signatures of the PHMs since the plasmon modes depend on both spin and valley degrees offreedom. Therefore, our proposed hybrid systems in this paper are expected to be useful in measuring the dielectricproperty of 2D material open systems (2DMOS) and spin-orbit coupling strength of individual 2D materials. Moreimportantly, we have demonstrated the possibility to design the plasmonic resonances at almost all frequencies andwave numbers for different types of newly discovered 2D materials. This was not feasible for either a free-standingsilicene layer or a graphene-based hybrid structure.Additionally, our model and numerical results for 2DMOS have confirmed a recently discovered phenomenon relatedto a significant damping of an acoustic-like plasmon branch as the separation to the conducting substrate becomesvery small. From our current studies, we have found that in silicene this critical distance can be modified by eitherapplying an external electric field or varying doping types and levels. The unique linear dispersion obtained under thelong-wavelength limit makes the damping from intraband PHMs possible in 2DMOS but not for a free-standing 2Dlayer. We have also noted that the plasma energies in clean metals are usually much larger than the Fermi energiesand bandgaps in graphene. As a result, the plasmon modes in graphene can not be coupled to surface plasmons inthe presence of a metallic substrate without suffering from the strong damping by interband PHMs. However, theuse of MoS with a large bandgap in 2DMOS is able to suppress this damping effectively. Alternatively, one couldalso use Bi Se material, which is a doped topological insulator with a surface-plasmon energy around 104 meV, ora highly-doped semiconductor in 2DMOS. Acknowledgments
D.H. would like to thank the support from the Air Force Office of Scientific Research (AFOSR). We would like tomention a great help from Joseph Sadler and especially Patrick Helles, IT managers at the Center for High TechnologyMaterials of the University of New Mexico, much beyond their official responsibilities. ∗ Electronic address: [email protected] F. Stern, Phys. Rev. Lett. , 546 (1967). E. H. Hwang and S. Das Sarma, Phys. Rev. B , 205418 (2007). B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics , 318 (2006). P. Pyatkovskiy, Journal of Physics: Condensed Matter , 025506 (2008). T. Stauber, J. Schliemann, and N. M. R. Peres, Phys. Rev. B , 085409 (2010). A. Scholz, T. Stauber, and J. Schliemann, Phys. Rev. B , 195424 (2012). S. Das Sarma and Q. Li, Phys. Rev. B , 235418 (2013). A. Scholz, T. Stauber, and J. Schliemann, Phys. Rev. B , 035135 (2013). K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos, and A. Firsov, Nature , 197(2005). A. K. Geim and K. S. Novoselov, Nature Materials , 183 (2007). A. C. Neto, F. Guinea, N. Peres, K. S. Novoselov, and A. K. Geim, Reviews of Modern Physics , 109 (2009). A. Politano and G. Chiarello, Nanoscale , 10927 (2014). V. Despoja, D. Novko, K. Dekani´c, M. ˇSunji´c, and L. Maruˇsi´c, Phys. Rev. B , 075447 (2013). N. Papasimakis, Z. Luo, Z. X. Shen, F. D. Angelis, E. D. Fabrizio, A. E. Nikolaenko, and N. I. Zheludev, Opt. Express ,8353 (2010). V. G. Kravets, F. Schedin, R. Jalil, L. Britnell, K. S. Novoselov, and A. N. Grigorenko, Journal of Physical Chemistry C , 3882 (2012). F. Xia, V. Perebeinos, Y.-m. Lin, Y. Wu, and P. Avouris, Nature nanotechnology , 179 (2011). J. Yan, K. S. Thygesen, and K. W. Jacobsen, Physical review letters , 146803 (2011). F. H. L. F. H. L. Koppens, T. Mueller, P. Avouris, A. C. Ferrari, M. S. Vitiello, and M. Polini, Nature Nanotechnology ,780 (2014). W. Han, R. K. Kawakami, M. Gmitra, and F. J., Nature Nanotechnology , 794 (2014). G. Gumbs and D. Huang,
Properties of Interacting Low-Dimensional Systems (John Wiley & Sons, 2013). M. Campisi, P. H¨anggi, and P. Talkner, Rev. Mod. Phys. , 771 (2011). M. Campisi, P. Talkner, and P. H¨anggi, Phys. Rev. Lett. , 210401 (2009). M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod. Phys. , 1665 (2009). G. E. Crooks, Journal of Statistical Mechanics: Theory and Experiment , 10023 (2008). S. Mukamel, Phys. Rev. Lett. , 170604 (2003). W. De Roeck and C. Maes, Phys. Rev. E , 026115 (2004). G. E. Crooks, Phys. Rev. E , 2721 (1999). F. Setiawan and S. D. Sarma, arXiv preprint arXiv:1509.05067 (2015). U. Weiss,
Quantum dissipative systems , vol. 10 (World Scientific, 1999). E. Illes, C. Roy, and S. Hughes, Optica , 689 (2015). M. Silaev, T. T. Heikkil¨a, and P. Virtanen, Phys. Rev. E , 022103 (2014). G. Schaller,
Open Quantum Systems Far from Equilibrium (Springer, Lecture Notes in Physics, 2014). M. Katsnelson, K. Novoselov, and A. Geim, Nature Physics , 620 (2006). M. Ezawa, New Journal of Physics , 033003 (2012). G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, Physical Review B , 073103 (2007). N. Kharche and S. K. Nayak, Nano Letters , 5274 (2011). Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng, and Z. X. Shen, ACS Nano , 2301 (2008). O. Kibis, Physical Review B , 165433 (2010). C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801 (2005). M. Ezawa, Phys. Rev. Lett. , 055502 (2012). C.-C. Liu, W. Feng, and Y. Yao, Phys. Rev. Lett. , 076802 (2011). L. Zhang, P. Bampoulis, A. van Houselt, and H. Zandvliet, Applied Physics Letters , 111605 (2015). A. Acun, L. Zhang, P. Bampoulis, M. Farmanbar, A. van Houselt, A. Rudenko, M. Lingenfelder, G. Brocks, B. Poelsema,M. Katsnelson, et al., Journal of Physics: Condensed Matter , 443002 (2015). L. Li, S.-z. Lu, J. Pan, Z. Qin, Y.-q. Wang, Y. Wang, G.-y. Cao, S. Du, and H.-J. Gao, Advanced Materials , 4820 (2014). M. D´avila, L. Xian, S. Cahangirov, A. Rubio, and G. Le Lay, New Journal of Physics , 095002 (2014). P. Bampoulis, L. Zhang, A. Safaei, R. Van Gastel, B. Poelsema, and H. J. W. Zandvliet, Journal of Physics: Condensedmatter , 442001 (2014). M. Derivaz, D. Dentel, R. Stephan, M.-C. Hanf, A. Mehdaoui, P. Sonnet, and C. Pirri, Nano Letters , 2510 (2015). A. Iurov, G. Gumbs, and D. H. Huang, Journal of Physics: Condensed Matter ((to appear)). F. dAcapito, S. Torrengo, E. Xenogiannopoulou, P. Tsipas, J. M. Velasco, D. Tsoutsou, and A. Dimoulas, Journal of Physics:Condensed Matter , 045002 (2016). C. J. Walhout, A. Acun, L. Zhang, M. Ezawa, and H. J. W. Zandvliet, Journal of Physics: Condensed Matter , 284006(2016). B. Aufray, A. Kara, S. Vizzini, H. Oughaddou, C. Leandri, B. Ealet, and G. Le Lay, Applied Physics Letters , 183102(2010). P. De Padova, C. Quaresima, C. Ottaviani, P. M. Sheverdyaeva, P. Moras, C. Carbone, D. Topwal, B. Olivieri, A. Kara,H. Oughaddou, et al., Applied Physics Letters , 261905 (2010). B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet, and B. Aufray, Applied Physics Letters , 223109(2010). C. J. Tabert and E. J. Nicol, Physical Review Letters , 197402 (2013). C. J. Tabert and E. J. Nicol, Physical Review B , 235426 (2013). C. J. Tabert and E. J. Nicol, Physical Review B , 085434 (2013). C. J. Tabert and E. J. Nicol, Phys. Rev. B , 195410 (2014). X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057 (2011). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045 (2010). J.-A. Yan, M. A. D. Cruz, S. Barraza-Lopez, and L. Yang, Applied Physics Letters , 183107 (2015). K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. , 136805 (2010). B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nature Nanotechnology , 147150 (2010). T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, and E. W. B. L. . J. F. Qian Niu, Pingheng Tan, Nature Communications , 887 (2012). H. Rostami, A. G. Moghaddam, and R. Asgari, Phys. Rev. B , 085440 (2013). D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. , 196802 (2012). O. L. Berman and R. Y. Kezerashvili, Phys. Rev. B , 245410 (2016). L. Zhang, P. Bampoulis, A. N. Rudenko, Q. Yao, A. van Houselt, B. Poelsema, M. I. Katsnelson, and H. J. W. Zandvliet,Physical Review Letters , 256804 (2016). G. Gumbs, A. Iurov, and N. J. M. Horing, Phys. Rev. B , 235416 (2015). N. J. Morgenstern Horing, E. Kamen, and H.-L. Cui, Phys. Rev. B , 2184 (1985). K. A. Kouzakov and J. Berakdar, Phys. Rev. A , 022901 (2012). N. J. M. Horing, Phys. Rev. B , 193401 (2009). C. Kramberger, R. Hambach, C. Giorgetti, M. H. R¨ummeli, M. Knupfer, J. Fink, B. B¨uchner, L. Reining, E. Einarsson,S. Maruyama, et al., Phys. Rev. Lett. , 196803 (2008). G. Gumbs, A. Iurov, and D. Huang, Coherent Phenomena , 1 (2014). A. Iurov, G. Gumbs, D. Huang, and V. Silkin, Physical Review B , 035404 (2016). N. J. Horing, A. Iurov, G. Gumbs, A. Politano, and G. Chiarello, in
Low-Dimensional and Nanostructured Materials andDevices (Springer, 2016), pp. 205–237. G. Gumbs, A. Iurov, J.-Y. Wu, M. Lin, and P. Fekete, Scientific Reports (2016). G. Gumbs, A. Iurov, D. Huang, and W. Pan, Journal of Applied Physics , 054303 (2015). G. Gumbs, N. Horing, A. Iurov, and D. Dahal, Journal of Physics D: Applied Physics , 225101 (2016). M. Ezawa, Phys. Rev. B , 161407 (2012). S. M. Badalyan, A. Matos-Abiague, G. Vignale, and J. Fabian, Phys. Rev. B , 205305 (2009). R. Sensarma, E. H. Hwang, and S. Das Sarma, Phys. Rev. B , 195428 (2010). A. Politano, A. R. Marino, V. Formoso, D. Far´ıas, R. Miranda, and G. Chiarello, Phys. Rev. B , 033401 (2011). A. Politano, A. R. Marino, and G. Chiarello, Phys. Rev. B , 085420 (2012). A. Politano and G. Chiarello, Applied Physics Letters , 201608 (2013). A. Politano, V. M. Silkin, I. A. Nechaev, M. S. Vitiello, L. Viti, Z. S. Aliev, M. B. Babanly, G. Chiarello, P. M. Echenique,and E. V. Chulkov, Phys. Rev. Lett.115