Terahertz radiative coupling and damping in multilayer graphene

The nonlinear interaction between intense terahertz (THz) pulses and epitaxial multilayer graphene is studied by field-resolved THz pump–probe spectroscopy. THz excitation results in a transient induced absorption with decay times of a few picoseconds, much faster than carrier recombination in single graphene layers. The decay times increase with decreasing temperature and increasing amplitude of the excitation. This behaviour originates from the predominant coupling of electrons to the electromagnetic field via the very strong interband dipole moment while scattering processes with phonons and impurities play a minor role. The nonlinear response at field amplitudes above 1 kV cm−1 is in the carrier-wave Rabi flopping regime with a pronounced coupling of the graphene layers via the radiation field. Theoretical calculations account for the experimental results.


Introduction
The electronic band structure of graphene is characterized by double cones in the vicinity of the K and K points (Dirac points) of the Brillouin zone. The energy increases for the upper or conduction band or decreases for the lower or valence band linearly with the distance from the nearest Dirac point. Accordingly, the electron velocity is independent of the energy and equal to the Fermi velocity [1] v F = 10 6 m s −1 . Another consequence of graphene's band structure is the 1 Author to whom any correspondence should be addressed. If an electron on such a trajectory reaches a region in k space where the THz field is in resonance with the interband transition frequency, combined inter-and intraband trajectories will occur (lower arrow). constant optical absorption over a large frequency range, having a value of π α = 0.023 per layer (α: fine structure constant) [2,3]. The constant absorption is a result of the cancellation of two effects, a density of states proportional to and a transition dipole moment inversely proportional to the transition frequency. In the terahertz (THz) range this results in a huge interband transition dipole moment, at 2 THz d = e × 80 nm (e: elementary charge).
Because of the huge interband transition dipole moment, one enters the regime of extreme nonlinear optics [4] or the nonperturbative regime in the interaction of electrons with the THz field already for moderate electric field amplitudes. One indication for this behaviour is the interband Rabi frequency = d E/h, which is equal to the THz carrier frequency at a THz field amplitude of 1 kV cm −1 . As a consequence strong radiative coupling should occur in multilayer graphene samples where the layer separation is much smaller than the THz wavelength. So far, radiative coupling in multilayer graphene has remained mainly unexplored.
Ideal graphene has a Fermi energy E F = 0, i.e. at zero temperature the lower band is completely filled and the upper band is completely empty. In this case THz absorption is connected with interband transitions, bringing electrons from the valence into the conduction band ( figure 1(a)). For a nonzero Fermi energy, interband transitions are blocked up to photon energieshω = 2|E F |. On the other hand, electrons in the conduction band (figure 1(b)) (or holes in the valence band) give rise to intraband transitions. Therefore, it is not a priori clear whether a finite Fermi energy will result in an increase or a decrease of the absorption in the THz range. The same problem arises when free carriers are generated by optical excitation. There have been several reports on the change of THz absorption after excitation in the near-infrared (NIR) [5][6][7], in the mid-infrared (MIR) [8], and in the THz range [9][10][11]. While measurements performed on single-layer graphene [6,7,10,11] with Fermi energies between −170 and −300 meV show an absorption decrease, measurements on epitaxial multilayer graphene [5,8,9] (Fermi energies around 25 meV) show an increase of absorption after excitation.
In this paper, we address the nonlinear THz response of multilayer graphene in pump-probe experiments with a field resolved detection of the pulses. We demonstrate an enhanced transient THz absorption after excitation by ultrashort pulses in the THz, the MIR and the NIR spectral range. The decay of enhanced absorption occurs on a time scale of a few picoseconds. We show that this decay is governed by the predominant interaction with the external THz field, i.e. by radiative coupling and damping. Theoretical calculations account quantitatively for this behaviour.

Experiment
The experiments are based on collinear two-dimensional (2D) THz spectroscopy, a technique introduced recently [12,13]. In the following, we present pump-probe results in the time domain measured in the following way. The output from a femtosecond Ti:sapphire oscillator-amplifier system is split and sent into two GaSe crystals to generate two phase-locked THz pulses A and B by optical rectification [14]. The two pulses with a spectral maximum at 2 THz and a mutual delay τ interact with the graphene sample. The transmitted electric-field transients are measured by electro-optic sampling [15] in a ZnTe crystal. Two choppers synchronized to the pulse repetition rate allow for measuring the electric-field transients E AB (t, τ ) (both pulses A and B are incident on the sample), E A (t) (only pulse A) and E B (t, τ ) (only pulse B). The nonlinear signal emitted from the sample, which is proportional to the nonlinear current in the sample [16], is the difference and depends on the delay τ and the real time t. The graphene samples grown by C-face epitaxy on SiC (from Graphene Works) contain 45 layers. Electronic coupling between the layers is negligible [17]. The spectrally integrated linear transmission of the sample at T = 300 K was 0.36, as measured with THz pulses of a low electric field amplitude of 0.2 kV cm −1 . A detailed analysis shows that the linear THz transmission spectra are dominated by intraband (Drude) absorption with an average Fermi energy of E F = 24 meV and a Drude scattering time of T s = 40 fs.

