Theoretical description of high-order harmonic generation in solids

We consider several aspects of high-order harmonic generation in solids: the effects of elastic and inelastic scattering; varying pulse characteristics; and inclusion of material-specific parameters through a realistic band structure. We reproduce many observed characteristics of high harmonic generation experiments in solids including the formation of only odd harmonics in inversion-symmetric materials, and the nonlinear formation of high harmonics with increasing field. We find that the harmonic spectra are fairly robust against elastic and inelastic scattering. Furthermore, we find that the pulse characteristics play an important role in determining the harmonic spectra.


I. INTRODUCTION
The rapid development of time-resolved pump-probe experiments such as angle-resolved photoemission spectroscopy (ARPES) 1-3 , optical reflectivity 4 and X-ray scattering 5,6 has expanded greatly the use of strong ultrashort laser pulses in the study of condensed matter systems. The advent of these new experimental techniques requires theory to properly describe the non-equilibrium physics observed in these measurements. The strong temporal dependence of these problems makes analysis in the frequency domain difficult which is more common for problems in equilibrium or within linear-response. A variety of methods have been developed or adapted to describe the behavior of these systems directly in the time domain including semiclassical approximations, 7-10 Floquet theory 11 and direct numerical integration of the time-dependent Schrödinger equation. 12 Each of these methods has limitations varying from lack of numerical feasibility to the inability to explicitly treat time-translation invariance breaking due to the pump and/or probe fields.
The non-equilibrium Keldysh formalism, developed by Kadanoff and Baym 13 and Keldysh, 14 can treat broken time-translation invariance explicitly making use of the so-called Keldysh contour (see Fig. 1). Quantities such as the single-particle Green's function are defined along this contour with multiple arguments that take appropriate real and imaginary times along the contour. Characteristics of the pump pulse, such as the driving frequency or pulse length and shape, as well as material specific tight-binding or more realistic band structure parameters for the equilibrium system can be included directly allowing for a closer match between theoretical calculations and realistic experimental conditions. The ab initio inclusion of nonlinear effects of the driving field, rather than treating it as a perturbation, is necessary for a proper description of the electronic system as it is driven out of equilibrium. This Keldysh formalism has been adapted and used to study the effect of electronic correlations on systems out of equilibrium within dynamical mean-field theory, 15-22 , continous-time Quantum Monte Carlo 21,23 and cluster perturbation theory 24 . A number of other methods exist to treat correlated and uncorrelated time-domain phenomena, including non-equilibrium exact diagonalization 25 , numerical renormalization group 26 and density matrix renormalization group, 27,28 with various benefits and limitations.
An ideal application for this Keldysh approach is in the area of high-order harmonic generation (HHG) where one need only consider the non-linear effects of the pump, and can neglect details associated with probe pulses and the often complicated evaluation of a proper pump-probe measurement cross-sections. High-order harmonics result from nonlinear light-matter interactions and in the present example from the manner in which the sample current varies in time due to the pulse characteristics of the driving field. Used in atomic gases for the production of attosecond pulses, 29 HHG opens avenues for time-domain experiments with high temporal resolution, and in particular are short enough to resolve electronic interactions which occur on these fundamental timescales. A second use for highorder harmonics is in the production of a wide spectrum of frequencies not available in the fundamental frequencies of standard lasers. By applying strong driving fields, the harmonics generated can be harvested for further use as a pump or probe.
Trains of attosecond pulses that can then be used for further experiments 30,31 can be generated through the interaction of a femtosecond laser pulse with a jet of gas. In this process the pulses of light are created in three-steps: (1) field-assisted tunnel ionization of an atom, (2) acceleration of the ionized electron back to the atom by the change in sign of the field as the driving pulse oscillates and (3) creation of a train of high-energy photons as the electron re-collides with the charged ion. 32,33 In addition to this generation in gaseous media, HHG can be accomplished through non-linear response in solids. The availability of ultra-short laser pulses with relatively strong fields has facilitated studies in solid state materials while mitigating the risk of sample damage. This capability has led to recent observations of HHG in ZnO crystals. 34 Their bulk crystalline nature makes HHG in these materials fundamentally different from what occurs in atoms or small molecules.
The electron ionization in an atom is replaced by electron-hole pair excitation between the conduction and valence bands, respectively. The charged particles are then accelerated in the field with a non-linear response due to the Bloch oscillations which occur when the field strength is large enough, as discussed in some detail by several authors. 9,10 This is fundamentally different from the single-atom picture, where the harmonics are generated by the stimulated ionization and recombination of an electron in the outer shell of an atom. In a solid, this would correspond to electron-hole pair excitation and de-excitation.
However, as we will show, high-order harmonics can be generated solely via the acceleration of the conduction electrons by the field, and the Bragg scattering that occurs repeatedly in high field. 34 The end result is HHG of the fundamental frequency of the applied pump pulse. Although this has been previously studied for non-interacting systems, 9,10,34 using the non-equilibrium Keldysh formalism, we can go beyond the semi-classical approach and consider the effects of elastic and inelastic scattering on the HHG spectra.
A better understanding of HHG in solids requires a proper description of the underlying electronic states in the material. In atomic gases, these states have been modeled in a number of ways, from partial-wave expansion of the electron cloud to other descriptions utilizing first-principles methods. 32 In solids, Bloch states are a basis in which one can describe the itinerant conduction band electrons and valence band holes. In the presence of a strong driving electric field, the states are modified Bloch states that, as the field energy becomes comparable to the bandwidth, become localized. 35 . The behavior covering the full range of fields, is captured by the inclusion of the field via the Peierls' substitution.
The paper proceeds as follows. In Section II, we present relevant details of the nonequilibrium Keldysh formalism as it pertains to HHG. Section III discusses the effects of impurity and phonon scattering. In section IV we consider the effect of different pulse characteristics on HHG and compare and contrast the resulting HHG spectra for two experimentally relevant cases. Section V focuses on the time structure of the individual harmonics.
Finally, section VI studies the HHG spectra from a more realistic, complex band structure appropriate for Si.

