Wannier quasi-classical approach to high harmonic generation in semiconductors

We develop a quasi-classical theory of high harmonic generation in semiconductors based on an interband current that has been transformed from Bloch to Wannier basis. The Wannier quasi-classical approach reveals a complete picture of the mechanisms shaping high harmonic generation, such that quantitative agreement with full quantum calculations is obtained. The intuitive picture revealed by quasi-classical wavepacket propagation will be helpful in the interpretation and design of high harmonic and attosecond experiments. Beyond that, the capacity to quantitatively model quantum dynamics with classical trajectories should prove useful for a wider spectrum of condensed matter research, including coherent control, transport theory, and strong field physics.

We develop a quasi-classical theory of high harmonic generation in semiconductors based on an interband current that has been transformed from Bloch to Wannier basis. The Wannier quasiclassical approach reveals a complete picture of the mechanisms shaping high harmonic generation, such that quantitative agreement with full quantum calculations is obtained. The intuitive picture revealed by quasi-classical wavepacket propagation will be helpful in the interpretation and design of high harmonic and attosecond experiments. Beyond that, the capacity to quantitatively model quantum dynamics with classical trajectories should prove useful for a wider spectrum of condensed matter research, including coherent control, transport theory, and strong field physics.
Although some experimental features can be reasonably well reproduced by numerical models [14][15][16][17], a thorough understanding of all the components shaping harmonic spectra is still missing. This inhibits progress in optimizing HHG as a radiation source and in further developing HHG as a diagnostic tool.
The principal mechanism of interband HHG has been clarified by saddle point integration of the interband current derived in the Bloch basis [18]. Electron and hole are born at the same lattice site in real space by tunnel ionization and quiver in the laser field. When they recollide at some lattice site, a harmonic photon is emitted. Its energy is equal to the bandgap at the crystal momentum of the electron-hole pair at recollision. Despite its merits, the Bloch quasi-classical model falls short of accounting for the lattice structure; quantum mechanics allows recombination of electrons and holes at different lattice sites, as was clearly demonstrated in recent work [10,19,20].
Here we develop a generalized quasi-classical approach that accounts for the lattice structure; this is achieved by transforming the interband current from Bloch to Wannier basis followed by saddle point integration. The basis change has a substantial effect. The resulting Wannier quasi-classical (WQC) model is found to be in quantitative agreement with quantum calculations. So far, quasi-classical k-space analysis has been used to qualitatively investigate strong field effects in gases and in the condensed matter phase; quantitative agreement has not been demonstrated yet. Whether quantitative agreement can be obtained in the Bloch basis remains to be seen, however the richer physics revealed by the WQC picture indicates that this might not be the case. The more refined WQC picture arises from the fact that the transition dipole moment enters the classical action in the exponent, and therewith the saddle point equations.
The quantitative agreement with full quantum calculations suggests that the physical picture for HHG in semiconductors revealed by the WQC analysis is complete. An electron and hole can ionize and recombine at different lattice sites with a probability determined by the tunneling exponent and Wannier dipole moments; birth and recombination sites are connected by classical trajectories; quantum effects are included by a quadratic expansion of the classical action about the classical trajectories. Beyond that, our WQC analysis allows unprecedented insight into the real-space aspects of tunnel ionization in solids; it gives access to the tunnel ionized wavefunction in real space and therewith, to the birth location of the electron hole pair.
More generally, our analysis opens an avenue for modeling quantum dynamics of wavepackets by propagating classical trajectories. This is potentially relevant for a wide spectrum of applications ranging from strong field physics to transport phenomena [21,22] and coherent control [23,24]. On a fundamental level, the WQC approach could open an alternative pathway to modeling noise and few electron-hole dynamics in solids; as propagation from initial to final Wannier wavepacket is done by classical trajectories, the space in between does not need to be resolved in contrast to a full quantum approach.

