Plasmon-assisted high-harmonic generation in graphene

High-harmonic generation in condensed-matter systems is both a source of fundamental insight into quantum electron motion and a promising candidate to realize compact ultraviolet and ultrafast light sources. While graphene is anticipated to efficiently generate high-order harmonics due to its anharmonic charge-carrier dispersion, experiments performed on extended samples using THz illumination have revealed only a weak effect. The situation is further complicated by the enormous electromagnetic field intensities required by this highly nonperturbative nonlinear optical phenomenon. Here we argue that the large light intensity required for high-harmonic generation to occur can be reached by exploiting localized plasmons in doped graphene nanostructures. We demonstrate through rigorous time-domain simulations that the synergistic combination of strong plasmonic near-field enhancement and a pronounced intrinsic nonlinearity result in efficient broadband high-harmonic generation within a single material. Our results support the strong potential of nanostructured graphene as a robust, electrically tunable platform for high-harmonic generation.


Introduction
High-harmonic generation (HHG) is an extreme nonlinear optical phenomenon first observed by driving atomic gases with intense ultrashort light pulses [1,2].The harmonic intensity remains surprisingly large up to a high order of the pulse carrier frequency, stimulating applications for HHG as a source of ultraviolet and x-ray radiation [3][4][5], as well as in the generation of attosecond pulses [6][7][8], which has enabled tomographic imaging of molecular orbitals [9] and the exploration of subfemtosecond dynamics in chemical reactions [10].
Recent observations of HHG from condensed-matter systems [11][12][13][14][15] are currently attracting much interest not only in the pursuit of new solid-state optical technologies, but also in the underlying physics of HHG in bulk crystals and its analogy with atomic gases.Indeed, while HHG from individual atoms is well-understood as the coherent emission produced by the optically induced tunneling ionization of an electron, its acceleration by the driving field, and the subsequent recollision with its parent ion [16,17], the picture becomes less clear in crystalline media, where collective effects associated with the high density of electrons and their interaction with the lattice significantly complicate the generation process.As expected, HHG in solids is found to depend strongly on the electronic band structure and the interplay between inter-and intraband transitions [11,12,14,18,19].
The linear, gapless dispersion relation of graphene elec-trons [20,21] garners strong interest in the nonlinear optical response of the atomically thin material, which recent experiments demonstrate to be intrinsically large [22][23][24][25][26][27].On the theory side, monolayer graphene is expected to produce intense HHG in the THz regime [28,29], attributed to complementary inter-and intraband charge carrier motion at low temperatures and doping levels.Unfortunately, recent experiments report either no evidence [30] or only a weak effect [31] associated with the generation of low-order harmonics from multilayer graphene for currently available THz illumination intensities.This situation could be improved by using more intense sources at higher frequencies, and further relying on enhanced graphene-light interaction mediated by localized plasmon resonances.
Graphene plasmons [32][33][34][35][36], which provide an efficient way to couple the carbon layer with impinging light, are capable of generating intense local electric fields that are essential to trigger nonlinear optical phenomena.This near-field enhancement, in combination with the highly anharmonic response of graphene [29,37,38], is predicted to give rise to large optical nonlinearies [39][40][41][42][43]. Importantly, these plasmons only exist in highly doped graphene, while their frequency is strongly dependent on the doping level [32][33][34][35][36]. Electrical gating thus provides a mechanism to tune the harmonic generation in graphene to the desired frequency range.
Here we predict that highly efficient HHG takes place in doped graphene nanostructures when the incident light is tuned to their localized plasmons.Specifically, we obtain harmonic intensities that are orders of magnitude higher than in other materials.Additionally, no sharp cutoff is observed with harmonic order.Our results are based on nonperturbative time-domain numerical simulations of the nonlinear optical response of graphene using arXiv:1609.09794v1[cond-mat.mes-hall]30 Sep 2016 two complementary approaches: a random-phase approximation (RPA) description of the single-particle density matrix within a tight-binding (TB) model for the electrons of ribbons and finite islands [39]; and the solution of the single-particle Bloch equations for massless Diracfermions (MDFs) in extended graphene, complemented by a classical electromagnetic (CEM) description of the self-consistent field produced by the illuminated nanostructure (see Methods).We find both approaches to be in excellent agreement at intensities below the saturable absorption threshold.Our prediction of highly efficient HHG assisted by coupling to graphene plasmons suggests applications to a wide range of nonlinear photonic technologies, including tunable sources of broadband attosecond light.

