Expanded view of electron-hole recollisions in solid-state high-order harmonic generation: Full-Brillouin-zone tunneling and imperfect recollisions
AAn Expanded View of Electron-Hole Recollisions in Solid-State High-HarmonicGeneration: Significance of the Full Brillouin Zone
Lun Yue ∗ and Mette B. Gaarde † Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803-4001, USA (Dated: February 25, 2021)We theoretically investigate electron-hole recollisions in high-harmonic generation (HHG) in band-gap solids irradiated by linearly and elliptically polarized drivers. We find that in many casesthe emitted harmonics do not originate in electron-hole pairs created at the minimum band gap,where the tunneling probability is maximized, but rather in pairs created across an extended regionof the Brillouin zone (BZ). In these situations, the analogy to gas-phase HHG in terms of theshort- and long-trajectory categorizations is inadequate. Our analysis methodology comprises threecomplementary levels of theory: the numerical solutions to the semiconductor Bloch equations, anextended semiclassical recollision model, and a quantum wave packet approach. We apply thismethodology to two general material types with representative band structures: a bulk system anda hexagonal monolayer system. In the bulk, the interband harmonics generated using elliptically-polarized drivers are found to originate not from tunneling at the minimum band gap Γ, but fromregions away from it. In the monolayer system driven by linearly-polarized pulses, tunneling regionsnear different symmetry points in the BZ lead to distinct harmonic energies and emission profiles.We show that the imperfect recollisions, where an electron-hole pair recollide while being spatiallyseparated, are important in both bulk and monolayer materials. The excellent agreement betweenour three levels of theory highlights and characterizes the complexity behind the HHG emissiondynamics in solids, and expands on the notion of interband HHG as always originating in trajectoriestunnelled at the minimum band gap. Our work furthers the fundamental understanding of HHG inperiodic systems and will benefit the future design of experiments.
I. INTRODUCTION
The recent experimental observations of high-harmonicgeneration (HHG) in solids [1–6] have contributed to therapid progress of attosecond physics in the condensedmatter phase [7–9]. Solid-state HHG carries excitingprospects for the engineering of compact attosecond lightsources [10–15] and ultrafast spectroscopy methods capa-ble of probing band structures [16, 17], impurities [18–21],and topological features [22–29]. The understanding ofthe HHG process has been aided by accurate computa-tional quantum theories such as time-dependent densityfunctional theory (TD-DFT) [30–32] and semiconductorBloch equations (SBEs) [33, 34]. These theories havehelped to establish that the high harmonics with energiesless than the band gap have large contributions from the intraband currents originating in the electron-hole mo-tion in their respective bands, while above-band-gap har-monics are generally dominated by the interband currentsoriginating in the coupling between the bands. While theSBEs and TD-DFT methods are able to accurately simu-late the HHG process, they can be regarded as numericalexperiments that contain all the relevant physics, and theunderlying physical pictures can be difficult to extract.For this reason the celebrated gas-phase three-step model[35, 36] has been generalized to solids [37–39], and hasbeen shown to provide an intuitive real-space picture forthe interband harmonics: an electron-hole pair is created ∗ lun [email protected] † [email protected] when the external field causes an electron to tunnel fromthe valence band to the conduction band at the mini-mum band gap; the electron and hole are driven apartby the laser; they can recollide when they spatially reen-counter each other, leading to the emission of harmonicradiation. The recollision picture has been instrumen-tal in our fundamental understanding of solid-state HHG[2, 17, 37, 38, 40, 41], as well as other related nonlinearphenomena such as high-order-sideband generation [42–46].The tunneling, propagation and recollision dynam-ics responsible for HHG in solids differ significantlyfrom their counterpart in the gas phase. In gases, thecontinuum-electron is free and its dispersion is quadraticsuch that its group velocity is always along the direc-tion of the canonical momentum. In crystalline solids,however, the quadratic dispersion only holds near certainhigh-symmetry points in reciprocal space. Consequently,the group velocities of the electron and hole are generallymuch more complex, and can even lead to imperfect rec-ollisions where an electron-hole pair recollide even thoughtheir centers are spatially separated [41, 47, 48]. The ex-ponential dependence of the tunneling rate on the bandgap [49] dictates that tunneling occurs with the highestprobability at the minimum band gap. However, dueto the complicated dispersions, electron-hole trajectoriesthat originate away from the minimum band gap couldhave higher chances of recollision and end up dominat-ing the emission process. Similarly, electron-hole pairscreated near different symmetry points in the Brillouinzone (BZ) could lead to drastically different harmonic en-ergies and emission time profiles. The full understanding a r X i v : . [ phy s i c s . a t o m - ph ] F e b of these complexities for HHG in solids are critical for theprobing of the full BZ, as well as the design of new ultra-fast light sources. The original semiclassical recollisionmodel in solids, however, assumes tunneling at the mini-mum band gap and with perfect recollisions, and cannotprovide a full framework for the HHG process in solidsapart from simple cases. Due to these limitations, forexample, the authors in [50] concluded that the recolli-sion picture would fail for solid-state HHG with circularlypolarized fields.In this manuscript, we conclusively show that in manycommon experimental scenarios, and for several typesof materials, the electron-hole pairs created away fromthe minimum band gap not only contribute to, but canstrongly dominate the recollision-driven harmonic emis-sion. In these cases, the understanding of the emissiondynamics in terms of short and long trajectories, well-known from gas-phase HHG, breaks down. This break-down can be due to either the band structure or the laserpolarization, and we provide two examples of current ex-perimental and theoretical interest: HHG in a genericbulk crystal induced by elliptically polarized fields, andHHG induced in a generic monolayer material by linearlypolarized fields. Our analysis comprises three comple-mentary levels of theory: the full numerical solution ofthe SBEs, an extended recollision model, and a modelbased on construction of recolliding electron-hole wavepackets. We find that the novel harmonic emission pro-files are due to the collective emission associated withtrajectories originating in extended regions near differ-ent symmetry points. We show that tunneling from dif-ferent regions can lead to different time-frequency emis-sion characteristics, with impact for HHG-based opticalprobing of the whole BZ. The familiar short and longtrajectories can be recovered in special cases: the bulkirradiated by linearly polarized field, and the monolayerirradiated by fields polarized along specific symmetry di-rections. We provide general rules for when one can ex-pect the underlying physics in solid-state HHG to sub-stantially deviate from that of the gas phase.This paper is organized as follows. Section II containsthe theoretical framework pertinent to this work: TheSBEs are given in Sec. II A, the semiclassical picture isdetailed in Sec. II B, and the electron-hole wave packetconstruction is described in Sec. II C. Section III treatsthe ellipticity dependence of HHG in the representativebulk system zinc oxide (ZnO): the model is defined inSec. III A, the full quantum result from the SBEs arepresented in Sec. III B, the semiclassical analysis withtunneling at Γ is discussed in Sec. III C, and the full-BZrecollision picture is given in Sec. III D. Section IV inves-tigates the HHG in the representative monolayer systemhexagonal boron nitride (hBN): The model is describedin Sec. IV A, the wavelength dependence is explored inSec. IV C, the orientation dependence in Sec. IV B, andthe quantum wave packet analysis in Sec. IV D. Section Vconcludes the paper and provides outlook. Appendix Aprovides details on the derivations of the saddle points equations, and Appendix B includes relevant supplemen-tal figures. II. THEORETICAL METHODS
In this section, we describe the theoretical frameworkpertinent to this work. The numerical solutions to theSBEs can be considered a numerical experiment and isour “full quantum” benchmark result, which the semi-classical recollision and wave packet methods will be com-pared to. In the following subsections, we assume thatrelevant quantities such as the band structures, transi-tion dipole moments, Berry connections and Berry cur-vatures are known in advance, either by employing modelsystems or using commercial solid-state structure codes[51, 52].
A. Semiconductor Bloch equations
The SBEs governing a solid driven by a strong laserreads [33, 34, 37, 53, 54]˙ ρ K mn ( t ) = − i (cid:20) E K + A ( t ) m − E K + A ( t ) n − i (1 − δ mn ) T (cid:21) ρ K mn ( t ) − i F ( t ) · (cid:88) l (cid:104) d K + A ( t ) ml ρ K ln ( t ) − d K + A ( t ) ln ρ K ml ( t ) (cid:105) , (1)with K the crystal momenta in a reciprocal referenceframe moving with A ( t ) ≡ − (cid:82) t F ( t (cid:48) ) dt (cid:48) , F ( t ) the elec-tric field, E k n the band energies, d k mn = i (cid:10) u km (cid:12)(cid:12) ∇ k (cid:12)(cid:12) u k n (cid:11) the dipole matrix elements, (cid:12)(cid:12) u k m (cid:11) the cell-periodic partof the Bloch function (cid:12)(cid:12) φ k m (cid:11) , ρ k mn the density matrix ele-ments, T the dephasing time, and A k n ≡ d k nn the Berryconnections.The total current can be split into the interband andintraband contributions by j ter ( t ) = − (cid:88) K (cid:88) m (cid:54) = n ρ K nm ( t ) p K + A ( t ) mn (2a) j tra ( t ) = − (cid:88) K (cid:88) n ρ K nn ( t ) p K + A ( t ) nn , (2b)with p k mn = (cid:10) φ k m (cid:12)(cid:12) ˆ p (cid:12)(cid:12) φ k n (cid:11) the momentum matrix ele-ments, and the summation over K is over the whole BZ(throughout the text). The HHG spectrum is taken asthe modulus squares of the Fourier transforms of the cur-rents, after weighting by a window function.Throughout this work, we make use of the two-bandapproximation with an initially filled valence band la-belled “v” and an empty conduction band labelled “c”.For notational convenience, we henceforth use d k ≡ d k cv for the transition dipole and ω k g ≡ E k c − E k v for the bandgap. B. Extended semiclassical picture
In this subsection, we go over the semiclassical modelsused in this work. We start by obtaining the saddle pointequations from the SBEs, and then describe our extendedrecollision model. The case for linearly-polarized fieldswas partly discussed in [48]. More info on the details ofthe derivation can be found in Appendix A.
1. Saddle-point equations
For the laser pulses and systems considered in thiswork, the conduction band population during the laseris small [ ρ k vv ( t ) − ρ k cc ( t ) ≈ j ter ( ω ) = (cid:82) ∞−∞ dte iωt j ter ( t ), where the Carte-sian ( µ = { x, y, z } ) current components in the fixed frameare (see the derivation in Appendix A 1) j ter ,µ ( t ) = (cid:88) k R k µ (cid:90) t T κ ( t,s ) e − iS µ ( k ,t,s ) ds + c.c. (3)with T κ ( t,s ) = (cid:12)(cid:12) F ( s ) · d κ ( t,s ) (cid:12)(cid:12) the transition matrix el-ement, R k µ = ω k g | d k µ | the recombination dipole, and κ ( t, t (cid:48) ) = k − A ( t ) + A ( t (cid:48) ). The times s and t can beinterpreted as the excitation and emission times, respec-tively. The accumulated phase in Eq. (3) is (dephasingignored) S µ ( k , t, s ) = (cid:90) ts (cid:104) ω κ ( t,t (cid:48) ) g + F ( t (cid:48) ) · ∆ A κ ( t,t (cid:48) ) (cid:105) dt (cid:48) + α k ,µ − β κ ( t,s ) (4)with ∆ A k ≡ A k c − A k v , α k ,µ ≡ arg( d k µ ) the transition-dipole phases (TDPs), and β κ ( t,s ) ≡ arg[ F ( s ) · d κ ( t,s ) ].The saddle point conditions for S µ ( k , t, s ) − ωt read ω κ ( t,s ) g + F ( s ) · Q κ ( t,s ) = 0 , (5a)∆ R µ ≡ ∆ r − D k ,µ + Q κ ( t,s ) = 0 , (5b) ω k g + F ( t ) · (cid:104) Q κ ( t,s ) + ∆ r (cid:105) = ω, (5c)where the electron-hole separation distance and groupvelocities are∆ r ≡ (cid:90) ts (cid:104) v κ ( t,t (cid:48) ) c − v κ ( t,t (cid:48) ) v (cid:105) dt (cid:48) (6a) v κ ( t,t (cid:48) ) n ≡∇ k E κ ( t,t (cid:48) ) n + F ( t (cid:48) ) × Ω κ ( t,t (cid:48) ) n , (6b)with the Berry curvature Ω k n ≡ ∇ k × A k n , and D k ,µ ≡ ∆ A k − ∇ k α k ,µ (7a) Q k ≡ ∆ A k − ∇ k β k . (7b) Note that in our notation, µ used as a subscript pointsto a scalar quantity, while µ used as a superscript cor-respond to a function variable: for example, D k ,µ is avector that depends on µ .Equations (5a)-(5c) can be interpreted by the followingthree steps in the interband HHG process: at time s , anelectron-hole pair is created by tunnel excitation at thecrystal momentum k ≡ κ ( t, s ); the laser accelerates theelectron and hole with the instantaneous group velocities v κ ( t,t (cid:48) ) n ; the electron and hole recollide at time t with finalcrystal momentum k and relative distance ∆ r , with thesimultaneous emission of high-harmonics with energy ω .The saddle-point equations first proposed by Vampaand co-workers [2, 37, 38] include only the first term onthe left-hand sides of Eqs. (5a)-(5c). The above equa-tions includes additionally (i) laser-dressing of the bandswith F · Q κ ( t,s ) in Eqs. (5a) and (5c); (ii) the anomalousvelocity term [55] in Eq. (6a) involving the Berry curva-tures; (iii) a shift of the recollision condition in Eq. (5b)and the possibility of an additional nonzero electron-hole-pair polarization energy at recollision ( eh -PER) F · ∆ r in Eq. (5c). Physically, the eh -PER constitutes the po-tential energy of the electric dipole comprised of thepositively-charged hole and negatively-charged electronat the time of recollision. We note that in systems withinversion and time-reversal symmetries, the Berry cur-vatures are zero. We mention that similar equations toEq. (5) appear in [56], but without inclusion of the eh -PER.Note that for the all physics to be relevant and consis-tent, under an arbitrary “structure”-gauge transforma-tion (cid:12)(cid:12) u k n (cid:11) → (cid:12)(cid:12) u k n (cid:11) e iϕ k n ( φ k n ∈ R ), Eq. (5) should remainunchanged. While the individual terms in the right-handsides of Eq. (7) generally depend on the gauge-choice, thecomposed quantities, D k µ and Q k µ are shown to be gaugeinvariant in Appendix A 2. The gauge invariance of thesaddle-point equations in Eq. (5) then trivially follows.In many studies of solid-state HHG, linearly polarizeddrivers F ( t ) = F ( t ) ˆe (cid:107) are used, in which case Eq. (5)reduces to ω κ ( t,s ) g + F ( s ) · D κ ( t,s ) (cid:107) = 0 , (8a)∆ R µ ≡ ∆ r − D k ,µ + D κ ( t,s ) , (cid:107) = 0 , (8b) ω k g + F ( t ) · (cid:104) D κ ( t,s ) , (cid:107) + ∆ r (cid:105) = ω, (8c)with µ = (cid:107) , ⊥ , ⊥ where ˆe ⊥ and ˆe ⊥ are unit vectorsperpendicular to ˆe (cid:107) , and we used β κ ( t,s ) = α κ ( t,s ) , (cid:107) +arg[ F ( s )] in the derivation.
2. Extended recollision model
We solve the saddle-point equations in Eq. (5) semi-classically. The presented method is an extension to the“original” recollision model [2, 37, 38] and will henceforthbe referred to as extended recollision model (ERM). Sincethe bandgaps in semiconductors and insulators are neverzero, we neglect to solve Eq. (5a). Instead, we chooseto consider the electron-hole creation at an initial crystalmomentum, k ≡ κ ( t, s ), taken to be close to or at ahigh-symmetry point. We then integrate the group ve-locities in Eq. (6b) to obtain the classical motions of theelectron and hole, with the time-dependent crystal mo-mentum given by κ ( t, t (cid:48) ) = k + A ( t (cid:48) ) − A ( s ) , t (cid:48) ∈ [ s, t ] . (9)During the propagation, we calculate the “generalized”electron-hole distance vector ∆ R µ in Eq. (5b), and a re-turning trajectory is said to have recollided at time t (cid:48) = t if: (i) (cid:107) ∆ R µ (cid:107) as a function of t (cid:48) is at a local minimumand (ii) the (cid:107) ∆ R µ (cid:107) < R requirement is fulfilled, where R is a preset recollision threshold value. We set R = 30unless indicated otherwise. If a trajectory has recollided,we record s , t and the recollision energy ω ( k , s, t ) inEq. (5c). We initiate trajectories for times spanning anoptical cycle (o.c.) s ∈ [ − T, t ∈ [ s, s + 2 T ]. For each trajectory, we count up to amaximum of 3 recollisions. However, unless indicatedotherwise, we present results for the first recollision. Of-ten, we perform calculations for all k s in a disc withradius ∆ k around a high-symmetry point. We note thatallowing electron tunneling at different k points awayfrom the minimum band gap is inherently distinct fromthe summation over all k points in the expressions forthe total current in Eqs. (2) and (3). In the former case,we attempt to find all the stationary-phase points { k , s , t } that most contribute to the integral (3).We formally define an imperfect recollision as havinga nonzero electron-hole distance (cid:107) ∆ R µ (cid:107) (cid:54) = 0, i.e. when-ever the electron and hole centers do not exactly spatiallyreencounter each other. Note that the electron-hole pairwill get driven further apart spatially whenever the di-rection of the time-dependent crystal momentum [andthus A ( t (cid:48) )] in Eq. (9) is not along the group velocitiesin Eq. (6b), and in such cases the imperfect recollisionscould be important. One can thus control and force suchrecollisions by either tuning the laser or studying differentmaterials, e.g. by using elliptical drivers or consideringmaterials with large Berry curvatures. Also, since thegroup velocities are the gradient of the band dispersions,tunneling at a reciprocal point that is not the minimumband gap can also lead to imperfect recollisions [48], evenwhen using linearly polarized driving fields.It should be mentioned that recent progress [17, 39, 57]has been made towards solving the saddle point equations(5) and performing the stationary phase approximationon the integral in Eq. (3). However, these studies presenta monumental task, even in reduced dimensionalities andwithout the extra terms involving Q κ ( t,s ) and D k ,µ . Italso remains to be seen whether such formalisms can treatelectron-hole-pair creation at different symmetry pointsin the BZ. C. Electron and hole wave-packet analysis
We present here a quantum wave packet method thatis able to provide additional details on the spatially ex-tended nature of the imperfect recollisions, by explicitlyconstructing the real-space electron and hole wave pack-ets for a specific semiclassical trajectory. We label thismethod, which was first applied in our previous work [48]for the electron wave packet, as the wave-packet trajec-tory (WPT) method.For concreteness, consider a semiclassical electron tra-jectory in the conduction band that tunneled at time s and reciprocal coordinate k = K − q A ( s ), where q = − k ( K ) is the crys-tal momentum in the fixed (moving) frame. We expandthe real-space electron wave packet in the Houston-statebasis Ψ e ( r , t ) = (cid:88) K a K e ( t ) h K c ( r , t ) , (10)where (cid:12)(cid:12) a K c ( s ) (cid:12)(cid:12) is chosen to be a Gaussian centered at K , with a full width at half maximum (FWHM) esti-mated by Zener-tunneling as described in Appendix A 3and Eq. (A8). The Houston states [58, 59] are related tothe accelerated Bloch states h K c ( r , t ) = e iq A ( t ) · r φ K − q A ( t ) c ( r ) . (11)Inserting Eq. (10) into the time-dependent Schr¨odingerequation, and neglecting coupling to the other bands,leads to the equations of motion: i ˙ a K e ( t ) = (cid:104) E K − q A ( t ) c − q F ( t ) · A K − q A ( t ) c (cid:105) a K e ( t ) . (12)We thus propagate Eq. (12) starting from time s , andat desired time intervals calculate the real-space wavepacket using Eq. (10). With access to the real-space wavepacket, the observables such as the expectation values (cid:104) r (cid:105) ( t ) and standard deviations σ µ = (cid:113) (cid:104) ˆ µ (cid:105) − (cid:104) ˆ µ (cid:105) , withˆ µ ∈ { ˆ x, ˆ y, ˆ z } , can be calculated. For better visualizationof the width, we define the FWHM-like ¯ σ ≡ √ σ x + σ y ).A hole is left behind in the valence band when an elec-tron tunnels from the valence to the conduction band.Seen as a quasiparticle, the hole has positive charge q h = 1, and satisfies K h = − K (total crystal momentumconservation) and E k h = − E k v . The corresponding equa-tions for the hole wave packet Ψ h ( r , t ) is then obtainedfrom Eqs. (10) and (12) by substituting in the equations a e → a h , h c → h v , K → K h , q → q h and E c → E h .Finally, it should be noted that the construction ofthe quantum wave packet in Eqs. (10) and (11) requiresthe knowledge of the Bloch wave functions φ k n ( r ) whichsometimes can be hard to obtain. III. ELLIPTICITY DEPENDENCY IN A BULKSOLID
Recently, the ellipticity-dependence of HHG in bulksolids has attracted both theoretical and experimentalattention [1, 3, 4, 22, 41, 60–62]. In contrast to HHGin gases, where the HHG yield falls off with increasingellipticity, HHG in solids exhibits nontrivial ellipticitydependence where the harmonic yield can increase withincreasing ellipticity. In this section, we investigate theellipticity-dependence of HHG in a generic bulk-solid sys-tem with the minimum band gap at the Γ point. Weconsider a model for bulk ZnO, using a two-band ap-proximation and neglecting the Berry connections, Berrycurvatures, and TDPs. As we will show below, this treat-ment allows a detailed and quantitative understanding ofthe recolliding trajectories and emission dynamics in ageneric bulk solid.
A. Generic bulk solid: ZnO model
For the band structure of wurtzite ZnO, we considerthe plane containing the Γ, K and M high-symmetrypoints. The band structure is obtained using the analyt-ical model E k n = u − (cid:104) t n (cid:112) f k + q n + t (cid:48) f k + p n (cid:105) , (13a) f k = 2 cos( ak y ) + 4 cos (cid:20) ak y (cid:21) cos (cid:34) √ ak x (cid:35) , (13b)with the fitted parameters t v = 2 . t c = − . q v =4 . q c = 3 . t (cid:48) = − . p v = − . p c = 10 . u = 27 .
1. Our model is adapted from Ref. [41], but nowusing the real lattice constant of a = 6 .
14 for ZnO. The k -dependent band gap shown in Fig. 1(a) is seen to exhibithexagonal symmetry with the minimum band gap at Γ ω Γ g = 3 . d k x = d k y = (cid:115) K ω k g ) (14)with the Kane parameter K = 0 . K ) point acts as a sink with the vec-tors pointing towards it, while the K (Γ) point acts as asource.Note that even though a hexagonal BZ is visualized inFig. 1, in the actual calculations we use a Monkhorst-Pack mesh spanned by the reciprocal vectors b =2 π (3 − ˆe x + ˆe y ) /a and b = 2 π (3 − ˆe x − ˆe y ) /a . FIG. 1. (a) Band gap of ZnO in the plane containing the highsymmetry points Γ, K and M . (b) Norm of the transitiondipole moment. (c) Group velocities of holes in the valenceband. (d) Group velocities of the electrons in the conductionband. The hexagon in the plots guides the eye and traces thefirst BZ. B. Driver ellipticity dependence of HHG in ZnO
We irradiate the bulk with elliptically-polarized vectorpotentials of the form A ( t ) = A g ( t ) √ (cid:15) [sin( ω t ) ˆe x + (cid:15) cos( ω t ) ˆe y ] , (15)where (cid:15) is the ellipticity, ω is the carrier frequency, F = ω A is the electric field maximum, and thepulse envelope is on the form g ( t ) = cos [ πt/ (2 τ )] with t ∈ [ − τ, τ ]. For our calculations in this section, we choose ω = 0 . λ = 3200 nm), A = 0 .
35 and τ = 106 . − M ( ˆe x ) direction.Simulations with the ellipse major axis along Γ − K yieldsnearly indistinguishable results from those in Fig. 2 andwill not be discussed further in this work.Figure 2(a) shows the HHG spectrum for three dif-ferent ellipticities: (cid:15) = 0 (linear polarization), (cid:15) = 0 . (cid:15) = 1 . ω Γ g and harmonic 27 (H27): a drop-off region,a plateau region and a cut-off region. For (cid:15) = 0 .
5, theharmonic intensity in the plateau region is reduced byup to 5 orders of magnitudes compared to (cid:15) = 0, whilein the cut-off region the harmonic yield is actually in-creased going from (cid:15) = 0 to (cid:15) = 0 .
5. The shape of the s p ec t r a li n t e n s i t y ǫ = 0 . . . drop-off plateau cut-off harmonic order0.00.20.40.60.81.0 e lli p t i c i t y -15 -10 -5 FIG. 2. High-harmonic spectra of ZnO driven by ellipticallypolarized light (major axis of ellipse along Γ-M). The laserparameters are ω = 0 . λ = 3200 nm), A = 0 .
35 and τ = 106 . T = 10 fs. The first verticaldashed line at ∼ H8.5 traces the minimum band-gap energy ω Γ g and separates the drop-off region from the plateau region; thesecond vertical line at H27 guides and eye and approximatelyseparates the plateau and cut-off regions. spectrum for circular polarization is qualitatively similarto the (cid:15) = 0 . (cid:15) ∈ [0 , ω Γ g where the interbandharmonics dominate. In the plateau region, a monotonicdecrease of yield with increasing (cid:15) is evident, with theyield almost vanishing at (cid:15) = 1. Such a behavior issimilar to the ellipticity dependence of HHG in gases.The cut-off region in Fig. 2(b), however, exhibits anoma-lous ellipticity-dependence, with relatively large yieldsbetween (cid:15) = 0 and (cid:15) = 1. Qualitatively similar ellipticity-dependencies were reported in Refs. [41, 50]. C. Emission profiles and semiclassical analysis for Γ The character and periodicity of the harmonic time-frequency emission profiles also depend strongly on theellipticity, as illustrated in Fig. 3. For (cid:15) = 0 in Fig. 3(a),the profile exhibits half-cycle periodicity, with the mostprominent feature exhibiting a peak at around H27 andemitted at − .
29 o.c. Overall, it resembles a typical time-frequency profile from HHG in gases, where every energybelow the maximum is emitted twice, corresponding tothe short and long trajectories. In contrast, the time-frequency profile for (cid:15) = 0 . triangular structure, and with theemission time shifted to ∼ − . Γ h a r m o n i c o r d e r Γ h a r m o n i c o r d e r -12 -10 -8 (a) ǫ = 0 A B Γ h a r m o n i c o r d e r Γ h a r m o n i c o r d e r -16 -14 -12 (b) ǫ = 0 . -0.50 -0.25 0.00 0.25 0.50emission time (o.c.) h a r m o n i c o r d e r Γ -0.50 -0.25 0.00 0.25 0.50emission time (o.c.) h a r m o n i c o r d e r -20 -18 -16 -14 (c) ǫ = 1 FIG. 3. Time-frequency emission profiles (colormaps) of theharmonics for three different ellipticities: (a) (cid:15) = 0; (b) (cid:15) = 0 . (cid:15) = 1 .
0. The black lines are semiclassicalERM results with tunneling initiated at k = Γ, counting thefirst two recollisions, and assuming the recollision thresholds R = 30 , ,
100 for (cid:15) = 0 , . ,
1, respectively. Note that semi-classical recollisions with tunneling from the Γ point are onlyobserved for the linearly polarized driver in (a). A and Bmark different classes of trajectories as discussed in the text. ted with energies corresponding to the plateau region (H9to H21) are shifted in time by a quarter cycle comparedto the triangular structure. The time-frequency profilefor (cid:15) = 1 in Fig. 3 shows six burst of light during eacho.c. – a clear reflection of the six-fold rotational symme-try of the BZ (see Fig. 1).We first analyze the emission profiles for the linearly-polarized case, by using the ERM in Sec. II B and as-suming that tunneling occurs at the minimum band gapΓ. The emission times for individual trajectories areshown in Fig. 3(a) by the gray dots. The agreementwith the colormap is quite good, with the semiclassicalresults reproducing the emission profiles during each half-cycle. Even the peculiar structure at ∼ H38 is capturedby the semiclassical model. The very different emissionprofiles of the trajectories labeled A and B in Fig. 3(a)suggest that they belong to different classes of trajecto- − M Γ M -0.30.00.3 elec.hole030609010203040 -0.9 -0.5 -0.1 -0.9 -0.5 -0.1 κ x ( a . u . ) Γ traj. A − M Γ M (a) Γ traj. B(b) v x ( a . u . ) -0.30.00.3 (c) (d) k ∆ r k ( a . u . ) (e) (f) e n e r g y ( ω ) time (o.c.) (g) time (o.c.) -0.9 -0.5 -0.1 (h) FIG. 4. Two different semiclassical trajectories for the lin-early polarized driver with Γ as the tunnel point. The rec-ollision event is marked with the filled circle in each panel.Left panels: trajectory that (perfectly) recollide at t = − . ω = 27 ω ; right panels: trajectorythat (imperfectly) recollide at t = − .
14 o.c. with recollisionenergy ω = 38 ω . The gray lines in (a) and (b) show A x ( t ).In (g) and (h), the gray lines show the energies without the eh -PER, while the purple lines show the total energies. ries. This is further illustrated in Fig. 4, where we in theleft panels consider trajectory A. Figure 4(a) shows thereciprocal-space motion: the electron-hole pair is createdat the Γ point at time s = − .
941 o.c., and afterwardis driven by the vector potential according to Eq. (9).The time-dependent crystal momentum initially movestoward − M , and later changes direction when the vec-tor potential changes direction; it never moves beyondthe BZ boundaries (at ± M ). In Fig. 4(c), the elec-tron group velocity is negative (positive) when k x < k x > (cid:107) ∆ r (cid:107) = 0) attime − .
29 o.c. in Fig. 4(e), and consequently the recol-lision energy in Fig. 4(g) with and without the eh -PERis the same. Trajectories of class A in Fig. 3(a) are thussimilar to the ones in HHG in gases, consisting of shortand long trajectories.Consider now the special trajectory labelled B inFig. 3(a). After tunneling at Γ, the crystal momentumgoes beyond the BZ-boundary [Fig. 4(b)], where electronand hole group velocities abruptly change sign [Fig. 4(d)], and undergo a Bragg reflection. Consequently, the elec-tron and hole only imperfectly recollide in real space,with (cid:107) ∆ r (cid:107) ≈
30 shown in Fig. 4(f). The resulting ex-tra eh -PER contributes ∼ ω which is added to thetotal energy of the emitted harmonics in Fig. (4)(h).Bragg reflections can thus lead to imperfect recollisionsin bulk solids even when using linearly-polarized drivers.Note that the effects of the Bragg reflection on HHG insolids have been investigated in several previous works[1, 41, 63–65].In contrast to the linear-polarization case, the ellipti-cally polarized fields do not lead to any recollisions initi-ated from the Γ point, as shown in Figs. 3(b) and 3(c).To ensure that this is not just due to a larger recollisiondistance, the recollision thresholds has been relaxed from R = 30 at (cid:15) = 0 to R = 60 at (cid:15) = 0 . R = 100 at (cid:15) = 1. Including only trajectories initiated at the mini-mum band is thus insufficient for the description of HHGwith elliptically polarized drivers in bulk solids. D. Full semiclassical picture - effect of the full BZ
We now extend our semiclassical ERM analysis to in-clude tunneling from a disc around Γ in reciprocal space,as described previously in Sec. II B. In our calculations,we choose the radius of the disc for (cid:15) = 0 , . , k = 0 . , . , .
25, respectively. In the left (right) pan-els of Fig. 5, the results for the ERM simulations without(with) the eh -PER are shown together with the quantumresults in the background. For (cid:15) = 0 in Figs. 5(a) and5(b), the ERM result is similar to the one in Fig. 3(a):the overall structure is broader, and the short and longtype of trajectories are now continuously connected tothe higher energy structure. The result with and with-out the eh -PER are also similar, with Fig. 5(b) havingsome trajectories with higher energy, forming a “boot”structure. We stress here that the goal of our analysis ispurely to reveal all possible semiclassical recolliding tra-jectories, and we currently do not weight each trajectorywith the tunneling and recollision probabilities.For (cid:15) = 0 . eh -PER. It is importantto notice that the triangular structure is not due to tra-jectories tunneled at a single k , but rather trajectoriesfrom different k points that collectively give rise to thefull triangular emission structure. This new finding is instark contrast to HHG in gases and updates the previ-ous conception for HHG in solids where the trajectoriestunnelled from the minimum band gap is the only onesthat mattered. Interestingly, comparing the (cid:15) = 0 casewith (cid:15) = 0 .
5, the origin of the triangular structure seemsto be due to the class B trajectories. no EH-PER with EH-PER h a r m o n i c o r d e r h a r m o n i c o r d e r (a) ǫ = 0 (b) ǫ = 0 h a r m o n i c o r d e r h a r m o n i c o r d e r (c) ǫ = 0 . (d) ǫ = 0 . -0.50 -0.25 0.00 0.25emission time (o.c.) h a r m o n i c o r d e r -0.50 -0.25 0.00 0.25emission time (o.c.) h a r m o n i c o r d e r (e) ǫ = 1 -0.50 -0.25 0.00 0.25emission time (o.c.)-0.50 -0.25 0.00 0.25emission time (o.c.)(f) ǫ = 1 FIG. 5. Semiclassical recollision energies versus recollisiontimes obtained using the ERM (gray dots) with tunnelinginitiated for all k in a disc around Γ. Left (right) panelsshow the semiclassical results without (with) the inclusion ofthe eh -PER F ( t ) · ∆ r in Eq. (5c), while different row panelscorrespond to different ellipticities (cid:15) . For (cid:15) = 0 , . ,
1, wechose ∆ k = 0 . , . , .
25 and R = 30 , , For the circularly polarized case in Figs. 5(e) and 5(f),the six-fold symmetry is clearly reproduced in the semi-classical calculations. Inclusion of the eh -PER repro-duces the almost vertical structures extending up to H60.For clarity, only the first recollision for each trajectory iscounted in Fig. 5; when more recollisions are counted,the ERM reproduces more features in the colorplots (seeFig. 13 in Appendix B).To explore the contributions from different initial tun-nel sites k to the time-frequency profiles and the HHGspectra, we show in Fig. 6 the maximum recollision en-ergy (colorbar) as a function of the tunnel site in recip-rocal space. Each subfigure in Fig. 6 uses the same ERMdata set as the corresponding subfigure in Fig. 5. Forexample, the ERM calculation in Fig. 5(a) contains allinitial k points inside the gray circle in Fig. 6(a) (withradius ∆ k = 0 . k . For the (cid:15) = 0 case shown in Figs. 6(a) and6(b), the recolliding trajectories clearly originate with k along the k x -axis (laser polarization direction): electron- no EH-PER with EH-PER -0.20.00.2 k y ( a . u . ) -0.20.00.2 k y ( a . u . ) h a r m o n i c o r d e r (a) ǫ = 0 h a r m o n i c o r d e r (b) ǫ = 0 -0.20.00.2 k y ( a . u . ) -0.20.00.2 k y ( a . u . ) (c) ǫ = 0 . (d) ǫ = 0 . -0.2 0.0 0.2 k x (a.u.) -0.20.00.2 k y ( a . u . ) -0.2 0.0 0.2 k x (a.u.) -0.20.00.2 k y ( a . u . ) (e) ǫ = 1 -0.2 0.0 0.2 k x (a.u.)-0.2 0.0 0.2 k x (a.u.)(f) ǫ = 1 FIG. 6. Maximum semiclassical recollision energy for differentinitial tunneling crystal momenta k in the BZ. Left (right)panels show the semiclassical results without (with) the in-clusion of the eh -PER F ( t ) · ∆ r in Eq. (5c), while differentrow panels correspond to different ellipticities (cid:15) . All k pointscontained inside the gray circles are propagated, and missingpoints indicate no trajectory recollision (or recollision energy ω < ω ). The k with the highest recollision energy in (d)is enclosed by a square. hole trajectories created too far away from the k x -axiswill be driven apart in the y -direction in real space andnever recollide [see Figs. 1(c) and 1(d)]. For (cid:15) = 0 . (cid:107) k (cid:107) (cid:38) .
15 can recollide, and the six-fold symmetry of the BZ is clearly visible. Note in Fig. 6that the larger the ellipticity, the larger the “hole” aroundΓ becomes, and the less the trajectories starting near Γcontribute to the interband emissions. The maximumrecollision energies including the eh -PER (right panelsin Fig. 6) are substantially higher than the calculationswithout (left panels), in agreement with the results inFig. 5. Figure 6 again reinforces our central finding thattaking into account only k = Γ is insufficient for thedescription of HHG in bulk solids with elliptically polar-ized drivers. While tunneling indeed occurs mostly at r ec o l. e n e r g y ( ω ) recol./tunnel time (o.c.)tunnelrecol. (a) k y ( a . u . ) k x (a.u.)recol. -0.50.00.5 -0.5 0.0 0.5 (b) k ∆ r k ( a . u . ) time (o.c.) (c) e n e r g y ( ω ) time (o.c.) (d) FIG. 7. Semiclassical analysis for (cid:15) = 0 . k = (0 . , . s = − . t = − .
081 o.c.): (b) Motion in recipro-cal space; (c) Electron-hole distance versus time; (d) Energyversus time, with the gray (purple) line showing the energieswithout (with) the eh -PER. Γ [e.g. the band gap at k = (0 , , .
1) is increased by13% compared to k = Γ], the dynamics imposed by thelaser and the dispersion relation is such that recollisionis prevented.For (cid:15) = 0 .
5, the k with the highest recollision energyis marked with a square in Fig. 6(d). Correspondingly,this k gives rise to the tip of the triangular structure intime-frequency profiles of Fig. 5(d). Figure 7(a) showsthe recollision energies versus the tunnel and recollisiontimes for all trajectories originating with this k . Clearly,all resemblances to the short and long trajectories fromgas-phase HHG are gone. Instead, the recollision ener-gies versus the recollision times exihibits a highly irregu-lar structure, with harmonics above order ∼
35 emittedapproximately at the same time.To give an example, we now focus on the trajec-tory with the highest recollision energy that tunnels at s = − .
83 o.c. and recollides at t = − .
081 o.c. Thetime-dependent crystal momentum shown in Fig. 7(b)extends beyond the first BZ, and the electron-hole tra-jectories are seen to recollide imperfectly with (cid:107) ∆ r (cid:107) ≈ ∼ eh -PER [Fig. 7(d)], whichleads to the correct reproduction of the triangular struc-ture in Fig. 5(d) and not in Fig. 5(c).To summarize this section, we have shown that for ageneric bulk solid with the minimum band gap at Γ, ellip- -0.50.00.5 -0.5 0.0 0.5 -0.5 0.0 0.5-0.50.00.5 k y ( a . u . ) ω k g (eV) Γ KM M M K ′ (a) band gap k d k k (a.u.) Γ KM M M K ′ (b) dipole coupl. k y ( a . u . ) k x (a.u.)(c) hole velocity k x (a.u.)(d) electron velocity FIG. 8. (a) Band gap of monolayer hBN. (b) Norm of thetransition dipole moment. Group velocities of (c) holes in thevalence band and (d) of the electrons in the conduction band,without the anomalous velocity. The hexagon in the plotsguides the eye and traces the first BZ. tical drivers enhance the harmonic emissions at high fre-quencies typically associated with the cut-off region of aharmonic spectrum, and greatly reduce the harmonic in-tensity in the plateau region. The time-frequency analy-sis reveals that the highest-order harmonics are not emit-ted from trajectories tunnelled at Γ, but rather due tocollective emissions originating from many k s near Γ. IV. RECOLLISIONS IN A MONOLAYERBANDGAP MATERIAL
In the previous section, we investigated HHG in ageneric model for bulk solids where the minimum bandgap is at the high-symmetry point Γ of the BZ. Inbandgap monolayer matarials, in contrast, the minimumband gap is usually located at the high-symmetry point K , with the maximum band gap at Γ. In this section, weinvestigate the wavelength and orientation dependence ofHHG in a typical topologically trivial monolayer system,using the formalisms presented in Sec. II. A. Typical monolayer system: hBN model
We use monolayer hBN as an example of a typicalmonolayer band-gap material. For the band structurecalculations, we employ the pseudo potential from [66],0 M KM K ′ -11 -9 -7 (a) harmonic order(a) -0.25 0.00 0.25emission time (o.c.)20304050 h a r m o n i c o r d e r -16 -14 -12 (b) 0 ◦ ( Γ − M ) FIG. 9. (a) HHG spectrum for the parallel-polarized har-monics as a function of the driver polarization orientation.The pulse parameter of the driver is ω = 0 . λ = 2 . µ m), F = 0 .
010 ( A = 0 . ∼ H
15 indicates the minimum band gap. (a) Time-frequency profile for the Γ − M driver polarization direction,with the semiclassical ERM result superimposed. For theERM results, the gray points are recollisions with tunnelingaround k = M with ∆ k = 0 . R = 15; the purple points are recollisions with tunnelingaround k = M , M , K, K (cid:48) . and we employ the twisted parallel transport gauge [67]to obtain BZ-periodic transition dipole moments andBerry connections [68]. Contrary to the bulk case, theband gap is smallest near the K and M symmetry pointsand largest at Γ, as shown in Fig. 8(a). Correspondingly,the norm of the dipole coupling is largest near the K and M points in Fig. 8(b). Figures 8(c) and 8(d) show that Γ( K ) acts as a source (sink) and K acts as a sink (source)for the hole (electron) group velocity. B. Orientation dependence of HHG in hBN
We irradiate hBN with linearly polarized infraredpulses and investigate the HHG process with respect tothe driver polarization angle θ . The chosen field param-eters are λ = 2 . µ m, F = 0 .
010 and τ contains 5.5o.c. The SBEs are solved with T = 10 fs, and the HHGspectra for the parallel-polarized harmonics are shown inFig. 9(a). The six-fold symmetry of the BZ is clearlyreflected in the spectrum, with stronger yields along theΓ − K directions compared to the Γ − M directions.Figure 9(b) shows the time-frequency profile for theΓ − M driver direction with the semiclassical ERM re-sult superimposed [69]. Part of the time profile resem-bles that due to the short trajectories in HHG in gases,with a single “arm” extending from H20 to H50 dur-ing each half cycle. In addition, the most intense partof the radiation is emitted between H15 and H25, attimes around t = 0 , ± . k = M with ∆ k = 0 . R = 15. Thesingle arm in the time-frequency profile is clearly repro-duced. Note that the current situation is similar to the atomic HHG case, as well as HHG in bulk solids drivenby linearly polarized pulses described in Sec. III. In allthese cases, the group velocities [Figs. 8(c) and 8(d)] ofthe trajectories are pointing along the vector potential A ( t ), leading to (almost) perfect recollisions. The intensefeatures at lower-order harmonics are due to recollisionsfrom trajectories that tunnel near the other symmetrypoints k = M , M , K, K (cid:48) . The ERM results originat-ing from these points are shown in Fig. 9(b) by the purplepoints and reproduce the intense features very well. Thefact that trajectories originating from M lead to muchhigher recollision energies compared to the other sym-metry points can be intuitively predicted by consideringthe band structure in Fig. 8(a): starting from the M symmetry point, the time-dependent crystal momentum κ moving along the Γ − M direction can get closer tothe large-band-gap region near the Γ point, compared toif one starts from a non- M symmetry point. Our resultshere show that tunneling from different regions in the BZcan lead to distinct regions in the emission profiles sep-arated in frequency and time. In such cases, the ERMprovides a full understanding of the emission dynamics.The case for θ = 30 ◦ , i.e. driver polarization along Γ − K , will be discussed in detail in the next subsection. Wehere only note that compared to the θ = 0 ◦ case discussedhere, additional complexities will arise by considering thedispersion relations in Figs. 8(c) and 8(d): the groupvelocities of the electron-hole pairs that start near k = M will no longer be along the vector potential direction(Γ − K ). C. Wavelength dependence of HHG in hBN
We irradiate hBN with linearly polarized pulses alongthe Γ − K direction, keeping the field maximum fixedat F = 0 .
010 and varying the wavelength λ from 1.6 µ m to 2.4 µ m, with the FWHM of the pulse τ chosento contain 5.5 o.c.. The HHG spectra calculated fromthe SBEs are shown in Fig. 10(a) and Fig. 10(b) for theparallel- and perpendicular-polarized harmonics, respec-tively. The HHG spectra extend toward higher harmonicenergies with increasing wavelengths, which can be qual-itatively understood simply by the larger A and conse-quently the larger excursion of the time-dependent crys-tal momenta in Eq. (9).The time-frequency profiles in the colorplots of Fig. 11reveal the emission dynamics of the HHG process. Forthe 1.6 µ m case in Fig. 11(a) the characteristic double-peak structure during each half-cycle is observed, whichwas studied in detail in Ref. [48]. When the wave-length is increased, during each half-cycle, the double-peak seems to split into two almost-vertical, downwards-sloping structures, as shown in Figs. 11(c). The ERMresults in the case of only taking into account the M and M symmetry points are overlaid on top of the col-orplots in Figs. 11(a) and 11(c). Clearly, in Fig. 11(a)they are unable to reproduce the double-peak structure1 s p ec t r a li n t e n s i t y µ m2.0 µ m2.4 µ m (a) k HH s p ec t r a li n t e n s i t y harmonic energy (eV) (b) ⊥ HH FIG. 10. High-harmonic spectra of hBN driven by linearlypolarized pulse along Γ-K, with (a) parallel harmonics and(b) perpendicular harmonics. The field amplitude is fixed at F = 0 . T = 10 fs. The vertical dashed linestrace the minimum band-gap energy. and in Fig. 11(c) the semiclassical results seem to be atodds with the colorplot, predicting recollisions at timeswhen there are actually no emissions.We extend the ERM analysis to include all tunnelpoints k in a disc of radius ∆ k = 0 . M and M , shown by the gray dots in Figs. 11(b) and 11(d).The recollision threshold is set to R = 15. For the 1600nm case in Figs. 11(b), it is seen that the double peakstructures are attributed to imperfect recollisions for tra-jectories tunnelled close to the M points [48]. For the 2 . µ m case in Fig. 11(d), we first focus on the downwards-sloping structures enclosed by the rectangular boxes inFig. 11(d). The ERM results are seen to reproduce thesedownwards-sloping structures quite well, although theyare not continuous, and resemble groups of horizontallines separated by vertical spacings. These gaps are dueto the density of discrete k points chosen in our simula-tions: when the density is increased, the empty spacingsin the semiclassical results get filled.We now turn our attention to the prominentdownwards-sloping structures in the time profiles thatare not reproduced, highlighted by the rectangular boxesin Figs. 11(e) and 11(f). We argue that these structuresare due to not just contributions from a number of dif-ferent regions in the BZ, as we have seen above, butalso how these contributions interfere with each other.We perform ERM calculations with R = 30 and tak-ing into account the first three recollisions (instead ofone). The gray dots in Fig. 11(f) show the ERM resultsfor trajectories with k near M and M . The prominentdownward-sloping structures in the boxes are mostly cov-ered by the ERM results. However, the time-frequencyprofile contains prominent holes in the left bottom part -15 -13 -11 (a) 1.6 µ m -15 -13 -11 (b) 1.6 µ m h a r m o n i ce n e r g y ( e V ) h a r m o n i ce n e r g y ( e V ) (c) 2.4 µ m (d) 2.4 µ m-0.50 -0.25 0.00 0.25emission time (o.c.) -0.50 -0.25 0.00 0.25emission time (o.c.) (e) 2.4 µ m -0.50 -0.25 0.00 0.25emission time (o.c.)-0.50 -0.25 0.00 0.25emission time (o.c.)(f) 2.4 µ m FIG. 11. Time-frequency profiles (colormaps) for the parallel-polarized harmonics obtained from the SBE simulations withthe semiclassical ERM results superimposed. Panels (a) and(b) are for a 1 . µ m driver, while panels (c)-(f) are for a 2 . µ m driver. For the ERM results, panels (a) and (c) showthe results for tunneling initiated only at k = M , M , whilepanels (b) and (d) are for tunneling in discs of radius ∆ k =0 . M , M symmetry points, with the recollisionthreshold set to R = 15. Panels (e) and (f) show ERMresults with the first recollisions and R = 30: (e) resultsfor tunneling from discs around k = M , K, K (cid:48) , while (f)results for tunneling from discs around k = M , M . Therectangular squares are guides to the eye for the discussionsin the text. of the boxes, indicating the absence of harmonic emis-sions, which are not reproduced by the ERM results. Toexamine further, we show in Fig. 11(e) the ERM resultsfor recolliding trajectories initiated with k near the K , K (cid:48) and M symmetry points. Recollisions are observedcovering the lower left parts of the boxes, exactly in theregions where emissions should be absent according tothe quantum results. Since emissions with the same timeand harmonic energy should be added coherently, andthe trajectories inititated near all the symmetry points( M , M , M , K , K (cid:48) ) overlap here with widely differentphases, they appear to destructively interfere and leadto the absence of emissions. The dominant structures inthe time profiles are then reproduced by the parts of theERM result in Fig. 11(f) that do not overlap with theERM result in Fig. 11(e). Note that a definite proof ofthe described destructive interference effect is beyond theERM and the scope of the current work. Still, the semi-classical method gives us insight on where in the BZ thedifferent trajectories originate, which can lead to, in ouropinion, a satisfying understanding of this interference2 Hole Electron semiclassicalwave packet-20-1001020 y ( a . u . ) (a) 1600 nm (b) 1600 nm-20 -10 0 10 20 x (a.u.)-20-1001020 y ( a . u . ) (c) 2400 nm -20 -10 0 10 20 x (a.u.) (d) 2400 nm ¯ σ ( a . u . ) time (o.c.) ¯ σ ( a . u . ) time (o.c.) ¯ σ ( a . u . ) time (o.c.) ¯ σ ( a . u . ) time (o.c.) FIG. 12. Quantum wave packet results for specific recollidingtrajectories that tunnel at k = M with initial FWHM 0 . µ m driver pulse and the specific trajectory that tunnelsat s = − .
931 o.c. and recollides at t = − .
321 o.c. The solidblack curves show the trajectory real-space motion obtainedusing ERM, while the dashed red curve show the expectationvalue (cid:104) r (cid:105) of the wave packets. The real-space position atthe time of recollision t is indicated by the circle, with thecolorplot showing the wave packet probability density at t .The insets show the widths of the wavepackets from s to t .(c),(d) Same as (a) and (b), but for the 2.4 µ m driver and thetrajectory that tunnels at s = − .
922 o.c. and recollides at t = − . effect and the final dynamics.Again, in this section we have found that the noveltime-frequency profiles for HHG in solids are due to thecollective emission of harmonics originating from differ-ent k points in the BZ. D. Quantum wave packet analysis
We have shown that the semiclassical ERM model isable to capture the emission dynamics of the HHG pro-cess in solids. The imperfect recollision and the originof the eh -PER can be interpreted in the context of spa-tially extended wave packets at the time of recollision.We employ the formalism described in Sec. II C to con-struct and visualize such wave packets. We assume tun-neling at k = M with initial FWHM 0 . µ m driver polarized along Γ − K , we consider the tra-jectory with the highest recollision energy which tunnelsat s = − .
931 o.c. and recollides at t = − .
321 o.c. The dashed red curves in Figs. 12(a) and 12(b) show thereal-space motion (cid:104) r (cid:105) of the hole and electron wave pack-ets, respectively. The quantum wave packet results areseen to agree perfectly with the semiclassical ERM re-sults shown by the solid black curves. The hollow circlesindicate the real-space position at the time of recollision t , and the colorplots show the wave packet probabilitydensity at t . The wave packets have a large width andextend over many lattice sites. The time dependence ofthe wave packet width ¯ σ from s to t is shown in the insetsof Figs. 12(a) and 12(b). Due to the lower effective massof the conduction band, the electron wave packet moves agreater distance compared to the hole and spreads more.Note that at the time of recollision, the electron and holewavepackets clearly occupy the same spatial region andoverlap.Figures 12(c) and 12(d) show the wave packet resultsfor the 2.4 µ m driver and the trajectory that attains thehighest recollision energy. The quantum wave packet andthe semiclassical motion are again in full agreement. Dueto the longer half-cycle and larger A compared to the 1.6 µ m case, the electron and hole are driven apart furtheralong the y -direction according to the group velocities inFig. 8 before the vector potential changes sign, leadingto a larger recollision distance. The large electron-holespatial separation at the time of recollision, however, doesnot prevent their spatial overlap, as evidenced by thethe wave packet densities. Due to longer time durationbetween tunneling and recollision for the 2.4 µ m, thewave packets spread more compared to the 1.6 µ m cases(insets of Fig. 12). V. CONCLUSION AND OUTLOOK
We have presented a recollision formalism for HHGin solids that conclusively shows that in many realisticsituations the harmonic spectrum and emissions are notdue to tunneling at the minimum band gap, but insteaddue to the collective effect of trajectories originating neardifferent symmetry points in the BZ. Indeed, for the ex-ample of HHG in a bulk solid with elliptical drivers, weshowed that the electron-hole pairs created at Γ do notrecollide at all and contribute nothing to the highest-order harmonics. For monolayer materials with hexag-onal symmetry, the highest order harmonics originatenot from the minimum band gap at the K symmetrypoints, but near the M points. In addition, we foundthat the HHG for different driver orientations results invery distinct time-frequency profiles, and we showed thatthis is due to collective emissions from many differentreciprocal-space tunneling sites. Interestingly, for cer-tain driver orientations, different parts of the emissionprofiles can be ascribed to electrons initially tunnelingnear different symmetry points in the BZ, allowing forfuture prospects of probing the BZ tunneling regions. Wealso showed that the imperfect recollisions leading to theelectron-hole polarization energies are important for the3correct description of the harmonic emissions, a resultwhich is further supported by our quantum wavepacketconstructions. Situations where solid-state HHG differssignificantly from gas-phase HHG can be summerized bytwo simple rules of thumb: (A) whenever the instanta-neous carrier group velocities are not along the laser po-larization direction; (B) when the time-dependent crystalmomentum goes beyond the BZ boundaries and inducesBragg reflections.Our work illustrates the complexity of HHG in solidscompared to HHG in the gas phase and broadens the no-tion of which parts of the BZ contribute to the emissionprocess - in particular that the most important symme-try point is not always at the minimum band gap. Thedetailed knowledge gained from the collective emissionsresponsible for the novel time-frequency profiles, asidefrom the fundamental perspective, will have impact onfuture experiments that involves phase matching and ul-trafast spectroscopy. The strong interest in the genera-tion of elliptically-polarized harmonics will also benefitfrom this work. Furthermore, the fact that tunnelingfrom different regions in the BZ leads to distinct har-monic time-frequency characteristics can potentially fa-cilitate the all-optical reconstruction of the band struc-ture not only near the minimum band gap as demon-strated in [16], but near all relevant symmetry points inthe BZ. ACKNOWLEDGMENTS
Appendix A: Derivations
In this appendix, we provide some more details onsome of the derivation steps in Sec. II.
1. Saddle-point equations
In the two-band approximation, the SBEs in Eq. (1)reduces to˙ ρ K vv = i F · d K + A ρ K vc + c.c. (A1a)˙ ρ K cc = − i F · d K + A ρ K vc + c.c. (A1b)˙ ρ K cv = (cid:104) − iω K + A g − i F · ∆ A K + A − T − (cid:105) ρ K cv − i (cid:0) ρ K vv − ρ K cc (cid:1) F · d K + A , (A1c)where ∆ A k ≡ A k c − A k v and for notational conveniencethe explicit time-dependencies in A ( t ), F ( t ) and ρ K mn ( t ) have been omitted. Now we make the approximation ofminimum population transfer for the conduction band, ρ K vv − ρ K cc ≈
1, the formal solutions to Eq. (1) read ρ K vv ( t ) = i (cid:90) t ds F ( s ) · d K + A ( s ) ρ K vc ( s ) + c.c. (A2a) ρ K cc ( t ) = − i (cid:90) t ds F ( s ) · d K + A ( s ) ρ K vc ( s ) + c.c. (A2b) ρ K cv ( t ) = − i (cid:90) t ds F ( s ) · d K + A ( s ) e − T − ( t − s ) (A2c) × e − i (cid:82) ts (cid:104) ω K + A ( t (cid:48) ) g + F ( t (cid:48) ) · ∆ A K + A ( t (cid:48) ) (cid:105) dt (cid:48) . (A2d)The expression for the interband current in Eq. (3)is then obtained by inserting the solution (A2) intoEq. (2a), and transforming into the fixed frame k = K + A ( t ).The saddle point conditions for the interband har-monics are obtained by taking the partial derivatives of S µ ( k , t, s ) − ωt . The derivative with respect to s reads(henceforth the dephasing time T is ignored) ∂ s [ S µ ( k , t, s ) − ωt ]= − ω κ ( t,s ) g − F ( s ) · ∆ A κ ( t,s ) − ∂ s β κ ( t,s ) = − ω κ ( t,s ) g − F ( s ) · Q κ ( t,s ) , (A3)where we used ∂ s β κ ( t,s ) = − F ( s ) · ∇ k β k | k = κ ( t,s ) , and Q κ ( t,s ) is defined in Eq. (7b). The derivative with respectto k reads ∇ k [ S µ ( k , t, s ) − ωt ]= (cid:90) ts (cid:110) ∇ k ω κ ( t,t (cid:48) ) g + ∇ k (cid:104) F ( t (cid:48) ) · ∆ A κ ( t,t (cid:48) ) (cid:105)(cid:111) dt (cid:48) + ∇ k α k ,µ − ∇ k β κ ( t,s ) = (cid:90) ts (cid:110) v κ ( t,t (cid:48) ) c − v κ ( t,t (cid:48) ) v (cid:111) dt (cid:48) − ∆ A κ ( t,t (cid:48) ) (cid:12)(cid:12)(cid:12) t (cid:48) = tt (cid:48) = s + ∇ k α k ,µ − ∇ k β κ ( t,s ) =∆ r − D k ,µ + Q κ ( t,s ) , (A4)with D k ,µ , ∆ r and v κ ( t,t (cid:48) ) n defined in Eqs. (5) and (6),and in the second equality we used the identity ∇ ( A · B ) = ( A · ∇ ) B + ( B · ∇ ) A + A × ( ∇ × B ) + B × ( ∇ × A ),such that (cid:90) ts ∇ k (cid:104) F ( t (cid:48) ) · A κ ( t,t (cid:48) ) n (cid:105) dt (cid:48) = (cid:90) ts (cid:110) [ F ( t (cid:48) ) · ∇ k ] A κ ( t,t (cid:48) ) n + F ( t (cid:48) ) × (cid:104) ∇ k × A κ ( t,t (cid:48) ) n (cid:105)(cid:111) dt (cid:48) = (cid:90) ts (cid:104) − ∂ t (cid:48) A κ ( t,t (cid:48) ) n + F ( t (cid:48) ) × Ω κ ( t,t (cid:48) ) n (cid:105) dt (cid:48) = − A κ ( t,t (cid:48) ) n (cid:12)(cid:12)(cid:12) t (cid:48) = tt (cid:48) = s + (cid:90) ts (cid:104) F ( t (cid:48) ) × Ω κ ( t,t (cid:48) ) n (cid:105) dt (cid:48) . (A5)4The derivative with respect to t reads ∂ t [ S µ ( k , t, s ) − ωt ]= ω k g + F ( t ) · ∆ A k − ∇ κ β κ ( t,s ) · F ( t ) − ω + (cid:90) ts ∂ t (cid:104) ω κ ( t,t (cid:48) ) g + F ( t (cid:48) ) · ∆ A κ ( t,t (cid:48) ) (cid:105) dt (cid:48) = ω k g + F ( t ) · (cid:104) ∆ A k − ∇ κ β κ ( t,s ) (cid:105) − ω + (cid:90) ts ∇ κ (cid:104) ω κ ( t,t (cid:48) ) g + F ( t (cid:48) ) · ∆ A κ ( t,t (cid:48) ) (cid:105) dt (cid:48) · F ( t )= ω k g + F ( t ) · (cid:104) ∆ A k − ∇ κ β κ ( t,s ) (cid:105) − ω + (cid:90) ts ∇ κ ω κ ( t,t (cid:48) ) g dt (cid:48) · F ( t ) − ∆ A κ ( t,t (cid:48) ) (cid:12)(cid:12)(cid:12) t (cid:48) = tt (cid:48) = s · F ( t )+ (cid:90) ts (cid:104) F ( t (cid:48) ) × (cid:16) Ω κ ( t,t (cid:48) ) c − Ω κ ( t,t (cid:48) ) v (cid:17)(cid:105) dt (cid:48) · F ( t )= ω k g + F ( t ) · (cid:16) ∆ A κ ( t,s ) − ∇ κ β κ ( t,s ) (cid:17) − ω + (cid:90) ts (cid:16) v κ ( t,t (cid:48) )2 − v κ ( t,t (cid:48) )1 (cid:17) dt (cid:48) · F ( t )= ω k g + F ( t ) · (cid:104) Q κ ( t,s ) + ∆ r (cid:105) − ω. (A6)The saddle-point conditions in Eq. (5) are then obtainedby setting Eqs. (A3), (A4) and (A6) to zero.
2. Structure-gauge invariance of D k ,µ and Q k Under the gauge transformation (cid:12)(cid:12) u k n (cid:11) → (cid:12)(cid:12) u k n (cid:11) e iϕ k n ,with n = v, c , the relevant quantities transform as A k n → A k n − ∇ k ϕ k n (A7a) d k → d k e − i ( ϕ k c − ϕ k v ) (A7b) α k ,µ → α k ,µ − ϕ k c + ϕ k v (A7c) β k → β k − ϕ k c + ϕ k v (A7d) D k ,µ → D k ,µ (A7e) Q k → Q k , (A7f)where the transforms in Eqs. (A7e) and (A7f) are ob-tained by using the definitions in Eq. (7).
3. Approximation of tunneling width in WPT
The Landau-Zener tunneling probability [70–72] reads P k ∝ exp (cid:34) − πω k g | E · d k | (cid:35) , (A8) with ω the laser carrier frequency and E chosen at a timewhen (cid:12)(cid:12) E · d k (cid:12)(cid:12) is maximal. In our WPT calculations, westart with a Gaussian wave packet in reciprocal space,with the FWHM width approximated by the FWHM ofthe above formula. ǫ = 0 ǫ = 0 . h a r m o n i c o r d e r h a r m o n i c o r d e r (a) max 1 recol. (b) max 1 recol. h a r m o n i c o r d e r h a r m o n i c o r d e r (c) max 2 recol. (d) max 2 recol.-0.50 -0.25 0.00 0.25emission time (o.c.) h a r m o n i c o r d e r -0.50 -0.25 0.00 0.25emission time (o.c.) h a r m o n i c o r d e r (e) max 3 recol. -0.50 -0.25 0.00 0.25emission time (o.c.)-0.50 -0.25 0.00 0.25emission time (o.c.)(f) max 3 recol. FIG. 13. Semiclassical recollisions energies versus recollisiontimes obtained with the ERM (gray dots). Left (right) panelsare for (cid:15) = 0 ( (cid:15) = 0 . Appendix B: Supplemental calculations [1] S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini,L. F. DiMauro, and D. A. Reis, “Observation of high- order harmonic generation in a bulk crystal,” Nat. Phys. , 138 (2011).[2] G. Vampa, T. J. Hammond, N. Thir´e, B. E. Schmidt,F. L´egar´e, C. R. McDonald, T. Brabec, and P. B.Corkum, “Linking high harmonics from gases and solids,”Nature , 462 (2015).[3] Y. S. You, D. A. Reis, and S. Ghimire, “Anisotropichigh-harmonic generation in bulk crystals,” Nat. Phys. , 345 (2017).[4] G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne,K. J. Schafer, M. B. Gaarde, and D. A. Reis, “Solid-state harmonics beyond the atomic limit,” Nature ,520 (2016).[5] M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Kloster-mann, A. Guggenmos, and E. Goulielmakis, “Multi-petahertz electronic metrology,” Nature , 359 (2016).[6] Z. Wang, H. Park, Y. H. Lai, J. Xu, C. I. Blaga, F. Yang,P. Agostini, and L. F. DiMauro, “The roles of photo-carrier doping and driving wavelength in high harmonicgeneration from a semiconductor,” Nat. Commun. ,1686 (2017).[7] S. Ghimire, G. Ndabashimiye, A. D. DiChiara,E. Sistrunk, M. I. Stockman, P. Agostini, L. F. DiMauro,and D. A. Reis, “Strong-field and attosecond physics insolids,” J. Phys. B , 204030 (2014).[8] S. Y. Kruchinin, F. Krausz, and V. S. Yakovlev, “Col-loquium: Strong-field phenomena in periodic systems,”Rev. Mod. Phys. , 021002 (2018).[9] J. Li, J. Lu, A. Chew, S. Han, J. Li, Y. Wu, H. Wang,S. Ghimire, and Z. Chang, “Attosecond science based onhigh harmonic generation from gases and solids,” NatureCommunications , 2748 (2020).[10] T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet,M. T. Hassan, and E. Goulielmakis, “Extreme ultravi-olet high-harmonic spectroscopy of solids,” Nature ,498 (2015).[11] M. Sivis, M. Taucer, G. Vampa, K. Johnston, A. Staudte,A. Y. Naumov, D. M. Villeneuve, C. Ropers, and P. B.Corkum, “Tailored semiconductors for high-harmonic op-toelectronics,” Science , 303–306 (2017).[12] S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, “High-harmonic generationby field enhanced femtosecond pulses in metal-sapphirenanostructure,” Nat. Commun. , 13105 (2016).[13] G. Vampa, B. G. Ghamsari, S. Siadat Mousavi, T. J.Hammond, A. Olivieri, E. Lisicka-Skrek, A. Y. Naumov,D. M. Villeneuve, A. Staudte, P. Berini, and P. B.Corkum, “Plasmon-enhanced high-harmonic generationfrom silicon,” Nat. Phys. , 659–662 (2017).[14] M. Garg, H. Y. Kim, and E. Goulielmakis, “Ultimatewaveform reproducibility of extreme-ultraviolet pulses byhigh-harmonic generation in quartz,” Nat. Photonics ,291–296 (2018).[15] Y. Yang, J. Lu, A. Manjavacas, T. S. Luk, H. Liu,K. Kelley, J.-P. Maria, E. L. Runnerstrom, M. B. Sin-clair, S. Ghimire, and I. Brener, “High-harmonic gen-eration from an epsilon-near-zero material,” Nat. Phys. , 1022–1026 (2019).[16] G. Vampa, T. J. Hammond, N. Thir´e, B. E. Schmidt,F. L´egar´e, C. R. McDonald, T. Brabec, D. D. Klug, andP. B. Corkum, “All-optical reconstruction of crystal bandstructure,” Phys. Rev. Lett. , 193603 (2015).[17] A. J. Uzan, G. Orenstein, ´A. Jim´enez-Gal´an, C. Mc-Donald, R. E. F. Silva, B. D. Bruner, N. D. Klimkin, V. Blanchet, T. Arusi-Parpar, M. Kr¨uger, A. N. Rubtsov,O. Smirnova, M. Ivanov, B. Yan, T. Brabec, and N. Du-dovich, “Attosecond spectral singularities in solid-statehigh-harmonic generation,” Nat. Photonics , 183–187(2020).[18] T. Huang, X. Zhu, L. Li, X. Liu, P. Lan, and P. Lu,“High-order-harmonic generation of a doped semiconduc-tor,” Phys. Rev. A , 043425 (2017).[19] S. Almalki, A. M. Parks, G. Bart, P. B. Corkum,T. Brabec, and C. R. McDonald, “High harmonic gen-eration tomography of impurities in solids: Conceptualanalysis,” Phys. Rev. B , 144307 (2018).[20] C. Yu, K. K. Hansen, and L. B. Madsen, “Enhancedhigh-order harmonic generation in donor-doped band-gapmaterials,” Phys. Rev. A , 013435 (2019).[21] K. Chinzei and T. N. Ikeda, “Disorder effects on the ori-gin of high-order harmonic generation in solids,” Phys.Rev. Research , 013033 (2020).[22] H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, andD. A. Reis, “High-harmonic generation from an atomi-cally thin semiconductor,” Nat. Phys. , 262 (2017).[23] T. T. Luu and H. J. W¨orner, “Measurement of the berrycurvature of solids using high-harmonic spectroscopy,”Nat. Commun. , 916 (2018).[24] D. Bauer and K. K. Hansen, “High-harmonic generationin solids with and without topological edge states,” Phys.Rev. Lett. , 177401 (2018).[25] R. E. F. Silva, ´A. Jim´enez-Gal´an, B. Amorim,O. Smirnova, and M. Ivanov, “Topological strong-fieldphysics on sub-laser-cycle timescale,” Nat. Photonics ,849 (2019).[26] A. Chac´on, D. Kim, W. Zhu, S. P. Kelly, A. Dauphin,E. Pisanty, A. S. Maxwell, A. Pic´on, M. F. Ciappina,D. E. Kim, C. Ticknor, A. Saxena, and M. Lewenstein,“Circular dichroism in higher-order harmonic generation:Heralding topological phases and transitions in chern in-sulators,” Phys. Rev. B , 134115 (2020).[27] C. J¨urß and D. Bauer, “Helicity flip of high-order har-monic photons in haldane nanoribbons,” Phys. Rev. A , 043105 (2020).[28] Y. Bai, F. Fei, S. Wang, N. Li, X. Li, F. Song,R. Li, Z. Xu, and P. Liu, “High-harmonic generationfrom topological surface states,” Nat. Phys. (2020),10.1038/s41567-020-01052-8.[29] D. Baykusheva, A. Chac´on, D. Kim, D. E. Kim, D. A.Reis, and S. Ghimire, “Strong-field physics in three-dimensional topological insulators,” Phys. Rev. A ,023101 (2021).[30] E. Runge and E. K. U. Gross, “Density-functional theoryfor time-dependent systems,” Phys. Rev. Lett. , 997–1000 (1984).[31] N. Tancogne-Dejean, O. D. M¨ucke, F. X. K¨artner, andA. Rubio, “Impact of the electronic band structure inhigh-harmonic generation spectra of solids,” Phys. Rev.Lett. , 087403 (2017).[32] C. Yu, H. Iravani, and L. B. Madsen, “Crystal-momentum-resolved contributions to multiple plateausof high-order harmonic generation from band-gap ma-terials,” Phys. Rev. A , 033105 (2020).[33] D. Golde, T. Meier, and S. W. Koch, “High harmon-ics generated in semiconductor nanostructures by thecoupled dynamics of optical inter- and intraband exci-tations,” Phys. Rev. B , 075330 (2008). [34] M. Kira and S. W. Koch, Semiconductor Quantum Optics (Cambridge University Press, 2012).[35] P. B. Corkum, “Plasma perspective on strong field mul-tiphoton ionization,” Phys. Rev. Lett. , 1994–1997(1993).[36] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier,and P. B. Corkum, “Theory of high-harmonic generationby low-frequency laser fields,” Phys. Rev. A , 2117(1994).[37] G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug,P. B. Corkum, and T. Brabec, “Theoretical analysisof high-harmonic generation in solids,” Phys. Rev. Lett. , 073901 (2014).[38] G. Vampa, C. R. McDonald, G. Orlando, P. B. Corkum,and T. Brabec, “Semiclassical analysis of high harmonicgeneration in bulk crystals,” Phys. Rev. B , 064302(2015).[39] A. M. Parks, G. Ernotte, A. Thorpe, C. R. McDonald,P. B. Corkum, M. Taucer, and T. Brabec, “Wannierquasi-classical approach to high harmonic generation insemiconductors,” Optica , 1764 (2020).[40] C. R. McDonald, G. Vampa, P. B. Corkum, andT. Brabec, “Interband bloch oscillation mechanism forhigh-harmonic generation in semiconductor crystals,”Phys. Rev. A , 033845 (2015).[41] X. Zhang, J. Li, Z. Zhou, S. Yue, H. Du, L. Fu, andH.-G. Luo, “Ellipticity dependence transition induced bydynamical bloch oscillations,” Phys. Rev. B , 014304(2019).[42] R. Liu and B. Zhu, “High-order thz-sideband genera-tion in semiconductors,” AIP Conf. Proc. , 1455–1456(2007).[43] B. Zaks, R. B. Liu, and M. S. Sherwin, “Experimen-tal observation of electron-hole recollisions,” Nature ,580 (2012).[44] F. Langer, M. Hohenleutner, C. P. Schmid, C. Poell-mann, P. Nagler, T. Korn, C. Sch¨uller, M. S. Sherwin,U. Huttner, J. T. Steiner, S. W. Koch, M. Kira, andR. Huber, “Lightwave-driven quasiparticle collisions ona subcycle timescale,” Nature , 225 (2016).[45] H. B. Banks, Q. Wu, D. C. Valovcin, S. Mack, A. C. Gos-sard, L. Pfeiffer, R.-B. Liu, and M. S. Sherwin, “Dynam-ical birefringence: Electron-hole recollisions as probes ofberry curvature,” Phys. Rev. X , 041042 (2017).[46] F. Langer, C. P. Schmid, S. Schlauderer, M. Gmitra,J. Fabian, P. Nagler, C. Sch¨uller, T. Korn, P. G. Hawkins,J. T. Steiner, U. Huttner, S. W. Koch, M. Kira, andR. Huber, “Lightwave valleytronics in a monolayer oftungsten diselenide,” Nature , 76 (2018).[47] J. A. Crosse, X. Xu, M. S. Sherwin, and R. B. Liu,“Theory of low-power ultra-broadband terahertz side-band generation in bi-layer graphene,” Nat. Commun. , 4854 (2014).[48] L. Yue and M. B. Gaarde, “Imperfect recollisions in high-harmonic generation in solids,” Phys. Rev. Lett. ,153204 (2020).[49] L. V. Keldysh, “Ionization in the field of a strong elec-tromagnetic wave,” Zh. Eksp. Teor. Fiz. , 1945 (1964),[Sov. Phys. JETP , 1307 (1965)].[50] L. Li, P. Lan, X. Zhu, T. Huang, Q. Zhang, M. Lein, andP. Lu, “Reciprocal-space-trajectory perspective on high-harmonic generation in solids,” Phys. Rev. Lett. ,193901 (2019). [51] P. Blaha, K. Schwarz, G. K. Madsen, D. Kvasnicka, andJ. Luitz, WIEN2K: An Augmented Plane Wave + Lo-cal Orbitals Program for Calculating Crystal Properties (Vienna University of Technology, Austria).[52] G. Kresse and J. Furthm¨uller, “Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set,” Phys. Rev. B , 11169–11186 (1996).[53] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek,C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira,S. W. Koch, and R. Huber, “Sub-cycle control of ter-ahertz high-harmonic generation by dynamical bloch os-cillations,” Nat. Photonics , 119 (2014).[54] I. Floss, C. Lemell, G. Wachter, V. Smejkal, S. A. Sato,X.-M. Tong, K. Yabana, and J. Burgd¨orfer, “Ab initiomultiscale simulation of high-order harmonic generationin solids,” Phys. Rev. A , 011401(R) (2018).[55] D. Xiao, M.-C. Chang, and Q. Niu, “Berry phase ef-fects on electronic properties,” Rev. Mod. Phys. , 1959(2010).[56] J. Li, X. Zhang, S. Fu, Y. Feng, B. Hu, and H. Du,“Phase invariance of the semiconductor bloch equations,”Phys. Rev. A , 043404 (2019).[57] F. Navarrete, M. F. Ciappina, and U. Thumm, “Crystal-momentum-resolved contributions to high-order har-monic generation in solids,” Phys. Rev. A , 033405(2019).[58] W. V. Houston, “Acceleration of electrons in a crystallattice,” Phys. Rev. , 184–186 (1940).[59] J. B. Krieger and G. J. Iafrate, “Time evolution of blochelectrons in a homogeneous electric field,” Phys. Rev. B , 5494 (1986).[60] N. Tancogne-Dejean, O. D. M¨ucke, F. X. K¨artner, andA. Rubio, “Ellipticity dependence of high-harmonic gen-eration in solids originating from coupled intraband andinterband dynamics,” Nat. Commun. , 745 (2017).[61] N. Yoshikawa, T. Tamaya, and K. Tanaka, “High-harmonic generation in graphene enhanced by ellipticallypolarized light excitation,” Science , 736–738 (2017).[62] ´O. Zurr´on, A. Pic´on, and L. Plaja, “Theory of high-orderharmonic generation for gapless graphene,” New J. Phys. , 053033 (2018).[63] S. Ghimire, A. D. DiChiara, E. Sistrunk, G. Nd-abashimiye, U. B. Szafruga, A. Mohammad, P. Agostini,L. F. DiMauro, and D. A. Reis, “Generation and prop-agation of high-order harmonics in crystals,” Phys. Rev.A , 043836 (2012).[64] P. G. Hawkins, M. Y. Ivanov, and V. S. Yakovlev, “Effectof multiple conduction bands on high-harmonic emissionfrom dielectrics,” Phys. Rev. A , 013405 (2015).[65] T.-Y. Du, D. Tang, X.-H. Huang, and X.-B. Bian, “Mul-tichannel high-order harmonic generation from solids,”Phys. Rev. A , 043413 (2018).[66] A. Taghizadeh, F. Hipolito, and T. G. Pedersen, “Linearand nonlinear optical response of crystals using lengthand velocity gauges: Effect of basis truncation,” Phys.Rev. B , 195413 (2017).[67] D. Vanderbilt, Berry phases in electronic structure the-ory: electric polarization, orbital magnetization and topo-logical insulators (Cambridge University Press, 2018).[68] L. Yue and M. B. Gaarde, “Structure gauges and lasergauges for the semiconductor bloch equations in high-order harmonic generation in solids,” Phys. Rev. A ,053411 (2020). [69] We mention that for the ERM calculations here, we haveneglected the terms in the saddle point equations (9) in-volving D k ,µ due to numerical complexities associatedwith their evaluation. However, due to the agreement be-tween the quantum and semiclassical calculations, as wellas the fact that these terms are small in the Γ − K case,we believe this is a good approximation. [70] E. Kane, “Zener tunneling in semiconductors,” J. Phys.Chem. Solids , 181 (1960).[71] F. I. Gauthey, B. M. Garraway, and P. L. Knight, “Highharmonic generation and periodic level crossings,” Phys.Rev. A , 3093 (1997).[72] M. Wu, D. A. Browne, K. J. Schafer, and M. B. Gaarde,“Multilevel perspective on high-order harmonic genera-tion in solids,” Phys. Rev. A94