A. Two Band WQC Model
Our formalism is developed for a 3D, two-band model. We first summarize derivation of HHG in the Bloch ba-sis [18]; it starts from the time-dependent Hamiltonian H(t) = H 0 + x · F(t); F(t) represents the laser field; H 0 is the unperturbed lattice Hamiltonian with Bloch eigenstates Φ m,k (x) = 1/ √ V u m,k (x) exp(ik · x) and with energies E m (k) in band m with crystal momentum k; the band index m = v, c refers to valence and conduction band, respectively; u m,k is the periodic part of the Bloch function, Φ m,k |Φ m,k = 1, and u m,k |u m,k = υ. Finally, V = N υ is the volume of the solid, with N and υ the number and volume of primitive unit cells. Hartree atomic units are used, unless otherwise noted.
In the presence of the laser field the wavefunction becomes time-dependent. In the length gauge it is represented as where a m (k, t) are the probability amplitudes and integration is over the full Brillouin zone (BZ). As initial conditions we choose an empty conduction band a c (k, t = 0) = 0, and a filled valence band, a v (k, t = 0) = 1/ √ V BZ , where V BZ is the Brillouin zone volume. The Ansatz (1) is substituted into the time-dependent Schrödinger equation, and the interband polarization and current are found to be [18] , with d mm ′ (k) = i u m,k |∇ k |u m,k the transition dipole moment. For a two-band system, we denote and we assume a centro-symmetric system for which the diagonal elements d mm (k) can be set to zero [26]. In the following we will translate HHG, as described by the interband current of (2), from k-space to real space by using Wannier functions. The Bloch and Wannier basis functions are connected by a Fourier transform according to [27] Here, w m (x − x j ) is the Wannier function of band m corresponding to the primitive unit cell at position x j .
By virtue of (4b), the initial wavefunction, corresponds to the Wannier function at position x j = 0. HHG can start from any other site x j . The initial Wannier function can be shifted to x j by setting a v (k, t = 0) = exp(−ik · x j ). As all lattice sites are identical, it is sufficient to investigate x j = 0. In order to translate the interband current (2) into real space, the Bloch functions in the transition dipole moment (3) are replaced by the Wannier functions with the help of relation (4a). This leads to where the second line was obtained by setting x k = x j + x l and by replacing summation index k with l in the first line. Also, note that performing j in the second line changes the integration volume from a unit cell to the whole crystal volume. The Wannier dipole moments are equivalent to the Fourier series expansion coefficients of the Bloch dipole moment d(k). Interpreted in real space, the Wannier dipole moment d l describes a transition where an electron is born l lattice cells away from the hole. Bloch and Wannier dipole moments are not unique; is also an eigenfunction for any real function α that is periodic in k-space. Although the full equations, including the diagonal dipole elements d mm , are gauge invariant [25,26], it is computationally advantageous to choose strongly confined Wannier basis functions [28,29] in order to keep the number of relevant lattice sites small. In the 1D examples discussed further down we chose maximally localized Wannier basis functions [28] for which d mm = 0. Inserting (6) into (2), the interband current follows as Here x l ; P jl (ω) represents the probability amplitude that the harmonic ω is generated by an electron-hole pair that is born with a relative distance |x l | between electron and hole and later recombines with relative distance |x j |, and the propagator T jl describes the evolution between d * l and d j .