Results
In practice, cumbersome laser amplification schemes are usually needed to reach the extreme electromagnetic field intensities required to generate high-order harmonics.To overcome this limitation, plasmonic nanostructures have attracted considerable interest as in situ electric field enhancers for HHG in gaseous media [45][46][47][48].As illustrated schematically in Fig. 1a, we propose that compact, efficient HHG can be realized in graphene by combining the intense near-field enhancement associated with graphene plasmons with the intrinsically high nonlinear optical response of this material.The appeal of graphene as a nonlinear optical material stems in part from its linear charge carrier dispersion with electron wave vector k at low energies, ε k = hv F |k|, where v F ≈ c/300 is the Fermi velocity.In the single-particle MDF description of doped monolayer graphene, neglecting interband electronic transitions, this linear dispersion relation leads to a maximum achievable surface current density J max = −env F sign{sin(ωt)} when illuminated by a monochromatic field E(t) = E 0 cos(ωt) in the E 0 → ∞ limit [28,37].The current is thus limited by the doping charge-carrier density n.This square-wave profile of the induced current density under intense illumination translates into efficient generation of odd-ordered harmonics (see Fig. 1b).Conversely, in conventional 2D media, for which charge carriers obey a parabolic dispersion relation ε k = h2 k 2 /2m * , the system responds harmonically at the frequency of the driving field, regardless of electron-electron interactions [44].While this comparison favorably portrays graphene as a highly nonlinear optical material, it is important to note that interband optical transitions compensating the large intraband anharmonicity become significant at high intensities, even when the system is driven at frequencies below the Fermi level [28].
Quantitative analysis of plasmon-enhanced HHG in a doped graphene nanoribbon is presented in Fig. 2. The linear optical absorption of the nanoribbon under consideration (20 nm width, E F = 0.4 eV Fermi energy) shows a prominent dipolar plasmon (Fig. 2a), as predicted by TB-RPA atomistic simulations and classical electrodynamics, in excellent mutual agreement.We thus consider HHG produced by incident pulses with central frequency tuned to that plasmon.We present HHG simulations obtained with the MDF-CEM and TB-RPA approaches (see Methods) in Fig. 2b, which shows the spectral decomposition (time-Fourier transform) of the radiative emission intensities for 100 fs incident light pulses with three different peak intensities.Each spectrum is normalized to the maximum value around the fundamental frequency.The corresponding temporal evolution of the graphene induced current is shown in Fig. 2c.Remarkably, high harmonics up to 13 th order are clearly discernible in the emission spectrum even at a relatively low incident peak intensity I 0 = 10 12 W/m 2 .The agreement between MDF-CEM and TB-RPA descriptions is then excellent both in the spectra (Fig. 2b, upper plots) and in the time-resolved induced current (Fig. 2c).The temporal evolution of the induced current tends to follow the profile of the incident Gaussian pulse, although a small time delay of the peak current is observed in the atomistic simulation due to the self-consistent Coulomb interaction, which persists beyond the duration of the pulse on a timescale determined by the inelastic relaxation time τ = 13.2 fs.By raising the peak intensity, the conversion efficiency of high-order harmonics drastically increases in the MDF-CEM picture, while a more modest, yet impressive, enhancement is predicted in the atomistic TB-RPA simulations.Finite-size effects that are included in the atomistic simulations but not in the MDF-CEM description (see Methods) contribute to this discrepancy.Additionally, the plasmonic local-field enhancement is self-consistently described in the TB-RPA approach, but not in the MDF-CEM method.For the high level of doping under consideration, intraband electronic transitions dominate the optical response, particularly at low intensities, while interband transitions reduce the level of anharmonicity, as observed in the temporal profiles of the induced current when comparing MDF-CEM simulations with (center plots) and without (left plots) inclusion of interband processes (Fig. 2c).
The dramatic increase in HHG from localized plasmons in graphene nanoribbons is clearly shown in Fig. 3 by mapping the emission intensity over a wide range of input pulse carrier frequencies, where at each input frequency the response is normalized to its respective maximum at the fundamental harmonic.Noticeable enhancement in harmonic generation appears when the excitation frequencies coincide with the plasmon resonance, which can be tuned actively via electrostatic gating and passively by selecting different ribbon widths.Although yet high-order harmonics appear in the spectra, we restrict our investigation to low photon energies where the tight-binding model for graphene remains valid (i.e., below the π plasmon near 5 eV).In Fig. 3a,b we present results for the doped 20 nm ribbon considered previously, based on atomistic TB-RPA and MDF-CEM simulations, FIG.1: High-harmonic generation (HHG) assisted by graphene plasmons (a) Schematic illustration of a doped graphene nanoribbon illuminated by an intense optical pulse that is resonant with the ribbon transverse dipole plasmon.The latter produces strong in-plane electric-field intensity enhancement (see color scale) that boosts the generation of high harmonics.(b) The low-energy band structures of graphene (upper left) and a conventional 2D semiconducting crystal (upper right) respond differently to a monochromatic light electric field E(t) = E0 cos(ωt): in graphene, the induced current J(t) (lower left) acquires a square-wave temporal profile, which contains all odd-order harmonics in its Fourier decomposition, while the semiconductor responds harmonically at the driving frequency ω.A 2D free electron gas also shows a harmonic response [44].respectively, for 100 fs pulses with 10 12 W/m 2 peak intensity as those considered in the upper panel of Fig. 2b.While atomistic simulations quickly become computationally unaffordable for ribbons wider than a few tens of nanometers, the MDF-CEM approach enables the exploration of HHG in much larger structures, such as the 100 nm-wide ribbon explored in Fig. 3c, which is found to generate plasmon-enhanced high-order harmonics with superior efficiency than the 20 nm ribbons.The red-shifted plasmon resonances found in larger graphene nanostructures naturally lead to higher optical nonlinearities due to their increased proximity to the Dirac point [28].
Although graphene possesses a centrosymmetric crystal lattice, the geometry of a finite nanostructure can be chosen in a manner that breaks inversion symmetry, enabling even-ordered nonlinear response in certain directions.In Fig. 4 we present atomistic TB-RPA simulations of HHG in an armchair-edged 15 nm equilateral graphene nanotriangle for incident light polarized normal to one of the triangle sides.When the nanotriangle is doped to a Fermi energy E F = 0.4 eV and illuminated with pulses resonant with the dominant, low-energy plasmon mode (Fig. 4b), high harmonics of both even and odd orders are generated with a similar efficiency to the previously considered graphene nanoribbon (cf.Figs.3a and 4b).Despite the inversion symmetry of the atomic lattice, a nonzero even-order nonlinear current is produced by a combination of the strong local-field-intensity gradient and the relatively high Fermi wavelength λ F ∼ 10 nm [49], which is commensurate with the size of the triangle.In contrast, only odd-ordered harmonics appear if the nanoisland is undoped (Fig. 4a), as both of these effects (field gradient and long λ F ) are then absent.
Ultimately, we are interested in producing intense high harmonics with moderate incident intensities.With this goal in mind, we analyze the performance of graphene for HHG in Fig. 5 and also compare the results with available experiments in solid state systems.As a first observation, even without the involvement of plasmons, the strong intrinsic nonlinearity of graphene is capitalized in a large relative intensity of high harmonics normalized to the response at the fundamental frequency (Fig. 5a): the relative harmonic emission reaches the values measured in GaSe samples, but using 3-4 orders of magnitude lower pulse fluence.It should be noted that a level of theory similar to the MDF model produces excellent agreement with experiment in GaSe (cf.open and solid triangles in Fig. 5a), thus supporting the predictability of our results, which is also emphasized by the agreement between MDF-CEM and atomistic simulations shown in Figs. 2  and 3.By patterning the graphene into ribbons and tuning the incident light to the dominant dipole plasmon energy, HHG is boosted even more, a result that is particularly evident when analyzing the absolute harmonic intensity of resonant ribbons and extended graphene (Fig. 5b).Incidentally, in contrast to the enhancement observed in doped ribbons by exciting the plasmons, doping is detrimental in extended graphene because the Fermi level is then situated in a region where the difference between parabolic and linear electronic band dispersions is reduced, and so is the nonlinear response.

