First-principles study of ultrafast dynamics of Dirac plasmon in graphene

Exploring low-loss two-dimensional plasmon modes is considered central for achieving light manipulation at the nanoscale and applications in plasmonic science and technology. In this context, pump–probe spectroscopy is a powerful tool for investigating these collective modes and the corresponding energy transfer processes. Here, I present a first-principles study on non-equilibrium Dirac plasmon in graphene, wherein damping channels under ultrafast conditions are still not fully explored. The laser-induced blueshift of plasmon energy is explained in terms of thermal increase of the electron–hole pair concentration in the intraband channel. Interestingly, while damping pathways of the equilibrium graphene plasmon are entirely ruled by scatterings with acoustic phonons, the photoinduced plasmon predominantly transfers its energy to the strongly coupled hot optical phonons, which explains the experimentally-observed tenfold increase of the plasmon linewidth. The present study paves the way for an in-depth theoretical comprehension of plasmon temporal dynamics in novel two-dimensional systems and heterostructures.


Introduction
Understanding, and thus mastering, temporal dynamics of charge carriers in graphene and related quasi-two-dimensional materials is pivotal, but highly challenging task in material science. Many recent studies were devoted to explore the time evolution of the laser-excited electrons by means of time-resolved photoemission [1][2][3][4][5] and pump-probe optical absorption spectroscopies [6][7][8][9][10][11][12][13] in graphene and graphite in order to reach the aforesaid goal. Precise time scales of ultrafast electron interactions were extracted, in particular, electron-electron scattering was shown to rule the dynamics below, while the coupling with the optical phonons (OP) above ∼50 fs [1,[3][4][5][6]. However, underlying microscopic processes still remain largely unexplored, mostly due to a lack of accompanying first-principles methodology that can quantitatively capture these features.
Photoinduced plasmon excitation, i.e. collective electron oscillations under highly non-equilibrium condition, [14,15] is one such ultrafast phenomena that requires further insights. In graphene and graphene-based heterostructures, two-dimensional plasmons show quite exceptional features, e.g. electrical tunability [16,17] and low losses [18][19][20], making these materials promising building blocks for optoelectronic and plasmonic devices. Recently, the relaxation dynamics of laser-induced graphene plasmon was monitored with unprecedented temporal and spatial resolution by using antenna-based near-field nanoscopy [21,22]. High electron temperatures achieved in these experiments are increasing the energy of graphene plasmon while concurrently increasing (decreasing) its linewidth (lifetime) [22]. The former was explained in terms of increase of the Drude weight, or equivalently increase of thermally excited electron-hole pairs, with elevated electron temperature, while the origin of the latter remains unresolved. Since the relaxation of equilibrium graphene plasmon was shown to be governed by the electron-phonon coupling [19,20,23,24], mainly coupling with graphene acoustic phonons (AP) [19,20,23], it was 2. Theoretical methods

Three temperature model with ab initio input parameters
In order to simulate the laser-induced electron dynamics (i.e. electron-phonon thermalization) the three temperature model [29] with ab initio input parameters is utilized (see also references [30,31]). Note that the recent time-resolved photoemission experiments [1,4,5] have demonstrated that the nascent electron distribution is formed into Fermi-Dirac distribution almost instantly, i.e. within 25-50 fs, which justifies the use of the effective electron temperatures for the pump-probe time delays larger than 50 fs. In general, for a more quantitative description of electron dynamics below 50-100 fs one would need to adopt a more rigorous approaches, such as the classical [35][36][37] and quantum [38] kinetic theories or the Keldysh Green's function technique [39,40]. Nevertheless, in graphene and graphite, due to mentioned quasi-instant electron thermalization as well as rapid electron-phonon coupling [41] both electron [1,2,42] and phonon [43][44][45] temperatures are well defined quantities already at the early stage of photoinduced electron dynamics. Within the present model, the temperature of graphene is divided among three subsystems, i.e. electron temperature (T e ), temperature of the strongly coupled OP (T OP ), and the remnant temperature that mostly belongs to AP (T AP ). The energy flow between electron and phonon degrees of freedom is then dictated by the electron-phonon coupling [27], while the thermalization between two phonon subsystems goes via anharmonic coupling. Since the quasi-instant formation of the Fermi-Dirac electron distribution is assumed, electron-hole recombination processes are not considered and only single (electron) temperature represents the dynamics of the electronic subsystem.
The out-of-equilibrium dynamics of the three subsystems can then be simulated by the following coupled equations: The corresponding specific heats C e , C OP , and C AP are defined as: where N(ε) is electron, while F OP (ω) and F AP (ω) are phonon density of states. The electron-phonon relaxation rates G ν are obtained from electron-phonon calculations via [28,30,31]: where λ ν are the electron-phonon coupling strengths for ν = OP and ν = AP, and ω 2 ν are the corresponding second moments of the phonon spectrum, Here α 2 F(Ω) is the Eliashberg function, which quantifies the amplitude of electron-phonon coupling for each phonon energy Ω. The separation between the strongly coupled OP and weakly coupled AP are defined by introducing a cutoff λ c for the mode-resolved electron-phonon coupling strength λ qμ . Namely, the modes that satisfy λ qμ < λ c belong to the weakly-coupled, while the modes with λ qμ > λ c to the strongly-coupled subsystem. For the clean separation λ c = 1 is used. Furthermore, τ = 3.3 ps is the anharmonic scattering time between the OP and AP modes [46]. In addition, I(t) describes a femtosecond Gaussian pump pulse with fluence F and duration t p , while β determines the energy density of the pulse [1,31] and accounts for the fact that only a fraction of the pump fluence is absorbed by the sample. Here β = 400, which gives very good agreement with the experiment. All the input parameters (except β) required in equations (1)-(3) have been determined entirely from first principles, i.e. by using density-functional theory and density-functional perturbation theory [47].
The shortcomings of the present effective temperature model are that the pump-pulse term I(t) does not depend on laser frequency, but only on the laser power, as well as that it cannot describe the early electron dynamics before the formation of the Fermi-Dirac distribution. However, recent experiments have clearly demonstrated that the electron thermalization due to electron-electron interactions is almost instant (i.e. it happens within the first 25-50 fs) [1,4,5], and thus the discrete vertical electron excitations excited with the photon energy are almost instantly distributed into the hot Fermi-Dirac distribution, which can be described with the single electron temperature T e . Therefore, the present model should be able to describe well the electron dynamics above 50 fs and for the small excitation energies.

Time-dependent electron excitation spectrum with phonon-assisted processes included
In order to accurately describe electron excitation processes and optical absorption in doped single-layer graphene, the present paper employs a first-principles current-current formalism within linear response [24,26]. Within this approach the bare current-current correlation function π αα (where α is the polarization direction) is screened with Coulomb interaction by the following Dyson equation π αα = π αα + π αα ⊗ D ⊗ π αα (where D is the photon propagator). The corresponding optical conductivity at the given frequency ω is σ αα (ω) = iπ αα (ω)/ω, while the excitation spectrum that includes collective modes is given with S (q, ω) ∝ Im π αα (q, ω) /ω [24,26]. To include electron-phonon coupling into excitation spectrum as well as plasmon-phonon coupling the current-current correlation function is first separated into interband and intraband contributions and then the electron-phonon coupling is incorporated into the latter channel [24]. D Novko The interband term has the following form where ε nk is the electron energy for band index n and electron momentum k, f(ε nk ; T e ) is the Fermi-Dirac distribution function at electron temperature T e , γ inter is the phenomenological relaxation parameter for interband transitions [26], and A is the area of the unit cell. The current vertex (i.e. electron-photon coupling function) is defined as, where φ nk are the Kohn-Sham ground state wavefunctions.
The intraband current-current correlation function with electron-phonon interaction included can be written as [24,48] I have obtained this result by applying the Holstein theory for normal metals, where the conductivity is calculated by means of a diagrammatic analysis and solving the Bethe-Salpeter equation for the electron-phonon interaction [49][50][51]. Here the electron-phonon coupling is incorporated through dynamical electron-hole pair energy renormalization and decay rate parameters, i.e. λ ep (ω; {T}) and γ ep (ω; {T}), both of which are dependent on electron and phonon temperatures, i.e. T e and T ph [24,25,50]. In particular, dynamical electron-hole pair decay rate due to electron-phonon coupling is [24,25,32] The energy renormalization parameter λ ep (ω; {T}) is obtained by Kramers-Kronig transformation of γ ep (ω; {T}). By dividing the electron and phonon degrees of freedom into three subsystems defined with T e , T OP , and T AP , as in the previous subsection, one can separate γ ep (ω; {T}) in the three parts [32] where The Eliashberg function is accordingly also divided into contributions coming from the AP and OP phonon subsystems, i.e.
Finally, the excitation spectrum with phonon-assisted decay channels and full temperature dependence included is calculated as where π αα is the full (intraband and interband) current-current response tensor that includes Coulomb screening and electron-phonon interactions as presented above.
The time dynamics of the electron excitations and plasmons is then obtained by correlating equation (19) with equations for the time evolution of the effective temperatures equations (1)-(3) [32,52].
Note that the present methodology includes both intraband and interband electronic transitions and thus accounts for both low-and high-energy long-wavelength excitations. Therefore, it could be applied both for lightly-and heavily-doped samples, as well as to any kind of plasmonic materials such as metals, semimetals, or even gapped systems [53,54].

Computational details
The ground-state calculations were done by means of the quantum Espresso (QE) package [55] with a plane-wave cutoff energy of 50 Ry. Norm-conserving pseudopotentials were used with the LDA exchange-correlation functional [56]. A 24 × 24 × 1 Monkhorst-Pack grid was used for sampling the Brillouin zone (with Gaussian smearing of 0.02 Ry). Electron and hole dopings were simulated by adding and removing, respectively, the electrons and introducing the compensating homogeneous charged background. Phonon energies and electron-phonon matrix elements, which are needed for obtaining the input parameters for the three temperature model equations (1)-(3) as well as for the decay rates due to electron-phonon coupling equations (14)- (18), are obtained by using density functional perturbation theory [47] as implemented in QE. The Eliashberg function α 2 F(Ω) is calculated on 400 × 400 × 1 and 40 × 40 × 1 electron and phonon momentum grids, respectively. The electron momentum summations in the intraband and interband current-current correlation functions are done on a 400 × 400 × 1 grid including up to 20 unoccupied electronic bands.

Results and discussion
Ultrafast electron dynamics is explored here for the hole-doped graphene where the Fermi energy is ε F = −250 meV, as it is the case, e.g. for graphene adsorbed on SiC surface [1] (note that the conclusions of the paper would be the same for lightly electron-doped graphene as in the graphene/SiO 2 system [22] due to electron-hole symmetry in the low-energy region of band structure). The corresponding electron and phonon band structures are shown in figures 1(a) and (b). Additionally, figure 1(b) shows the results for the Eliashberg function that measures the degree of the electron-phonon coupling with energy resolution. As is well know [57], electrons are strongly coupled to OP at K and Γ point of the Brillouin zone, i.e. ω 0.16 eV where Eliashberg function shows two prominent peaks, while only weakly coupled to the rest of the modes (mostly low-energy AP). Accordingly, the plasmon excitations below 0.16 eV are weakly coupled to AP [23], while at higher energies plasmon decay is mostly due to the strong coupling with OP [24]. The graphene plasmon dispersion ω pl and the corresponding decay regions due to intraband and interband transitions (Landau damping) [58] are shown in figure 1(c). Figures 2(a) and (b) further depict the energy renormalization 1 + λ ep (ω) and decay rate γ ep (ω), respectively, of the electron-hole pairs due to coupling with phonons as a function of temperature [25,48,50,51]. The Fermi energy is here chosen to be ε F = 300 meV for the sake of comparison with the experiment [20] (the rest of the results are for ε F = −250 meV as mentioned earlier). For plasmon energies, i.e. ω = ω pl , the quantities 1 + λ ep (ω pl ) and γ ep (ω pl ) are equivalent to the plasmon energy renormalization and plasmon decay rate due to coupling with phonons [24]. The results show that the plasmon energy is insignificantly renormalized due to electron-phonon coupling, with very small temperature modifications. The corresponding plasmon broadening is as well small (especially for ω pl < 0.16 eV, where electrons mostly couple to AP) but notable. For the temperature range presented here, the temperature-induced change in plasmon broadening is more pronounced for the lower energies since AP are more easily excited D Novko  [20]. The corresponding decay rate without the influence of the boron nitride is shown with blue circles. with temperature than the high-energy OP (i.e. ω OP k B T). Figure 2(c) shows the temperature dependence of the plasmon decay rate when ω pl = 110 meV. The results are in very good agreement with the decay rate extracted from the recent experiment done on high-mobility graphene [20]. Further analysis demonstrates that the decay rate of the equilibrium graphene plasmon and its temperature dependence predominantly comes from coupling with AP, in agreement with the previous studies [19,20,23].
I turn now to the study of the non-equilibrium condition, i.e. of the ultrafast electron dynamics in graphene by means of the three temperature model with ab initio input parameters. The resultant time evolution of these temperatures is shown in figure 3(a), where the laser with fluence of F = 8 J m −2 and duration of 30 fs excites the system at time delay t d = 0. Right after the laser excitation, the electron temperature T e abruptly increases up to almost 3000 K and subsequently decays due to coupling with OP, in good agreement with the experiment [2]. Consequently, the temperature of the strongly coupled OP elevates above 1000 K, while the AP remain almost at the same temperature. Since the energy exchange rate between phonon and electron baths is proportional to λ ω 2 /C ph (where λ is the electron-phonon coupling strength, ω 2 is the second moment of the phonon spectrum, and C ph is the heat capacity of the relevant phonon subsystem) [27], the obtained dramatic difference between the OP and AP temperatures comes not only because the OP are coupled more strongly to the electrons than the AP (λ OP ω 2 OP > λ AP ω 2 AP ), but also because the strongly coupled OP subsystem consist of only very few modes around Γ and K points of the Brillouin zone (see figure 1(b)) and thus C OP C AP . Such laser-induced hot OP scenario was already discussed in various spectroscopy studies [1,3,6]. Figure 3(b) shows the ensuing charge density modifications δn as a function of time delay. As laser excites the system, electron-hole pair concentration increases both in intraband (electron and hole are in the same band) and interband (electron and hole are in different bands) channels. In fact, for the lightly-doped graphene with ε F = −250 meV it turns out that δn intra δn inter . The laser-induced changes of the electron-phonon coupling are shown in figures 3(c) and (d). In particular, time-dependant modifications of the electron-phonon decay rate δγ ep and energy renormalization parameter 1 + λ ep are depicted for three different excitation energies ω. Since 1 + λ ep is already small for the equilibrium situation, it is not surprising that 1 + λ ep experiences only minor changes as a function of pump-probe time delay and excitation energy. One can also note that laser excitation reduces the energy renormalization parameter as a function of time. On the other hand, for the same laser conditions, the photoinduced decay rate modifications δγ ep are relatively high [i.e. δγ ep γ ep (T = 300 K)] and actually follow the variations of both T e and T OP (see the discussion below). Also, the overall intensity of δγ ep over time is bigger for smaller excitation energies (i.e. for ω = 0.1 eV and 0.2 eV, compared to ω = 0.4 eV). This is because for ω 0.2 eV the equilibrium value of γ ep is small and includes mostly the contributions from the weakly-coupled AP modes, while when both T e and T OP are elevated the probability of scattering on the OP, which are strongly coupled to electrons, increases significantly. A more rapid increase of electron-OP scattering probability when T 300 K can, for example, be seen in figure 2(c). However, when ω 0.2 eV the probability of scattering on the OP is less altered for the present laser conditions, since ω > T e , T OP . Therefore, the relative modification of damping rate is less pronounced for ω 0.2 eV compared to ω 0.2 eV. The time dynamics of charge density and decay rate are relevant for comprehending the transient optical absorption, i.e. photoconductivity, as well as the relaxation mechanisms of non-equilibrium plasmons. Figure 3(e) depicts the time evolution of optical absorption σ 1 (ω) (i.e. the real part of optical conductivity σ) up to 0.7 eV. The observed modifications are due to photoinduced variations of T e , and also due to changes in decay rate δγ ep (note that σ 1 (0) ∝ n/γ ep and σ 1 (ω > 0) ∝ nγ ep /(ω 2 + γ 2 ep ) [24,48,51]). In particular, the low-energy part of σ 1 (ω) decreases significantly around t d = 0 due to changes in δγ ep and then it starts to increase back to its equilibrium value (see the inset in figure 3(e)). Contrary, for higher excitation energies up to around interband threshold 2|ε F | the photoinduced modifications of σ 1 (ω) are first increasing and then decreasing. The time evolution of photoconductivity Δσ 1 below, around, and above interband threshold 2|ε F | is also depicted in figure 3(f). In all these three energy regimes, the photoconductivity Δσ 1 shows different time behaviour, which is simply due to modifications of the interband onset when T e is altered.
Note that the experimentally observed negative and positive values of graphene photoconductivity over different values of ω are widely discussed and analyzed in literature [6][7][8][9][10][11][12][13]. This shows that the present methodology could be also useful for studying the non-equilibrium properties in pump-probe optical spectroscopy and non-equilibrium transport via σ αα (ω; {T}) and σ αα (ω = 0; {T}), respectively. However, here the focus is more on the ultrafast plasmons and the corresponding dynamics in different energy regimes. The time modulations of graphene plasmon properties (e.g. energy loss and Drude weight) under optical pumping were in fact discussed recently [59][60][61][62][63], however, the detailed ab initio study on the corresponding plasmon loss channels is still lacking.
Figures 4(a)-(f) show the variations in the plasmon energy δω pl , plasmon linewidth δγ pl , and the spectral function S(q, ω) as a function of time delay t d for two different values of wavevector q, i.e. for two different energy regimes (note that ω pl ∝ √ q). Note that the time variations of the plasmon energy of 2D system can be calculated as [64]: while the corresponding time evolution of the plasmon linewidth as Also note that damping due to electron-phonon coupling enters the intraband, while the Landau damping the interband part of time-dependent optical conductivity. The obtained time dynamics of the photoinduced graphene plasmon is in good agreement with the experimental observations [21,22]. Namely, the results show that the graphene plasmon is blueshifted and broadened upon the laser excitation. Also, δω pl and δγ pl are notably more pronounced for larger plasmon energies ω pl (larger wavevectors q). In order to understand these transient features of graphene plasmon, it is necessary to dissect different contributions to ω pl and γ pl . First of all, note that in general both intraband and interband excitations determine the total value of plasmon energy. Here, the laser-induced modifications of plasmon energy predominantly come from the increase of electron-hole pair concentrations in the intraband channel n intra (see figures 4(a) and (d)). In other words, laser elevates T e , i.e. increases the number of thermally excited electrons in the intraband channel (see also figure 3(b)), which in turn increases the Drude weight and thus the plasmon energy [21,22]. As already discussed above, the renormalization of the plasmon energy due to electron-phonon coupling is insignificant (see also figure 3(d)).
The processes underlying the plasmon decay rate δγ pl as a function of time delay are a bit more complex than the processes ruling δω pl . For plasmon energies ω pl ≈ 0.1 eV, it turns out that the increase of the plasmon broadening δγ pl is entirely ruled by electron-phonon coupling. What is intriguing and actually at odds with the previous assumptions [22], is that the photoinduced plasmon decays mostly due to scatterings with the OP modes (cf green and purple dashed lines in figure 4(b)). This is also in contrast with the decay mechanisms of the equilibrium plasmon, for which the scattering with the AP modes is the main loss channel (see figure 2(c)). The present analysis also shows that part of the electron-OP scattering contribution to the δγ pl is induced by the elevated T e and part by the elevated T OP (cf purple and black dashed lines in figure 4(b)). Namely, the laser heats the electrons (i.e. elevates T e ), which in turn increases the electron phase space for the electron-OP scattering, especially when T e > ω OP . In addition, hot electrons transfer the excess energy to the strongly-coupled OP, i.e. T OP rises substantially, which increases the number of thermally excited OP and therefore increases δγ pl . All in all, the results show that the plasmon-OP scattering is behind the experimentally-observed plasmon broadening under non-equilibrium conditions. The coupling with the OP modes is generally much more stronger than with the AP modes. However, for equilibrium case the energy conservation condition ω pl ω OP must be met in order to activate this strong coupling. On the other hand, under the strong laser excitations, this energy conservation condition loosens up considerably and graphene plasmons with energy ω pl < ω OP couple strongly with the D Novko hot OP modes. When the plasmon energy is more closer to the interband threshold 2|ε F |, i.e. ω pl ≈ 0.2 eV (see figure 4(f)), the decay is mostly due to the Landau damping (electron-hole pair interband excitations) [65], which increases with T e (see figures 3(e) and (f)).
Finally, I would like to emphasize that the presented results are valid for the freely suspended graphene samples. Quantitatively, but not qualitatively, different results are expected for graphene-substrate systems [66][67][68][69][70][71]. Namely, for graphene-dielectric system one expects only the slight redshift of the plasmon energy [66], while for the graphene-metal contacts the plasmon dispersion obtains linear momentum-dependence and several new features such as higher intensity and lower damping rates [69,70] due to surface-modified electron-electron interactions [71]. Nevertheless, the present results regarding the non-equilibrium plasmon decay rate are expected to hold also in these graphene-substrate systems.

Conclusion
Temporal dynamics of non-equilibrium plasmon under intense laser excitation was explored in lightly-doped graphene by means of density functional and density functional perturbation theories. By considering plasmon-phonon coupling, decay rates of graphene Dirac plasmon were studied under both equilibrium (i.e. electron and nuclear degrees of freedom are thermalized) and ultrafast (i.e. electrons and phonons are thermally excited and have disparate energies) conditions. Due to available phase space and energy constraints, equilibrium plasmon with energy ∼0.1 eV and at ambient temperature is mostly damped due to scatterings with the acoustic phonon modes. Interestingly, the photoinduced non-equilibrium graphene plasmon is, on the other hand, underlain by completely different damping mechanism. Namely, the pump laser pulse increases the population of hot electrons, which in turn transfers the large portion of energy to the strongly coupled OP, creating hot OP bath. Under such out-of-equilibrium condition phase space for electron scatterings and population of optical phonon modes are increased, as well as the energy constraints are loosened. Consequently, the broadening of the non-equilibrium plasmon is increased immediately after the laser excitation, mostly due to coupling with the hot OP. The corresponding strong plasmon energy renormalization is explained in terms of photoinduced Drude weight increase. For energies ∼0.2 eV damping pathways of non-equilibrium plasmon are again different and are due to photoinduced interband excitations (Landau damping). Present study might also help elucidate energy-transfer mechanisms under optical pumping in novel quasi-two-dimensional materials that support collective plasmon modes, like metallic or doped semiconducting transition metal dichalcogenides [53,[72][73][74], topological insulators [33,34], borophene [75], buckled-honeycomb lattices [76,77], or layered electrides [78,79].