B. Saddle Point Integration
The integrals in (7b) are solved by saddle point integration. The saddle point equations, result from ∂ϕ/∂µ = 0 with µ = t ′ , t, k, respectively. The field quiver motion between times t ′ and t is given Note that the classical action depends on the difference between conduction and valence band. As a result, the above quantities represent the difference between electron and hole band velocity and excursion distance. Finally, the ∓ in (8b) accounts for the complex conjugate term in (2b).
The set of equations (8) are solved for a linearly polarized laser field F = Fx; further A = Ax and k x = k. The solutions of the saddle point equations are denoted by where we have approximated the bandgap as with i, j = x, y, z; β ij (k) = ∂ 2 ε/∂k i ∂k j the inverse mass tensor; and E g the minimum bandgap. The positive sign in (9) is chosen to obtain an exponentially decaying tunneling rate. The two remaining saddle point equations (8b) and (8c) determine t b and t r . They have to be solved numerically for each possible birth site x l and recombination site x j ; for instance, by running through t b and finding all t r (t b )'s that fulfill (8c). From those, the pairs [t b , t r ](ω) are selected that produce a given harmonic ω via (8b). The physical implications of the saddle point equations are discussed at the end of this subsection.
Next, the integrand of (7b) is evaluated at the saddle point, where the small imaginary birth time determines the tunneling exponent. Further, the phase ϕ is expanded to second order, which gives the multivariate Gaussian integral where q = (t ′ , t, k), and H is the Hessian H ij = ∂ 2 ϕ/∂ i ∂ j with i, j ∈ q. The full expression for the determinant of the Hessian is provided in appendix A. Putting everything together, we obtain the WQC propagator where g = ωF(t b + iδ)(2π) 5/2 / |H| and to leading order the determinant from the Gaussian integral |H| ≈ v x (k s )f (t b + iδ, t, k s ) [13], see appendix A. Further, it is convenient to split the phase in (12c) into dτ + ωt r contains the classical action and the harmonic frequency Fourier term. The second term is the Fourier term of the recombination dipole moment, χ 2 = k s · x j . The total probability amplitude is governed by the prefactor g, the ionization amplitude d * l e −tx , the quantum mechanical phase factor e −iχ1 acquired along the classical trajectory, and the recombination amplitude d j e −iχ2 . For each possible birth site x l and recombination site x j in the lattice, the summation runs over all birth and recombination times t r , t b that satisfy the saddle point conditions for a particular harmonic frequency ω.
The propagator (12) together with the saddle point equations (8) and the interband current (7a) represent the WQC description of HHG in semiconductors. They reveal a complete and detailed picture of the physical mechanisms driving HHG in real and reciprocal space, summarized in figures 1(a) and (b), respectively. The empty circles in figure 1(a) represent the centers of the atomic unit cells x l , where l = (l x , l y ) in the 2D schematic. A Wannier basis function is located at each center. Initially, all Wannier sites of the valence band are filled. As all lattice sites are identical, it is sufficient to investigate x l = 0, see below (5). Following the notation of our calculation we chose indices l, j to represent birth and recombination sites, respectively. HHG proceeds in three steps.
Step 1 -creation of electron-hole pair by ionization. At birth time t b , a valence band electron localized at lattice site x 0 transitions to the conduction band, and is localized at lattice site x 0 + x l . The tunneling probability is determined by the tunneling exponent t x and by the Wannier dipole moment d * l , see figure 1(a). The potential energy experienced by the created electron-hole dipole in the laser field makes the effective ionization potential E g + F (t b )x l birth site dependent, see (8a) and (12b). In reciprocal space the electron transitions from valence to conduction band at the Γ-point at time t b , see figure 1(b).
Step 1 is of quantum mechanical nature.
Step 2 -electron-hole evolution in laser field. The electron-hole pair quivers in the laser field. In real space it follows the classical trajectory ξ(t b , t r ) in figure 1(a) until electron and hole revisit each other and are separated by |x j | at time t r , see (8c). The propagation step is dominantly classical; of quantum mechanical nature are the phase χ 1 (t b , t r ) picked up between birth and recombination time, and the quasiclassical factor g coming from the quadratic expansion of the classical action S about the classical trajectory. The shaded green area about the classical trajectory in figure 1(a) indicates the quantum correction up to second order. In reciprocal space in figure 1(b) the electron-hole pair evolves from initial crystal momentum zero to saddle point crystal momentum k s (t b , t r ), defined below (8).
At time t r electron and hole recombine with probability amplitude d j e −iks(t b ,tr)·xj , see figure 1(a). The harmonic energy is given by the bandgap energy at k s (t b , t r ), see figure 1(b), plus the energy of the electron hole dipole in the field F (t r ), see (8b). Due to the second term, harmonics with energies somewhat larger than the maximum bandgap can be generated.