Conclusion
In summary, we predict that the combination of high intrinsic nonlinearity and strong plasmonic field confinement provided by doped graphene nanostructures under resonant illumination leads to unprecedentedly high HHG conversion efficiencies.Despite the fact that this material is only one atom thick, we show that it outperforms other solid state systems, such as GaSe, for which HHG measurements have been reported.It should be noted that our results are based on a conservative value of the phenomenological electronic relaxation time τ .The availability of high-quality graphene samples, in which τ is an order of magnitude longer, should boost HHG in this material even further.We have focused on relatively low fundamental frequencies, so that the high harmonic energies under consideration still lie within a range for which the optical response is dominated by the π band of graphene.At low intensities, the response is well described by the low-energy, linear-dipersion region of the FIG.4: Nonlinear optical response of triangular nanographenes.We show the emission intensity from doped graphene nanotriangles (armchair edges, equilateral 15 nm side length) as a function of incident and emitted photon energies, calculated within the TB-RPA approach upon pulse irradiation (100 fs FWHM duration, 10 12 W/m 2 peak intensity).We consider both (a) undoped and (b) doped (EF = 0.4 eV) triangles.The upper plots show the linear absorption spectrum and the orientation of the normally incident light polarization (insets).electronic band, which explains the agreement that we find between continuum DFM-CEM and atomistic TB-RPA descriptions.Although future work is required to extend these results to higher photon energies, which will require the involvement of deeper electron bands, we conclude the HHG conversion efficiencies associated with localized plasmons in graphene nanostructures appear to be remarkably high for an atomic layer, indicat-ing a strong potential for developing electrically tunable, ultra-compact nonlinear photonic technologies.

