Theory of high-order harmonic generation for gapless graphene

We study the high-harmonic spectrum emitted by a single-layer graphene, irradiated by an ultrashort intense infrared laser pulse. We show the emergence of the typical non-perturbative spectral features, harmonic plateau and cut-off, for mid-infrared driving fields, at fluences below the damage threshold. In contrast to previous works, using THz drivings, we demonstrate that the harmonic cut-off frequency saturates with the intensity. Our results are derived from the numerical integration of the time-dependent Schrödinger equation using a nearest neighbor tight-binding description of graphene. We also develop a saddle-point analysis that reveals a mechanism for harmonic emission in graphene different from that reported in atoms, molecules and finite gap solids. In graphene, the first step is initiated by the non-diabatic crossing of the valence band electron trajectories through the Dirac points, instead of tunneling ionization/excitation. We include a complete identification of the trajectories contributing to any particular high harmonic and reproduce the harmonic cut-off scaling with the driving intensity.


Introduction
High-order harmonic generation (HHG) is a remarkable process resulting from the interaction of physical systems with intense electromagnetic radiation. In contrast to most of the conventional photon up-conversion mechanisms, HHG is not based in the multiphoton excitation of atomic bound-state transitions, but on the dynamics of the unbound electrons. The non-perturbative character of HHG has as a distinguishable signature: the emergence of a plateau-like structure in the harmonic intensity spectrum, followed by an abrupt cut-off [1]. This plateau is characterized by a dependence of the harmonic intensity with the harmonic order, q, much weaker that the qth-power predicted by the perturbation theory. This structure extends the harmonic emission up to thousands of harmonic orders, therefore conveying the possibility of generating coherent extreme ultraviolet (XUV) or even soft x-ray radiation [2,3].
High-order harmonics have been observed from a wide variety of targets, including gases, solids, and plasmas. In the latter case HHG arises from the collective response at relativistic intensities [4]. For the atomic or molecular gases and finite-gap solids, high harmonic radiation is produced by intense pulses, yet below the relativistic limit. In these systems, HHG shares some common basic principles [5]. In particular, harmonics are generated by electrons initially bounded, that are promoted into free or quasi-free states by tunneling excitation. Once freed, the electrons are accelerated by the electric field to, subsequently, release the acquired kinetic energy in the form of high-frequency radiation. HHG leads to a wide range of applications, from imaging and spectroscopy with sub-femtosecond resolution to sources of XUV/soft x-ray coherent pulses [6,7]. Recently, HHG from solids has burgeoned a great interest, mainly motivated by the quadratic scaling of the harmonic conversion efficiency with the density of the target, as a result of the coherent nature of the process [5,[8][9][10][11][12][13][14].
Despite those similarities, the first experiment of HHG in solid state [8] noticed substantial differences in the laws governing the spectral plateau and cut-off frequency, compared to the atomic case. In atoms and molecules, the cut-off frequency scales with the product of the laser intensity and the squared wavelength, w l µ I of an ionized electron recolliding with the parent ion. In contrast, the cut-off frequency has been found to scale linearly with the driving field amplitude for semiconductors [8]. This kind of behavior may remind the HHG from simple two-level systems, where the cut-off energy reflects the maximum instantaneous Stark shift [15,16]. However, two-level systems are not sufficient to explain the complex dynamics introduced by the intraband contributions. Several theoretical models have been used to explain this linear response, including the intraband dynamics of the quasi-free states, such as semiclassical models for intraband currents [17], semiconductor optical Bloch equations [18][19][20][21], and Bloch equations using Wannier localized functions [22]. HHG in 2D materials with field polarization perpendicular to the plane has been recently explored theoretically [23].
Remarkably enough, electron tunneling plays a fundamental role as the first step in the HHG process for systems presenting an energy gap, or ionization potential, larger than mid-IR photon energy. In this sense, the gapless structure of graphene presents a new scenario, and the actual details of the mechanism underlying HHG remain unaddressed.
The gapless band structure of single-layer graphene (SLG) allows the optical resonant excitation at all frequencies, up to the vacuum ultraviolet. This conveys graphene particular optical properties, as a strong broadband linear response with a comparatively large optical absorption (>2%) of visible light [24], and a strong nonlinear response for THz radiation [25,26]. While the generation of the second harmonic is forbidden in the dipole approximation, due to the centrosymmetric structure of ideal SLG, it can be observed in stacked samples [27]. Third-order nonlinearities are found to be also remarkably strong in SLG, with nonlinear susceptibilities several orders of magnitude above those of transparent materials, and of the same order as in other resonant materials, such as metal nanoparticles. The third harmonic has been observed in few-layer graphene for transitions occurring near the K and M points of the Brillouin zone (BZ) [28][29][30]. Nonlinear effects arising from induced plasma have been predicted to enhance HHG [31,32]. HHG in gapless graphene has not been reported until recently [33,34], observing the generation of up to the 9th harmonic using a mid-infrared driving laser. This demonstrates the experimental feasibility of producing HHG in monolayer graphene.
Despite these promising characteristics, current theoretical models for gapless graphene are not complete in describing the nonlinear dynamics induced by short intense laser pulses. Some models have been used for dopped graphene with a small energy gap [35], which do not include the important effects induced by the singular coupling near the Dirac cone [36]. Other models accounting for the Dirac singularity do not consider the entire BZ, therefore limiting the cut-off scaling predictions [37]. Questions as the role of the spectrally ubiquitous resonance, the singular transition matrix elements at the Dirac points, the cut-off scaling, or the underlying physical mechanism in HHG remain unveiled.
In this paper we present new theoretical results of HHG in SLG induced by few-cycle laser pulses at infrared wavelengths. On one side, we develop an exact computational method that overcomes the numerical instabilities associated with the singular dipole transition elements near the Dirac points. We use this method to integrate the time-dependent Schrödinger equation (TDSE) for different driving field wavelengths and intensities, in order to characterize the emergence of the non-perturbative spectral features (a plateau in the harmonic intensities, followed by a cut-off). Finally, we develop an approximated model that retains the main contributions to HHG and reveals the fundamental mechanism of harmonic emission: excitation near the Dirac points and creation of an electron-hole pair that subsequently emits a high-frequency photon upon recombination, when the electron and hole overlap in direct space. Therefore, as a main conclusion, we demonstrate that graphene presents its particular mechanism of HHG, in which the first step is initiated by the non-adiabatic crossing of the valence electron trajectories through the Dirac points. This first step is radically different from the tunneling ionization/excitation process found in atoms, molecules and finite gap materials [19]. Our trajectory analysis accounts exactly for the cut-off scaling with the intensity.
The manuscript is organized as follows. In section 2 we present our theoretical framework. In section 3 we present our method for the numerical integration of the TDSE in SLG, showing HHG spectra and the scaling of the cut-off harmonic frequency with the intensity for different wavelengths. In section 4 we develop our approximate description and derive the physical conditions for the harmonic emission, as well as the rules governing the harmonic cut-off frequency. Finally, we present our conclusions in section 5.

Theory
We consider the standard tight-binding description of SLG [38], where the wavefunctions of the conduction (+) and valence (−) bands can be written as linear combinations of Bloch states from two neighboring sublattices, show the structure of SLG in the direct and reciprocal space, respectively. The Dirac points are located in points K and K′ in the reciprocal space. K is related to K′ by inversion symmetry. The Γ point, at the center of the BZ, corresponds to the maximum gap ;17.8 eV in the tight-binding description. The driving electric field is assumed polarized in the y direction, and aimed perpendicularly to the graphene layer. Dirac cones result from the overlap of atomic p z orbitals, and are coupled when the field polarization is parallel to the graphene sheet. The interaction with a field perpendicular to the layer has been studied recently elsewhere for the case of 2D materials, and found to resemble HHG in atoms [23]. The vectors in equation (1) are eigenstates of the system's Hamiltonian H G with energies . The interaction of the laser pulse F(t) with the system is described by the time-dependent Hamiltonian H(t)=H G +V i (t), where V i (t)=−q e F(t) y is the electric field coupling within the dipole approximation.
During the interaction, the time-dependent wavefunction can be expressed as a superposition of the eigenstates (1): , and projecting, we find the following set of coupled two-level equations [39]. For ultrashort pulses (<30 fs full-width at half maximum, FWHM), it is possible to remain in the Schrödinger representation instead of using the density-matrix formalism, as the carrier collisional diffusion has characteristic times of several tens of femtoseconds [40][41][42]. Furthermore, it is also possible to draw a parallelism between equation (4) and the standard strong-field formulation for atoms [43], replacing the band structure, E + (k), by the parabolic shape of the free electron energy, including the energy term −F(t)D y as a phase-shift, and taking C − (k, t)=1.
Equations (4) and (5) include interband and intraband couplings, the intraband contribution corresponding to derivatives in the reciprocal space. The computational complexity introduced by these gradients can be removed introducing the kinetic momentum: Figure 1. Scheme of graphene's structure in (a) direct and (b) reciprocal space. The unit cell is composed by two carbon atoms (red and blue). In the reciprocal space, points K and K′ correspond to the location of the Dirac cones. In our study the laser pulse is linearly polarized along the y direction and it is aimed perpendicularly to the graphene layer (x-y plane).
t e ( ) t A being the vector potential. The equations (4) and (5) are then reduced to: This set of equations corresponds to a driven two-level parametric oscillator.

Numerical results
The numerical integration of equations (7) and (8) is of difficult convergence due to the singular values of D y (k) near the Dirac points. It is possible to recast these equations and renormalize the terms containing D y (k) using the following transformation, which effectively undoes the diagonalization of H G : According to Larmor's semiclassical formula, the harmonic emission is given by the dipole acceleration. The mean value of the dipole can be written as the sum of two contributions, corresponding to intraband ( « d ) and Using (9) and (10), the total dipole can also be written as We have integrated numerically equations (11), (12) and (16), considering few-cycle driving pulses at different intensities and wavelengths. The driving pulses are modeled using an 8 cycle (full width) sin 2 temporal envelope, 2.9 cycles FWHM. The total field is defined as 8 , and 0 elsewhere, where F 0 is the field amplitude, T is the period, and ω 0 =2π/T the field frequency. Figure 2 shows the calculated HHG spectra at two different intensities for a driving field of 3 μm wavelength. For the lowest intensity, figure 2(a), the HHG spectrum shows the typical perturbative behavior: a monotonous decrease of efficiency with increasing harmonic order. This behavior changes drastically for driving fields at higher intensity, figure 2(b), with a clear emergence of a spectral plateau, extending up to a maximum (cut-off) frequency. It has been reported that the inclusion of additional energy bands results in the appearance of secondary plateaus, with higher cut-off frequencies, but with efficiencies smaller in various orders of magnitude [20], not affecting therefore our conclusions. It is interesting to note also that the spectral content in figure 2 is quite rich, even for the lower intensity case. This results from the interference of two different contributions to the same radiated energy, namely, the intraband and interband components of the emission dipole, and from the different electron-hole pair's trajectories leading to the same harmonic, as we shall describe below. For the sake of comparison, we show in the inset of figure 2(b) the HHG spectrum corresponding to the hydrogen atom, driven by a 800 nm laser with the same pulse width. We will show that the mechanism of electron excitation in graphene differs substantially from the tunnel ionization in atoms. The efficiency in the latter case is much lower, and HHG typically occurs for laser intensities one or two orders of magnitude above those required in graphene. Note that, although the spectra is quite complex in both cases, atoms present a more regular structure after the cut-off frequency. This difference will be discussed in the following section.
We show in figure 3 the scaling of the cut-off frequency with the driving field intensity for laser wavelengths 800 nm, 1.6 μm and 3.0 μm. Note that for the three cases the cut-off frequency saturates at the largest intensities, and the saturated cut-off corresponds to a photon energy of ;17.8 eV, i.e. the maximum gap in the BZ. This saturation has also been described in gap semiconductors [5]. Unlike the other physical systems mentioned in the introduction (atoms, molecules or two-level systems), there is no simple law (linear or quadratic) describing the scaling of the cut-off frequency with the field amplitude.
The use of arbitrary large intensities is precluded by the damage of the sample. We indicate with a colored area in figure 3 the intensities for which the material is expected to have damage, assuming a damage fluencethreshold of 150 mJ cm −2 [44]. Observe that for longer wavelengths it is possible to reach the saturated cut-off at intensities below the damage threshold. Note however that, as the laser period increases, the decoherence due to carrier collisions, typically of several tens of fs, becomes a limitation. This tradeoff indicates that the HHG production is more favorable at mid-IR wavelengths in the short-cycle regime.
Our numerical findings reflect some fundamental properties of the harmonic emission: (i) emergence of a spectral plateau at large intensities, (ii) a nontrivial cut-off dependence with the intensity, and (iii) a limiting photon energy close to the maximum energy gap of the material. In the following section, we develop a saddlepoint approximation model (SPAM) that provides a simple description of HHG in SLG. The SPAM allows to identify the fundamental mechanism and predicts the spectral cut-off energies, in excellent agreement with our numerical calculations.

Underlying mechanism
Our description starts with the integral form of the dynamical equations (7) and (8): where we define the action term In SLG, k ( ) D y t presents singularities near the Dirac points. We model this behavior defining k ( ) D y t as an impulse function , where k D is the coordinate of the Dirac point at K. The following discussion applies as well to the K′ Dirac point, due to the inversion symmetry, assuming a change of sign in the electric field amplitude. Using the impulse ansatz, we focus on the dynamics of the valence electron initially located at the point k, that follows the quiver trajectory in the BZ given by k t . Using the impulse function, the probability amplitudes (17) and (18) can be then written as . High-harmonics in solids, including graphene, are dominated by interband dynamics [20,34]. Using equations (19) and (20), the complex amplitude of the interband dipole (15) is given by , . The amplitude of the qth harmonic of the interband dipole is then given by the Fourier transform of the k-space integral  Following [19,43] (in the context of solids), we can identify that the main contributions to the integral in (22) are the stationary phase points. Therefore, the harmonic emission takes place predominantly at those times that fulfill the equation, This condition can be interpreted as the emission of a photon by interband relaxation at the point k t , where the gap is resonant with the emitted photon. This condition, together with equations (19) and (20), leads to a simple picture of the qth order harmonic emission in SLG: an electron, initially at point k in the valence band, quivers in momentum space with amplitude k t due to the interaction with the field; at time t D,k the electron's trajectory approximates the Dirac point, and the electron is promoted to the conduction band through the non-adiabatic interaction with the singular matrix element D y , leaving a hole; finally, the electron and the hole oscillate until time t, when they recombine emitting a photon resonant with the band gap at k t . This mechanism for SLG differs from the one described for finite-gap semiconductor solids [19], as the excitation in the latter is dominated by tunneling, rather than by the non-adiabatic crossing. The SPAM predicts that only those states with initial momentum vertically aligned near the Dirac points can be promoted to the conduction band, as these are the only trajectories that cross the Dirac points. As these states are also vertically aligned with the Γ point, the maximum energy of the photon emission will be limited by the gap at this point (17.8 eV), which matches with the value of the saturated cut-off energy mentioned above. We obtain an additional condition by applying the saddle-point analysis in momentum space to (22), which leads to ò ò where the terms represent the velocity of the electrons. The time integral of the velocities corresponds to classical electron trajectories. The interpretation in terms of semiclassical trajectories is at the core of the present understanding of HHG from atoms [43], and it is also extended to solids in [5] or used for the description of the dynamics of hot free carriers in graphene [30]. In this semiclassical framework the interpretation of equation (24) is straightforward: once the electron-hole pair is created at time t D,k near the Dirac point, both are driven by the field; condition (24) states that the pair recombines emitting a photon at time t, only when their trajectories intersect in real space. Figure 4(b) shows the classical trajectories of two electrons, with different initial positions in the BZ (k A and k B ), which cross the Dirac point at times t D k , A and t D k , B . The trajectories correspond to a 3 μm wavelength pulse, with intensity 10 11 W cm −2 . Note that for the trajectory A, the electron and the hole created at t D k , A meet in direct space at the final time t. In contrast, in B, the electron and the hole do not meet at t. According to equation (24), the photon emission at t is effective only for the case A. Figure 4(a) shows a map of the energy gap at the moment t where the photon may be emitted, as a function of the time in which the electron-hole pair is created, t D,k . The points A and B in the map correspond to the cases shown in figure 4(b). We have colored in red the points (t D,k , t) corresponding to electron-hole pairs, created at t D,k , that are driven to overlap in the same direct-space unit cell at the emission time t. These colored areas are, therefore, in compliance with both conditions (23) and (24), and correspond to the situations when the harmonic photon is effectively emitted. According to this map, point A corresponds to an electron-hole pair created at , which emits a w 10 0 photon when recombined at  t T 0.9 . Note that, as it happens in B, other points in the map may potentially emit higher frequency harmonics (up to w 14 0 , i.e. the maximum gap energy attained during the electron-hole excursion). However, the photon emission is not possible since the electron-hole pairs are spatially apart at t. The maximal photon energy is, therefore, given by the energy maxima in the map of figure 4(a), constrained to the colored zone, which corresponds to the trajectories ending with an electron-hole intersection and therefore, having the possibility of recombination. In figure 4, the cut-off harmonic corresponds to point A, with an energy w 10 0 , smaller than the maximum gap ( w 14 0 at the Γ point). In figure 3(a) we compare this cut-off prediction (red circle) with the results of the exact integration of equations (4) and (5) (blue diamonds). As it can be noticed, this cut-off prediction is in excellent agreement with the results of the numerical calculations.
We show in figure 5 the SPAM energy maps at 3 μm wavelength for the four additional intensities used in figure 3(a). We observe that for increasing intensity, the topology of the colored region, i.e. the condition for the electron-hole overlapping at the emission time, becomes more complex. It is interesting to note that the same harmonic may be emitted by two or more trajectories within the same driving field's half cycle. These interfering contributions are also found in HHG by atoms and molecules [45], although in SLG their characteristics, i.e. the number of trajectories, initial and final times, etc reveal a more intricate scenario. Note, in particular, that figures 5(c) and (d) reveal three or more path contributions for the maximum (cut-off) frequency. In contrast, there is a single electron trajectory emitting cut-off harmonics in atoms and molecules. As a result, in these latter systems, the HHG spectrum becomes regular at the cut-off, while in graphene remains, still, complex (see figure 2(b) and its inset). As described, we can use these maps to identify the cut-off harmonic order at each intensity. The predicted values are plotted in figure 3(a), and show an excellent agreement with the results of the numerical computations. Figures 3(b) and (c) show the cut-off extracted from the SPAM maps in comparison with the TDSE results, for different intensities and 1.6 μm and 800 nm wavelengths, demonstrating a similar accuracy of the SPAM for the whole range of parameters.
Before concluding, let us comment the harmonic emission associated with the intraband transitions from this simplified perspective. Substituting equations (17) and (18) in (14), and discarding the terms that vanish after the integration over the BZ (due to the graphene's inversion symmetry), we find, This corresponds to the time-dependent dipole associated to the position of a valence electron in direct-space, initially with pseudomomentum k, that undergoes an excitation to the conduction band at time t D,k . The intraband dipole, therefore, is also a source of harmonics due to the non-parabolic shape of the electronic bands. Previous studies in finite gap crystals show that the intraband contribution dominates the harmonic emission of low order harmonics, while the interband dominates for the higher frequency spectrum [19,20].

Conclusions
We have studied the HHG in SLG irradiated by intense few-cycle infrared laser pulses. The study is two-fold. First, we integrate exactly the TDSE using a non-diagonal basis to circumvent the numerical instabilities associated to the Dirac points. Our method allows, therefore, to account the singular non-adiabatic coupling near the Dirac cones without approximation. Our numerical results demonstrate the emergence of the nonperturbative signatures in the harmonic spectrum (intensity plateau, followed by an abrupt cut-off) at midinfrared wavelengths, below the damage threshold. In contrast to the atomic case, we show that in graphene the harmonic cut-off frequency saturates with increasing intensity. In a second part, we develop a SPAM to unveil the basic mechanism of HHG in graphene. According to the SPAM, harmonics are generated by electron-hole pairs produced during the non-adiabatic crossing of the electron trajectories near the Dirac points. This mechanism differs from the tunnel excitation giving rise to HHG in atoms, molecules and solids with finite gaps. Once generated, the electron and hole are driven by the field until their trajectories intersect in real space, allowing recombination, and thus emitting a high-frequency photon. This electron-hole recombination mechanism is analogous to the rescattering process found in atoms and molecules. In this latter case, however, the hole remains static in the ion, while in solids, it quivers with the electric field. Our SPAM reproduces the scaling of the harmonic cut-off frequency with the driving field intensity with excellent agreement with the numerical results, and fully characterizes the trajectories, their number, and the electron-hole creation and recombination times leading to the emission of harmonics. The theory presented in this manuscript can be extended to other two-dimensional materials. Our SPAM results are encouraging for the development and design of new materials, and of methods to control the electron dynamics to tailor the HHG emission [46]. This work paves the way for further investigations in the understanding of the nonlinear optical response of gapless materials and in the manipulation of electron carriers at the petahertz domain. The formalism presented here is also specially suitable to study the interplay of nonlinearly polarized pulses and the role of strong electron-hole correlations in the production of HHG [47].