Experimental results
Time-resolved THz transients measured with excitation pulses in the THz, MIR and NIR spectral range are summarized in figure 2. The electric field of the THz probe pulse and the nonlinear signal E NL are plotted as a function of real time t, the pump-probe delay was τ = 0. Independent of the pump wavelength the nonlinear signal displays a phase opposite to the driving electric field, indicating induced absorption.
In the following we concentrate on the dependence of the induced absorption on the delay τ between the two THz pulses. To extract from our 2D data the customary pump-probe signals, we perform a 2D Fourier transform of the nonlinear signal E NL (t, τ ) and filter out the pump-probe signal in the frequency domain [9,13]. After a Fourier transform back into the time domain we get the pump-probe signal E PP (t, τ ). With this the differential transmission is given by The probe field E pr can be either E A or E B , depending on which signal was chosen in the frequency domain. Results for the differential transmission T /T 0 (τ ) obtained from (2) are shown for different sample temperatures in figure 3(a) and for different amplitudes of the THz field in figure 3(b). In all cases the transmission decreases (induced absorption) upon pumping. The transmission recovers with time constants τ d between 1.8 and 4.5 ps, as obtained from single-exponential fits (dashed lines in figure 3). These times are surprisingly short compared to the experimental data in [5,8]. The decay times increase with decreasing temperature and with increasing THz field amplitude. The magnitude of the differential transmission increases with decreasing temperature.