II. METHOD
We shall consider a system of electrons in the time-domain using the non-equilibrium Keldysh formalism. The appropriate expressions were derived in detail 15,36 previously and we will present the relevant equations here for the specific case of a spatially uniform field.
The system starts in an initial equilibrium state at t = t min , at an initial temperature T , and evolves until a final time t = t max . We describe the system on the Keldysh contour shown in Fig. 1  is done in steps of δ t and iδ τ on the real and imaginary parts of the contour, respectively. This is different from the previous implementation by Turkowski and Freericks, who evaluated the temporal integrals for the Dyson equation 37 directly. 15 However, the direct evaluation of integrals is limited to the case of simple band structures, or, the infinite-dimensional case considered in the dynamical mean-field theory treatment of the Falicov-Kimball model. 38 Here, we numerically evaluate the integrals on a fine grid.
The bare non-equilibrium Green's function G 0 under the influence of a driving field can be obtained through the Peierls' substitution 39 in the standard definition of the double-time , T is the temperature (we work in units where the Boltzmann constant k B = 1), t and t lie on the Keldysh contour, θ is the contour-ordered version of the standard Heaviside function, ( k) is the single-particle dispersion, A(t) is the vector potential, which is related to the field by where we use e = c =h = 1 and work in the Hamiltonian gauge. The integral in the exponent runs along the Keldysh contour.
To illustrate the influence of the electronic self-energy on HHG in solids, we include the effect of impurities and electron-phonon scattering. For impurities, we consider dilute point-like scatterers, for which the self-energy is 40 where n i and V i denote the impurity concentration and potential, respectively. For electronphonon scattering in the Migdal limit of the Holstein model 41 , we consider a single optical phonon mode with frequency Ω coupled with a constant g, for which the self-energy is 37 where n B (Ω) = 1/(exp(Ω/T ) − 1). The full Green's function is then self-consistently calculated through Dyson's equation expressed on the Keldysh contour as After discretization, this equation can be cast in the form of a matrix equation where all times including those of the differentials ∆ t 1,2 lie on a branch of the Keldysh contour. This matrix equation can be solved using standard linear algebra techniques.
The process by which high-order harmonics are generated in undoped semiconducting crystals begins with tunnel ionization of electrons from the valence band into the conduction band by the pump laser. Next, this same pulse accelerates the electrons, creating an oscillating current. Here, the details of the tunnel ionization process are neglected and instead we focus on the motion of electrons only after they have been promoted to the conduction band. In principle, there is a time dependence due to the tunneling process since the tunneling rate depends on the field, which is time dependent. However, since the tunneling rate depends exponentially on the field, the majority of electrons injected into the conduction band will tunnel there when the field is largest, and will be essentially instantaneous. 34 The subsequent driving force, in the absence of scattering will give rise to a similar harmonic spectrum as if the electrons had been there to begin with, albeit with somewhat broader harmonics and some higher frequency content due to the nearly instantaneous injection of the electrons into the valence band. Finally, as pointed out by Golde et al., 9 the interaction between the excited electrons and the valence band holes leads to a contribution to the harmonic spectrum, in particular at very high harmonic orders. Here, we shall focus on the intra-band contribution, which is several orders of magnitude stronger.
The instantaneous current j(t) created by the applied field is related to the lesser Green's function by where v( k) = ∂ ( k)/∂ k is the group velocity of a wave packet. We assume that the accelerated electrons radiate proportional to the Fourier transform of the instantaneous current, with the HHG spectrum proportional to ˆd t e iωt d dt We take the Fourier transform of the instantaneous current after interpolating with a shapepreserving Akima spline 42 . The self-energy and current are calculated on a 100 × 100 grid of momentum points in the two-dimensional square Brillouin zone, and 24 × 24 × 24 grid in the three-dimensional Brillouin zone appropriate for the Si band structure, sufficient to give converged results for the current.
Since we are considering a laser pulse, a DC contribution to the electric field is excluded, as required by Maxwell's equations; any DC component does not propagate. Note that this restriction applies to the specific case of a propagating pulse; a DC contribution would be allowed for statically applied fields. A corollary to the absence of a DC component of the field is that the current is also zero at long times. This fact, combined with the relatively simple form of the current for a non-interacting system allows us to limit the length of the Keldysh contour to just the length of the pulse. Our discretization on the real-time branches of the Keldysh contour varies depending on the section of the paper, but is always kept above or equal to δ t = 0.015t −1 h , allowing up to the 30 th harmonic to be visible (according to the Nyquist theorem) 43 . We always use a discretization of 10 −2 T −1 along the imaginary-time spur, and set the simulation temperature to 0.1 eV. All times are measured in units of t −1 h , and applied fields in units of t h . This can be converted to real units using appropriate values for the material under study.
For sample values of the hopping t h = 1 eV and material lattice constant or characteristic length scale a 0 = 5Å, a field strength of E = t h corresponds to approximately 3.2 MV/cm (32 mV/Å); the temporal unit in this work, t −1 h would correspond to approximately 4.14 fs. Note that these values are not uniformly applicable; t h and a 0 need to be adjusted for the material under study.