A. TB-RPA simulations
We follow a previously-reported atomistic approach [39,43,50] to simulate the nonlinear optical response of graphene nanostructures via direct time-domain integration of the single-electron density matrix equation of motion, where H TB is a tight-binding Hamiltonian describing the one-electron states of the π band of graphene (one outof-plane p orbital per carbon site with nearest-neighbor hopping energy of 2.8 eV), φ is the self-consistent electric potential including both the impinging light and the Hartree interaction, and a phenomenological relaxation is assumed to bring the system to the relaxed state ρ 0 at a rate τ −1 with hτ −1 = 50 meV (i.e., the relaxation time is τ ≈ 13.2 fs).The density matrix ρ = jj ρ jj |j j | is expressed in the basis set of one-electron eigenstates of H TB , where ρ jj are the sought-after time-dependent expansion coefficients.In particular, we have ρ 0 jj = δ jj f j for the relaxed state, where f j are Fermi-Dirac occupa-tion numbers at the initial temperature T = 300 K.For ribbons, the states are treated as Bloch waves, arranged in bands as a function of their momentum along the direction of translational invariance, and the calculation is simplified by the orthogonality of different bands [43].The induced charge density at each carbon atom position R l is then constructed as ρ ind l = −2e jj ρ jj a jl a * j l , where the factor of 2 originates in spin degeneracy, while the coefficients a jl represent the change of basis set between state j and site l representations.Finally, the timedependent induced dipole and surface current are given by d(t) = l R l ρ ind l (t) and J(t) = ḋ(t), respectively.For ribbons, we normalize these quantities per unit of ribbon length [43].