III. RESULTS
For the remainder of the paper, the WQC approach and its physical significance are explored within a 1D model system. In this case the interband current, WQC propagator, and probability amplitude reduce to scalars; namelyj er , T jl , and P jl . Specifically, we use a 1D delta function model potential, V (x) = Ω ∞ n=−∞ δ[x − (n + 1/2)a] with unit cell size a and barrier penetration parameter Ω. Details of the delta function model are given in appendix B. For the investigated parameters the bandgap is well approximated by the nearest neighbour dispersion ε(k) = E g + ∆[1 − cos(ka)], where E g is the minimum bandgap and 2∆ represents the bandwidth. We chose a = 7 and considered two values Ω = 0.5, 1.5 to model a weakly and tightly bound semiconductor, respectively. The corresponding bandgap parameters are E g = 0.141, 0.269; ∆ = 0.269, 0.17. Finally, for all runs we use a dephasing time T 2 = T 0 /2 so that only returns within a single cycle are relevant.
For the WQC calculation we assume the continuous wave limit, F (t) = F 0 sin(ω 0 t), in order to facilitate interpretation of the results. Equation (12a) has been derived for finite pulses employing the Fourier transform. For a transition to the cw limit, the Fourier transform has to be replaced by a Fourier series; as a result, ω → nω 0 , pre-factor g → g/(2πT 0 ), where the 1/(2π) comes from the 1D nature of our model. The harmonic yield becomes |h n | 2 = |j er (nω 0 )| 2 with T jl given by the WQC propagator (12a).
In figure 2 the blue empty circles (exact) and blue filled circles (WQC) refer to results for the weakly bound model semiconductor, with Ω = 0.5, ω 0 = 0.01425. Red empty squares (exact) and red filled squares (WQC) refer to the tightly bound semiconductor, with Ω = 1.5, ω 0 = 0.0285. Plots with the same symbols in figures 2(a) and (b) correspond to the same values of Ω and ω 0 , but differ in F 0 .
The WQC approach agrees well with the exact solution, with most data points being off by less than a factor 2. Even the first 1-2 cutoff harmonics are described fairly well, which demonstrates that they are of quasiclassical origin. The good agreement allows us to interpret semiconductor quantum dynamics such as ioniza- tion, electron/hole transport, and HHG in terms of clas- sical trajectories. The quantum contributions to HHG are captured by the tunneling exponent t x , by the preexponential factor g in (12a), and by the Wannier dipole moments in (7).
A few points disagree by a larger factor of up to 6. In particular, figure 2(a) shows that the WQC result for harmonic n = 15 exhibits larger discrepancy for the weakly bound semiconductor (Ω = 0.5) compared to the more tightly bound semiconductor (Ω = 1.5). The reason for this behaviour is identified in figure 3 and will be discussed later.
Numerical solution of the full saddle point equations reveals two distinct classical trajectories that contribute to the probability amplitude P jl ; one long trajectory and one short. Moreover, each solution exists for only certain combinations of birth (l) and recombination (j) lattice sites. Figure 3 shows the contributions arising from the different classical trajectories for the fifteenth harmonic (n = 15) with Ω = 0.5, F 0 = 0.0025, corresponding to the filled blue circles in figure 2(a). Figure 3(a) depicts the regions in the j-l plane where each trajectory contributes to |P jl (n = 15)|. No solution exists for the dark region in the top-right, and the probability amplitude here is zero. Figures 3(b) and (c) show the individual contributions to the probability amplitude from the long and short trajectories, respectively. The long trajectory is dominant, as the electron-hole pair is born close to the field peak, whereas the short trajectory is born closer to the nodal point. This outweighs the effect of the short dephasing time, which favors the short trajectory. As a result, the contribution of each data point to the WQC propagator is dominantly determined by a factor ∼ ge −tx of a single (long) trajectory. The full probability amplitude |P jl (n = 15)| is essentially identical to figure 3(b).
In figure 4 the total probability amplitude for the fifteenth harmonic |P jl (n = 15)| is plotted as a function of birth and recombination site indices l, j for Ω = 1.5, F 0 = 0.008, which corresponds to the filled red squares in figure 2(a). For this system the long trajectory is also dominant, and analysis of the individual contributions would reveal a picture qualitatively similar to figure 3.
In both figures 3 and 4, harmonic n = 15 has been selected, as the WQC result for the weakly bound semiconductor exhibits a more pronounced difference, while it agrees well for the tightly bound semiconductor. For both systems, the maximum probability is shifted towards negative birth site indices; it is more likely for electron and hole to be born apart than at the same site. Tunnel ionization probability is determined by e −tx and by birth dipole moment d * l . The tunnel exponent t x depends on the ionization potential E g + F (t b )x l , see (12a). Thus, for positive field the electron-hole pair gains energy when born at increasingly negative distances which reduces t x . When −x l = E g /F (t b ), t x vanishes; in other words, the valence and conduction band levels separated by −l sites align, and the electron hops from the valence to the conduction band site. The penalty to be paid is a rapidly dropping dipole moment d l . As such, the birth site index at which ionization is maximum is determined by a tradeoff between tunnel exponent and Wannier dipole moment. The dipole elements for the parameters of figure 3(a) drop more slowly with increasing |l| than for (b); see appendix B. Therefore, the site of highest ionization probability is shifted more strongly towards negative l. Recombination is most probable for j = 0 in figures 3(a) and (b) which is consistent with previous findings [19]. The drop in probability for increasing j is due to d j , which is why |P jl | extends to larger j in figure 3 The results in figures 3 and 4 are displayed for birth times in the positive field cycle 0 ≤ t b ≤ T 0 /2; the negative half cycle would show the same picture, but mirrored about the x− and y−axis (j, l → −j, −l).
Recall that exact and quasiclassical results do not agree well for harmonic n = 15 in figure 2(a) (Ω = 0.5). The reason is found in figure 3(b); disagreement is due to the point (j, l) = (4, −2) that exhibits unusually high probability. We find that at this point k s is approximately zero, and therewith |H| ≈ 0. Since g ∝ 1/ |H|, this leads to a large value of the prefactor g. This behaviour indicates that the quadratic saddle point expansion is no longer sufficient and the next higher order term(s) must be included. In contrast, agreement for harmonic n = 15 in figure 2(a) for Ω = 1.5 is good. This is consistent with the fact that in figure 4, k s ≈ 0 does not occur in areas of high probability.
Finally, the WQC method hinges on saddle point integration which works well when the exponent is rapidly oscillating. This is fulfilled for wide-band semiconductors with large bandwidth (∆) and in the long wavelength limit. When transitioning to smaller ∆ (dielectrics) and shorter wavelengths, saddle point integration is expected to fail at some point. This will be subject to further research. Also, it is generally possible for transitions involving higher conduction bands to contribute to the harmonic spectrum, but this is beyond the scope of the two band model considered here.

IV. CONCLUSION
In summary, we have shown that the full quantum dynamics driving HHG in wide band materials, such as semiconductors, can be quantitatively explained in terms of quasi-classical trajectory propagation. The physical insight offered by trajectory analysis will prove useful for optimization and design of strong field and attosecond experiments and for the development of novel diagnostic applications of HHG, such as reconstruction of the dipole moment [30]. We believe that our approach presents a versatile tool for investigating open issues in strong field solid state physics, such as the role of noise and manybody effects in strong field processes. Beyond that, quantitatively accurate quasi-classical analysis should be of interest for a wider range of topics in material science. P. B. Corkum acknowledges the support of AROSR grant number FA9550-16-0109. G. Ernotte was supported by the Vanier Canada Graduate Scholarship program. treated here, set F y = F z = 0. To leading order where h, f are minors of |H|. For completeness, we have included time derivatives of the laser field which are however small in the long wavelength limit. As a result the leading order term is |H| = v x (k)f (t ′ , t, k).