III. SCATTERING BY IMPURITIES AND PHONONS
Motivated by the transition metal oxides, we consider a square lattice with a model band structure (in the absence of an electric field) given by where t h = 0.3t h and µ = −0.255t h . We treat the system as isolated -unconnected to a bath -so the chemical potential remains at the initial equilibrium value. The filling mainly affects the amplitude of the current, and does not affect qualitatively the results . The system is driven by a Gaussian-shaped pump pulse centered at t = 0 with a duration of 4.24t −1 h , and a single underlying fundamental frequency ω a = (9/5)π which corresponds to roughly 9 full cycles under the Gaussian envelope. Since it is a pulse with finite width, the Fourier transform will contain a range of frequencies centered around ω a , as will be discussed below.

A. Harmonic generation from non-interacting electrons
We calculate the current as a function of time for an oscillating field with a maximum field amplitude E max , as shown in Fig. 2. At the beginning of the pulse, the applied field amplitude is small and the field at that amplitude does not act very long; thus, the generated current is proportional to the field (up until t ∼ −2t  through impurity scattering, although oscillations have been observed in ultrasmall Josephson junctions, 44 semiconductor superlattices, 45,46 and optically trapped atoms. 47 However, the HHG spectrum is a direct consequence of the Bloch oscillations and provides an indirect measurement of them in solid-state systems, as discussed by various authors. 9,10,34 As the field strength grows, the effects from Bloch oscillations intensify, as seen near the center of the pump pulse which displays several oscillations of the current in a region where the field does not change direction. A time-varying current due to the applied field leads to radiation, we thus perform a Fourier transform of the current to obtain the HHG spectrum (shown in Fig. 2(c)). 48 On a semi-logarithmic scale, we clearly observe peaks near odd multiples of the central frequency of the applied field (ω a ). The fact that odd multiples dominate can be obtained from simple symmetry arguments and agrees with experimental observations. For inversion-symmetric systems, the velocity is anti-symmetric under parity, which leads to a cancellation of even harmonics of the current. 11 The peaks have a finite width resulting from the finite pulse length where as a continuous pulse would lead to sharp peaks at exactly odd harmonic frequencies.   This ratio increases as a function of field with the largest fields showing a ratio approaching one. For the largest fields, the 1 st is not necessarily the largest, which can be understood from the time structure of the individual harmonics discussed below.

B. Elastic scattering
To study simple effects from elastic scatterers on HHG, we consider a pump pulse whose field amplitude (A max =20t h ) produces a sufficient number of odd harmonics so that variations in the HHG spectrum as a function of impurity scattering strength can be tracked reasonably well. Figure 4 shows the HHG spectra for increasing impurity strength, where we define a In addition to impurity scattering, inelastic scattering can play a role in determining the current, and thus the HHG spectra. In particular, in semiconductors scattering due lattice vibrations (i.e. phonons) can play an important role. 49 Here, we consider the effect of a single phonon mode on the generated spectra (as in the Holstein model). 41 Due to computational restrictions, the pulse strength is limited to A max =5t h . We use a phonon frequency of Ω = 0.2t h , and define the electron-phonon coupling strength as λ = 2g 2 /ΩW , where W is the bandwidth and g is the electron-phonon coupling constant. With these parameters, we find the current and HHG spectra shown in Fig. 5. Similarly to the results for elastic scattering, although the overall current is strongly renormalized, the scaled HHG spectra only differ qualitatively at the highest frequencies, where the individual harmonics can no longer be resolved.  The HHG spectra can be affected by the particulars of the applied pulse. To illustrate this, we consider two types of pulses: the multi-cycle pulse discussed above, and a singlecycle pulse similar to those used in THz experiments. 50 For both pulses, although there is no DC component to the electric field, there is still a choice of the carrier-envelope phase.

IV. EFFECTS OF PULSE CHARACTERISTICS
The results here have an arbitrarily chosen phase; we present a detailed dependence of the results of said phase in Section V.
The field profiles for the two pulses considered are shown in Fig. 6. The inset shows the Fourier transform of the vector potential obtained from the field profile. Due to the relatively large number of oscillations in the multi-cycle pulse, the Fourier transform is strongly peaked at the driving frequency ω a ≈ 9π/5. On the other hand, the single-cycle pulse is temporally short, leading to a large number of Fourier components with the largest contribution at Figure 7 shows the HHG spectrum as a function of field amplitude and frequency. There is a clear distinction between the spectra associated with the multi-cycle and single-cycle pulses. As discussed above, the current induced by the multi-cycle pulse responds in odd harmonics of the fundamental frequency which is relatively well defined as seen in the inset of Fig. 6. However, for the single-cycle pulse, the response is more complex and we find two distinct regimes. When the field amplitude is small, the electronic system responds perturbatively as in the small amplitude regime for the multi-cycle pulse. As such, the dominant contribution to the HHG spectrum comes from that frequency which gives the largest contribution to the field, ω 0 . However, this linear response is limited to the fundamental, at a field amplitude where no harmonics can yet be observed. As the field grows larger, the electrons undergo Bragg reflection as discussed above, which leads to non-linear response in the current . Since the Fourier spectrum of the single-cycle pulse is very broad, the peaks in the Fourier-transformed current do not simply occur at odd harmonics of the fundamental applied frequency, but rather at frequencies that depend on the field strength and that are not linearly spaced. Additionally, the peaks observed are much broader than those of the multi-cycle pulse, which reflects the broad Fourier profile of the single-cycle pulse. Lastly, we note that the single-cycle pulse has a clear asymmetry about the pulse center. We have repeated the calculations with a symmetrized pulse, and find no qualitative difference in the HHG spectra.
In addition to the field frequency and amplitude, it is possible to experimentally control the carrier-envelope phase. This is directly addressable within the formalism presented, with the small caveat that there cannot be a DC component to the electric field. Thus, we choose to vary the phase φ in the vector potential: (9) Figure 8 shows the variation in the HHG spectra due to the carrier-envelope phase for both the 9-cycle and single-cycle pulse. The development of clear harmonic orders hinges on the applied field being sufficiently wide to clearly define the carrier frequency. Thus, the 9-cycle pulse does not show significant variation with the phase. However, the Fourier transform of the single-cycle pulse is very wide; the pulse is not long enough for a clear development of the carrier frequency. This causes a strong dependence on the carrier-envelope phase.

V. TIME STRUCTURE OF THE INDIVIDUAL HARMONICS
The HHG process in solids is due to the repeated Bragg scattering as the accelerated electrons hit the zone edge.This leads to a simple time structure in the generated harmonics (see also Refs. 10,34,48). Figure 9 shows the windowed Fourier transform (also known as a Gabor transform 51 ) for the 9-cycle pulse. The low order harmonics are mainly generated in the low-field regions of the pulse. As the pulse grows stronger, the low order harmonics are obscured by the higher orders, as the Bragg scattering occurs more rapidly. The higher orders are thus temporally localized to the center of the pulse. This can be clearly seen from both the time-dependent power spectra as well as the individual harmonics' contribution to the current.
Through temporal localization one can understand the increase of harmonic amplitudes above the 1 st . The low order harmonics are generated near the beginning and end of the pulse, and their duration is inversely related to the maximum field strength. Thus, for large field strengths, the time during which the 1 st order harmonic has a strong contribution is smaller than the next orders.   Fig. 3). As in the case of the simple 2D band structure, the observed intensity increases with the field amplitude. Furthermore, since the bandwidth of the tight-binding model is different from the 2D model band structure, the field strength at which the nonlinearity onsets is similarly different. According to Ghimire et al. 34 , inclusion of more complex terms beyond first order in the band structure leads to an enhancement or reduction of the harmonic content in the HHG spectra. Here we consider the multi-cycle pulse discussed above, applied along the (111) and (100) directions of bulk Si. The band structure in both directions is sufficiently complex that, after averaging over the full Brillouin zone, there is no singular qualitative difference which can be attributed to a particular feature of the band structure. There is an overall difference in the structure of the low field harmonic spectra, which appear to decay faster in the (100) direction compared to the (111) direction.

VI. SILICON BAND STRUCTURE
The lack of a strong distinction in the spectra agrees with more detailed calculations done based on density functional theory 53 .

VII. SUMMARY
In summary, we have considered several aspects of harmonic generation in solids. Harmonic generation from solids shows promise as a future tool as a tunable light source for further experiments. We find that the harmonic spectra are sensitive to the particular details of the driving pulse, which indicates a need to carefully control the parameters of said pulse.
On the other hand, the harmonic spectra are fairly robust against the influence of scattering, with the main feature being a decrease in the overall amplitude of the current. Finally, we find that for a realistic band structure, the HHG spectra do not depend sensitively on direction due to the complexity of the bands in any direction; however, this may not be the case in materials with simpler band structures. Overall, the most important message from this work is the very robust nature of the HHG response in solids, which does not seem to depend very strongly on the band structure, the orientation of the field, the scattering, or the pulse characteristics, except in the special cases detailed above.