B. MDF-CEM simulations
In a complementary approach, we model electron dynamics in graphene within the MDF picture by adopting a non-perturbative semi-analytical model [51], in which light-matter interaction is introduced through the electron quasi-momentum π = p+(e/c)A, where p is the unperturbed electron momentum, A(t) = −c t −∞ E(t )dt , and E is the classically-calculated in-plane electric field (see Sec. C).Electron dynamics is governed by the Dirac equation for massless fermions, which can be recast in the form of Bloch equations as [28,29,51] where n p (R, t) and Γ p (R, t) represent the population inversion and the interband coherence, respectively [51].Here, the damping energy hτ −1 = 50 meV is the same as in the TB-RPA approach.These equations describe both inter-and intraband transitions.We solve Eqs.(1) nonperturbatively under the slowly-varying-envelope approximation [51] by expanding Γ p (R, t) = N j=0 Γ + p,j (R, t)e i(2j+1)ωt + Γ − p,j (R, t)e −i(2j+1)ωt and n p (R, t) = n 0 (R) + N j=1 Re n + p,j (R, t)e 2ijωt + n − p,j (R, t)e −2ijωt in harmonic series up to N = 15.The current is then parallel to the local electric field E(R, t) x, while its amplitude is calculated as an integral over momentum-resolved contributions, Finally, the far-field power spectrum of the emitted light is proportional to |ω J(R, ω) | 2 , where J(R, ω) is the time-Fourier transform of J(R, t), and . . .denotes the space average over the graphene structure under examination.

C. Classical electromagnetic simulations
The classical response of graphene nanostructures is simulated by numerically solving Maxwell's equations us-ing the boundary-element method [52] for ribbons and a finite-element method (COMSOL) for triangles.We describe graphene as a thin film (thickness s = 0.5 nm) and permittivity (ω) = 1 + 4πiσ(ω)/ωs, where σ(ω) is the local-RPA conductivity [36,53,54].We thus obtain the linear optical extinction and the near-field distribution.Given the small lateral size of the ribbons and triangles compared with the light wavelength, we adopt a quasistatic eigenmode expansion [55] and only retain one term corresponding to the dominant plasmon in each case.The incident light pulses are taken to have a large duration compared with the optical cycle, so we approximate them by a single carrier frequency (the pulse peak frequency) times a Gaussian envelope.We also use this approximation for the input near-field E of the MDF-CEM approach, with the carrier component classically calculated as explained above.

FIG. 3 :
FIG. 3: Dependence of HHG on incident photon energy in graphene nanoribbons.We show the emission intensity from doped graphene nanoribbons under transverse normal illumination as a function of the incident and emitted photon energies as calculated within the (a) TB-RPA and (b,c) MDF-CEM approaches.The incident pulse has a FWHM duration of 100 fs and a peak intensity of 10 12 W/m 2 , while the Fermi energy is 0.4 eV in all cases.The ribbon width is 20 nm in (a,b) and 100 nm in (c), giving rise to the plasmon energies indicated by the vertical dashed lines.

FIG. 5 :
FIG. 5:Comparison of graphene plasmon-assisted HHG with extended graphene and measured bulk semiconductors.(a) We show the emission intensities of high harmonics normalized to the intensity of the fundamental peak as a function of pulse fluence.The harmonic index is indicated by the color-coordinated symbols and numbers.We present MDF-CEM simulations for a 100 nm-wide graphene nanoribbon doped to EF = 0.4 eV (filled diamonds) along with undoped (EF = 0, open hexagons) and doped (EF = 0.4 eV, filled hexagons) extended graphene, which we compare with theoretical (open symbols) and experimental (solid symbols) results for bulk GaSe taken from Refs.[12] (triangles) and[15] (circles).The pulses (100 fs FWHM duration, 0.158 eV peak energy) are resonant with the ribbon plasmon.(b) Comparison of HHG transmission intensity normalized to the incident intensity for a doped-graphene ribbon array (100 nm width, 200 nm period, EF = 0.4 eV, solid diamonds) and undoped extended graphene (open hexagons).The intensities of transmitted (∝ |Et| 2 ) to incident (∝ |E0| 2 ) light are taken at their peak frequencies.