Discussion
We first present a numerical estimate underlining the importance of radiative coupling in multilayer graphene. Considering only interband transitions in a fully coherent picture and assuming that all polarizations are parallel, the amplitude of the coherent interband current at a point in k space is e v F . Integrating this current over the part of the Brillouin zone with interband transition frequencies between ω 1 and ω 2 yields the total one-dimensional (1D) current density per layer of J = e(ω 2 2 − ω 2 1 )/(16π v F ). This current density corresponds to an emitted field of E em = N Z 0 J/(1 + n sub ) with N = 45 the number of layers, Z 0 = 377 the vacuum impedance and n sub = 3.1 the refractive index of the substrate [18]. The emitted field leads to an energy loss rate of L = (1 + n sub )E 2 em /Z 0 . On the other hand, the energy content, i.e.hω/2 integrated over the same part of the Brillouin zone is Thus, for a fully coherent emission the energy content has a radiative lifetime T rad = W/L. Using the frequency range relevant for our experiment (ω 1 = 2π × 1 THz, ω 2 = 2π × 3 THz), we find T rad = 1 ps. This value is in the same order of magnitude as the observed decay times of induced absorption in figure 3, suggesting a strong influence of radiative coupling on the measured transients.
We now address a key result of our experiments, the observation of induced absorption after THz, MIR and NIR excitation (cf figures 2 and 3). This finding has to be reconciled with the observation of an absorption decrease in samples with |E F | > 100 meV. NIR pump pulses at 1.5 eV excite electrons from the valence into the conduction band for all samples studied so far, since Pauli blocking of optically coupled states is absent even for the highest |E F | ≈ 300 meV [5][6][7]. For samples with E F ≈ 25 meV [8], this is also true for MIR excitation. The generated electrons and holes thermalize on a time scale of the order of 100 fs into quasi-Fermi distributions with an elevated carrier temperature [19], a process generating a free carrier population close to the K points. This distribution leads to a reduced interband absorption because of Pauli blocking and to an enhanced intraband absorption because of enhanced electron and hole densities. The higher carrier temperature, on the other hand, has the opposite effect, namely an increased absorption on interband transitions for which Pauli blocking is reduced, and a reduced absorption from intraband transitions as-at higher temperature-carriers move into regions of k space where they have larger effective masses 2 . It is not a priori clear which of these effects dominates. The difference in the sign of the transient absorption change for different E F can be explained by considering that intraband transitions are the dominant cause of THz absorption. For samples with large |E F |, i.e. a large density of free carriers, the main effect of pumping is the temperature increase, leading to reduced absorption, whereas in samples with low |E F | the main effect is the increase in carrier density, leading to increased absorption. Our sample with an average |E F | = 24 meV exhibits the latter behaviour.
Our data display an amplitude of induced absorption that increases with increasing amplitude of the pump field (figure 3(b)) and decreasing lattice temperature ( figure 3(a)). With increasing pump field, the free carrier densities in the valence and conduction band and, thus, the intraband absorption strength are enhanced. At low lattice temperature, the initial density of free carriers is small and pumping induces a strong increase in carrier density and differential absorption. At higher lattice temperature, the substantially higher initial free carrier density and the partial (Pauli) blocking of interband THz excitation of extra carriers result in a smaller differential signal.
Before addressing the picosecond decay times of the nonlinear THz signal, we briefly consider a potential influence of scattering processes on the observed kinetics. The energy range considered here is far below the optical phonon frequency of 45 THz [20] and, consequently, optical phonon emission plays a minor role. For sample temperatures up to 300 K, the thermal optical phonon population is very small, ruling out optical phonon absorption as well. Acoustic phonons have much lower energies, but their coupling to electrons is quite weak, resulting in decay times of a THz driven current on a time scale of 100 ps [21], much longer than the decays we observe. The recently predicted supercollisions [20,22], i.e. disorder-assisted electron-acoustic-phonon scattering, lead to a strong temperature dependence and to decay times of the same order of magnitude as our results. However, this mechanism leads to a decrease of the decay time τ d with increasing pump field, in contrast to our observations. If electron-electron scattering would be the dominant scattering process, τ d should decrease with increasing pump field and carrier density, again in contrast to our experimental results.
To account for the picosecond decays of the THz signals, radiative coupling within the multilayer sample was considered within a model which incorporates the electronic bandstructure and the coupling of carriers to the THz field but neglects scattering processes. To describe radiative coupling, we include not only the action of the electric field on the electrons, but also the influence of the electrons on the electric field. It turns out that this simple model is sufficient to describe all our results. While the cone approximation [23], in which the valence and conduction bands are mirror images of each other, is sufficient to account for the observed induced absorption and for the absence of photon-echo signals [9], it fails in predicting the picosecond decay times. A quantitative treatment of the latter requires an electronic band structure that correctly includes the threefold symmetry around the K points and the differences between upper and lower bands.
To describe the threefold symmetry of the band structure around the K point correctly, one has to consider the hexagonal lattice symmetry of graphene. One possibility to do this is a tight-binding calculation [1,24], in which the Bloch functions are composed of atomic orbitals centred at the positions of the atoms. In the simplest way just one type of atomic orbital is used, which leads to two bands as the primitive unit cell contains two C atoms. Another possibility is the pseudopotential method [25][26][27]. The pseudopotential form factors can be obtained either from first-principles calculations or from a fit to experimentally determined energies at certain points in the Brillouin zone (empirical pseudopotential). For our calculations we used such empirical pseudopotential band structure (see the appendix for more details). The pseudopotential method yields THz transmission values in agreement with the experiment (at low temperatures the THz transmission is determined by the quantum conductivity [2,3,28] σ 0 ). A comparison of the different band structure calculations is shown in figure 4. They agree at the K point and show a difference in carrier energy of less than 0.2 eV up to a tenth of the Brillouin zone around the K point. At larger k vectors, the energy differences are significant.
Since the stack of graphene layers is much thinner than all optical wavelengths involved, the local electric field is the same in all layers. The local electric field consists of the incident electric field E in (t), of the electric field emitted by the electrons E em (t) (in this way the radiative coupling is included), and of the electric field reflected at the interface between the substrate and the graphene E sub (t). In this geometry the local electric field is equal to the transmitted field [16] E tr (t) = E in (t) + E em (t) + E sub (t).
(3) We use the long-wavelength approximation [29], i.e. neglect the influence of the magnetic field and assume that the electric field in the sample is spatially homogeneous, i.e. only timedependent. This gives a spatially constant vector potential For a sample consisting of N graphene layers the emitted field is given by [30] E em (t) = − N Z 0 2 J (t).
The 2D current density J (t) is the current density of a single graphene layer. If the graphene layers are on a substrate (in our case hexagonal SiC [17]), the refractive index jump at the interface also emits a field determined by the refractive index n sub of the substrate [18]. Both the field emitted by the stack of graphene layers and by the substrate modify the driving vector potential of the system. As a result the influence of E em (t) on A(t) causes coherent radiative damping of the entire system. For the linear THz response of graphene, i.e. for low THz fields, the present model based on pseudopotential band-structure calculations gives identical results to those of the widespread conductivity model [2,5,28,[31][32][33][34][35]. To account for the finite broadening of intra-and interband transitions and to achieve irreversible absorption, in the linear regime we introduced a phenomenological friction force. In the nonlinear regime, there is no need for such a phenomenological damping, since losses due to radiative coupling and damping become dominant.
Intraband transitions in the absence of scattering are governed by the acceleration theorem, i.e. d k/dt = −e E(t)/h (figure 1). A THz pulse polarized in x-direction 3 moves an electron back and forth in k space along a trajectory given by k x (t) = k x (0) − e A x (t)/h determined by the transient vector potential A x (t) (k y remains constant). For constant k y the band structure E(k x , k y = const. = 0) as a function of k x − K x follows two hyperbolas (valence and conduction band). For an electron residing initially at k with |k y − K y | > ω THz /2v F resonant interband transitions are not possible so that the intraband current (motion of the carrier along the upper double arrow in figure 1(c)) dominates the response to the applied THz field.
This behaviour changes for electrons residing initially at k vectors with |k y − K y | < ω THz /2v F . In this case the change in k x introduced by the electric field brings the electrons to points in k space where the difference between valence and conduction band is resonant to the THz frequency. As a result, they perform trajectories as visualized by the lower double arrow in figure 1(c). Because of the extremely strong transition dipole the light-matter interaction is inherently in the nonperturbative regime and dominates over the interactions underlying scattering processes of the carriers.
The trajectories shown in figure 1 are only possible if the initial state is filled and the final state is empty. Thus, if the valence band is completely filled, the intraband current is zero. The same is true for an empty conduction band. Considering the case of a finite Fermi level E F at zero temperature, for low field amplitudes interband transitions are not possible if 2|E F | >hω THz . However, because the intraband motion happens simultaneously, even with a finite E F interband transitions become possible provided the field strength is large enough (e v F max(|A(t)|) > |E F |). For the field strengths used in our experiment this means that Fermi levels |E F | ≈ 25 meV, as obtained from linear absorption measurements, yield essentially the same results as E F = 0.
Our theory provides k-resolved populations and polarizations (figure 5). These results illustrate in detail why the radiative damping is faster when less energy is deposited in the sample, as is the case at higher lattice temperatures. In figure 5(a), the calculated electron distribution ρ 22 (k x , k y ) in the upper Dirac cone around K is shown after interaction with a THz pulse polarized along the x-direction with a very small electric field amplitude of 30 V cm −1 . For symmetry reasons the distribution functions in the cone around K are mirror images of those around K . The corresponding carrier distribution ρ 22 (k x , k y ) is essentially the product of the spectral density of the THz pulse with the modulus squared of the dipole matrix elements. After interaction with the driving pulse, the electron distribution in the upper cone ρ 22 (k x , k y ) has a very short lifetime, since radiative damping in the 40-layer graphene sample is extraordinarily strong. Thus, for very low THz fields the model predicts a decay within the pulse duration.
The situation changes in the nonperturbative regime. In figure 5(b) we plot ρ 22 (k x , k y ) after interaction with a THz pulse with an amplitude of 3 kV cm −1 . The generated electron distribution ρ 22 (k x , k y ) shows a speckle-like structure due to the chaotic phase distribution in k space of the concomitantly generated coherent interband polarizations shown in figure 5(d). The reason for the chaotic phases are the different Rabi frequencies at each point in k space. While for low field amplitude (figures 5(a) and (c)) all electrons in the upper band have the same phase and thus can emit coherently, for high field amplitude (figures 5(b) and (d)) this is only true for the electrons within a speckle. Our theory predicts a radiative damping rate that decreases with the number of speckles within this chaotic ensemble. The number of speckles increases with the stored energy. Since the signal strength is proportional to the stored energy, one expects the behaviour experimentally observed in figures 3(b) and 6. With increasing temperature, the picosecond decay of the enhanced absorption becomes faster and the amplitude of the pumpinduced absorption decreases ( figure 3(a)). The experimentally observed temperature-dependent damping rate originates from the temperature-dependent initial electron populations. The sample converts the initial coherent THz radiation at 2 THz into broad-band chaotic radiation (still coherent) with frequencies up to a cutoff frequency given by the highest transition frequency encountered during the intraband trajectory. The new frequency components can in turn be reabsorbed and reemitted via interband transitions that were initially not in resonance with the THz pulse. This multi-frequency character of the overall THz field emitted by the sample explains the absence of photon-echo signals in the nonlinear response [9]. The interband Rabi frequency is different at each point in k space, due to both the k-dependent interband dipole moment and the different electric field amplitudes E(ω) at different interband transition frequencies. Thus, a coherent rephasing of the different contributions by interaction with the THz driving field is impossible and the total photon-echo signal is negligible.
In the nonperturbative regime, the coupling of the carriers to the external THz field is much stronger than the interaction among the carriers and interactions of carriers with the lattice. As a consequence, carrier dynamics are governed by the highly nonlinear interaction with the THz field rather than by scattering processes. Under such conditions, the single-particle picture of carrier dynamics applied here represents a realistic description of the experimental scenario.
Finally, we compare the picosecond decays of enhanced absorption measured at different lattice temperatures ( figure 3(a)) and the experimental results of [8] with the predictions of our theory. For a direct comparison, the transient change of the intraband oscillator strength σ intra (τ, ω) dω/ σ 0 dω was calculated with (4) of [5]. In figure 6, this quantity is plotted as a function of pump-probe delay for different THz pump pulses and lattice temperatures. In panel (a) we compare the calculated change of the intraband oscillator strength for a 1 ps pump pulse athω = 8 meV and a field amplitude of E pu = 5 kV cm −1 with our experiments performed at temperatures of 50 and 300 K. 4 In panel (b) we calculated such a transient for a longer (10 ps) pump pulse at a higher photon energy,hω = 20 meV, with a field amplitude of E pu = 2.4 kV cm −1 , corresponding to the = 0.4 µJ cm −2 transient from figure 3(b) of [8]. These results show that the decay depends strongly on the energy deposited in the sample. Longer pulses (figure 6(b)) deposit more energy in the sample, leading to longer decays than shorter pulses (figure 6(a)). The good agreement of theory and experiments reveals the prominent role of radiative damping for energy relaxation of the the carriers while scattering processes play a minor role under the present conditions.

Conclusions
Using collinear 2D THz spectroscopy we have studied the nonlinear dynamics of electrons in multilayer graphene. Instead of the expected bleaching due to the saturation of interband absorption, the pump-probe signals show induced absorption, which decays with time constants of a few picoseconds. The time constants decrease with increasing temperature and increase with increasing energy density deposited. Our results are fully explained by theoretical calculations using a pseudopotential band structure for graphene and including radiative coupling and damping. Under the present conditions of strong radiative coupling in a multilayer sample, fast decoherence processes or other scattering mechanisms play a minor role for the nonequilibrium response to the external THz field.
The sum extends over all reciprocal lattice vectors G. The solutions of (A.1) have according to the Bloch theorem [36] the form ψ n, k ( r ) = exp(i k r ) u n, k ( r ) (A. 3) with u n, k ( r ) having the periodicity of the lattice. Thus, similar to (A.2), ψ n, k ( r ) can be written as The eigenvalues of H ( k) for all k in the first Brillouin zone then yield the electronic band structure.
For every numerical calculation one limits the infinite sum in (A.4) to a finite subset of reciprocal lattice vectors (total number Z ), so that H ( k) becomes a Z × Z matrix. In this case at most Z Fourier components of the potential can contribute.
Since graphene is atomically thin, a 2D treatment is sufficient. With the in-plane lattice constant a the primitive reciprocal lattice vectors can be taken as (A.6) The simplest pseudopotential calculation that yields the threefold symmetry around the K [ K = (2 a * + b * )/3] point needs three reciprocal lattice vectors, namely g 1 = 0, g 2 = a * and g 3 = a * + b * , which are on an equilateral triangle with K in the centre. For this calculation the three possible Fourier components of the potential have the same value because of lattice symmetry. This value is determined so that the band structure calculation yields the same energy for the lower band at the point as the experiment [37]. In this band structure calculation the Fermi velocity v F is nearly independent of the form factor. It is determined essentially by the free-electron mass m 0 and graphene's lattice constant a, v F = h/(3 m 0 a) = 10 6 m s −1 . The results shown in figure 4 were obtained from a calculation using 21 reciprocal lattice vectors, also oriented symmetrically around the K point. One advantage of this calculation is that it is very easy to incorporate the matter-light coupling by replacingh k in (A.5) withh k + e A(t) (see section 4).