Wannier-Bloch approach to localization in high harmonics generation in solids

Emission of high-order harmonics from solids provides a new avenue in attosecond science. On one hand, it allows to investigate fundamental processes of the non-linear response of electrons driven by a strong laser pulse in a periodic crystal lattice. On the other hand, it opens new paths toward efficient attosecond pulse generation, novel imaging of electronic wave functions, and enhancement of high-order harmonic generation (HHG) intensity. A key feature of HHG in a solid (as compared to the well-understood phenomena of HHG in an atomic gas) is the delocalization of the process, whereby an electron ionized from one site in the periodic lattice may recombine with any other. Here, we develop an analytic model, based on the localized Wannier wave functions in the valence band and delocalized Bloch functions in the conduction band. This Wannier-Bloch approach assesses the contributions of individual lattice sites to the HHG process, and hence addresses precisely the question of localization of harmonic emission in solids. We apply this model to investigate HHG in a ZnO crystal for two different orientations, corresponding to wider and narrower valence and conduction bands, respectively. Interestingly, for narrower bands, the HHG process shows significant localization, similar to harmonic generation in atoms. For all cases, the delocalized contributions to HHG emission are highest near the band-gap energy. Our results pave the way to controlling localized contributions to HHG in a solid crystal, with hard to overestimate implications for the emerging area of atto-nanoscience.


I. INTRODUCTION
Recently, the techniques of attosecond science, traditionally applied to atoms and molecules in the gas phase [1], have been extended to solid state [2][3][4][5][6][7][8][9][10][11][12]. A crucial difference between solid and gas targets is the localization of the initial state electron wave function, which is spatially confined in isolated atoms and molecules, but can be delocalized in a solid. The effect of wave-function localization on key aspects of light-solid interaction remains intensely debated. Hence some attosecond experiments [11,12] on photoemission from metal surfaces suggest that the localization of the core-band electrons results in relatively large ionization delays, attributed to transport [13], compared to photoemission from delocalized conduction-band states. Other experiments probing photoemission from the same initial state at different photon energies found that larger ionization delays came from resonant excitation into bulk excited states, rather than from the initial localization of the wave function [14,15].
In this work, we investigate electron localization and the underlying microscopic nonlinear response by focusing on the process of high harmonic generation (HHG) in a crystal solid.
A key feature of HHG in atoms is the recombination of the ionized electron with its parent ion, making it a highly localized process. This localization both dramatically limits HHG efficiency and leads to an exponential decline of HHG yield with increasing ellipticity of laser light (ellipticity creates a drift, which exponentially suppresses the return of ionized electrons to the parent ion [17]).
In contrast, the HHG process in a solid can be delocalized, since an electron ionized from one site of the crystal lattice may recombine with any other. However, little is understood about the specifics of this process. For instance, Ghimire et al. [5] found a much weaker dependence of high harmonic yield on ellipticity in solid ZnO than would be expected for a gaseous medium, suggesting a highly delocalized process. A significantly stronger ellipticity dependence (although still weaker than in atoms) in the same target was subsequently found in a theoretical work [18], which shows a 2-3 orders of magnitude drop in HHG yield for ellipticity of 0.5 (compared to only a factor of 5 drop measured in Ghimire [5]). At the same time, a recent experiment on solid Argon found the same dependence on ellipticity as in gas-phase Argon, suggesting the electron recombines with the same lattice site that it was emitted from [9]. The extent of spatial localization, measured experimentally by ellipticity dependence, is believed to be important for attosecond pulse generation and imaging of the electronic wave functions in the solid state [9].
Here, we investigate the spatial dependence of the HHG process in ZnO by introducing an analytic model which uses localized Wannier wave functions in the valence band and delocalized Bloch functions in the conduction band. Prior seminal work [18][19][20] used delocalized Bloch functions both in the valence and conductions bands, and hence was not able to extract spatial information. In addition to accurately calculating the total HHG yield, this Wannier-Bloch approach allows us to separate the contributions of individual lattice sites to each harmonic, and hence determine the degree of localization of the HHG process as a function of experimental parameters. We find that this localization varies significantly both with the harmonic order and with the orientation of a crystal. Our results point to a possibility of controlling the spatial localization of the HHG process (by varying, for instance, laser ellipticity or crystal orientation), with hard to overestimate implications for attosecond pulse generation, HHG efficiency, imaging of attosecond electron dynamics in condensed matter, and for the emerging area of atto-nanoscience as a whole [21].

II. WANNIER-BLOCH DESCRIPTION OF THE HIGH-ORDER HARMONIC GENERATION IN SOLIDS
Analogous to the three-step model in atoms [22], the HHG in a crystal solid via interband transitions can be described as follows [18]: (i) electron (hole) tunneling excitation from the valence band to the conduction one, (ii) electron (hole) acceleration in the conduction (valence) band, and (iii) electron-hole recombination, resulting in an emission of a high harmonic that is a multiple of the frequency of the driving laser.
In most recent experiments the laser field strength, E 0 , across the lattice constant a is comparable to the band-gap energy E g of a typical semiconductor (E 0 a E g few eV). As a consequence, the field can not be considered as a small perturbation [23]. In our model we therefore assume that this condition is satisfied, but the laser field amplitude is below the damage threshold. In addition, the photon energy of the laser field should be much smaller than the typical bandgap energy. This means we restrict our studies to the photon energies in the MIR domain ( ω 0 ≤ 0.5 eV), which implies that the central-laser wavelength λ 0 is much larger than the typical lattice constant, a. Thereby, the dipole approximation is valid for our description of laser-solid interaction. Since the laser field is linearly polarized in the x direction, we adopt here a one-dimensional description. The Hamiltonian of a single electron in a crystal under the action of a laser field is given by where is the laser-free Hamiltonian, with U (x) the lattice periodic potential. In Eq. (1), U int (x, t) = where E 0 is the electric field peak amplitude, ω 0 the carrier frequency and ϕ CEP the carrierenvelope phase (CEP) of the laser field, while N is the number of laser-period cycles.
Unlike the prior work [18,19], we describe the system within a mixed representation: Wannier states in the valence band and Bloch states in the conduction band. In contrast to the Bloch functions, the Wannier functions are spatially localized "elements" of an L 2 space. In terms of localized wave functions, they provide thereby an analogous insight into HHG mechanism as the usual approach used in atomic and molecular systems. Furthermore with the initial condition a j (0) = δ j,j 0 , i.e. the electron starts the dynamics at the site j 0 .
Here j runs over all atomic sites in the crystal. The Bloch functions of an m-th band (m = v for valence band and m = c for conduction band) have a form where u m,k is a periodic function with the same periodicity as the crystal. The wave functions in Eq. (4) can be equivalently represented by a set of Wannier functions wherew m (k) is a product of a normalisation constant and a phase depending on electron momentum k. It has been shown in [24] that for a 1D lattice thew m are independent of k provided the Wannier functions are real and symmetric under appropriate reflection and fall off exponentially with distance. So, to calculate the emitted harmonics firstly we compute the time-dependent dipole moment where d jc (k) is a dipole transition matrix between Wannier w v,j (k) and Bloch φ c,k states.
The physical meaning of this equation can be summarized as follows: at the observed time, t, the electron previously promoted to the conduction band recombines with the valence band via d jc (k) and emits a photon with an amplitude which depends on the amplitudes a j (t) and a c (k, t). Secondly and similar to Vampa et al. [18], the harmonic emission is obtained by modulus squared of the Fourier transform of Eq. (6) According to Vampa et al. [18], at long-laser wavelengths, i.e. between 1.0 and 5.0 µm, the main contribution to the harmonic spectrum is from interband transitions. Consequently, we neglect the intraband contribution terms w j,k |x|w j,k and φ c,k |x|φ c,k in Eq. (6). Our main task then consists of computing the dipole transition d jc (k) and the transition amplitudes a c (k, t) and a j (t). The dipole moment d jc can be further expressed as a product of two terms: one dependent and one independent of j. First, following [18] we approximate the matrix elements as follows The transition dipole moment from conduction to valence band is then expanded The replacement of x by (x − x j ) in the above formula is justified by the fact that w v,j |x − x j |φ c,k = w v,j |x|φ c,k . In addition, to obtain the coefficients a j (t) and a c (k, t) we employ the time-dependent Schrödinger equation, Eq.
(1), with the wave functions defined in Eq.
We use the tight-binding approximation for the description of the band structure and assume the dispersion relations for the valence/conduction bands where I v and I c are hopping parameters in the valence and conduction bands respectively, a is a lattice constant, E c = E g + 2I c − 2I v and E g is the bandgap energy of the solid material.
The matrix elements for the unperturbed Hamiltonian in both the valence and conduction bands read respectively. Thereby, with the previous definitions and after introducing the wave functions Eq. (3) into Eq. (11), we end up with a system of coupled differential equations for a j (t) and a c (k, t) Here we assume only nearest neighbour hopping in the tight-binding approximation for the valence band (only a j−1 and a j+1 appear in the formula forȧ j ). In solving Eq. (14) we take into consideration dynamics only due to the hopping in the lattice (the first two terms in Eq. (14)) and the laser electric field (the third term). Additionally, we neglect the valence electrons dynamics due to the valence-conduction dipole interaction (the last term in Eq. (14)).
We note from Eq. (15) that the coefficients a c (k, t) are directly related to a j (t), which, in turn, denote the localized place from which the electron will be excited from the valence to the conduction band. This provides a localized picture quite distinct from the one obtained using the Bloch-Bloch approach [18]. By neglecting the last term in Eq. (14) and solving it explicitly following [25], we obtain where q runs over all atomic sites, a is lattice constant, J q are Bessel functions of q-th order and We assume that the electron is initially localized at one atomic site j 0 , i.e. a q (0) = δ qj 0 .
Later, due to the interatomic hopping and the acceleration due to the laser electric field, the electron's wave function spreads in the lattice following Eq. (16). The width of the electron's wave function spread in the lattice at the end of the laser pulse depends on the hopping amplitude I v , the lattice constant a, the laser electric field strength E 0 and pulse duration. We allow all coefficients a j (t) to acquire non-zero values during the laser pulse.
Equation (15) is solved by following [18,26] and thereby we write Where, ϕ(p, t, t ) = t t E c (p + A(t ))dt is the accumulated phase of the electron in the conduction band. Equation (22) describes the harmonic emission originating from a single electron in a lattice. Its interpretation is similar to harmonic emission in atoms. In particular, the dipole radiation contains all the "relative trajectory contributions" of the electron wave function to the emitted harmonics from a solid. As depicted in Fig. 1, first, the electron located at the j -th atomic site is excited from the valence to the conduction band. Second, it is accelerated within the conduction band by the laser field. Finally, the electron has some probability of recombining with a different atom, located at j-th site (see the arrow pointing down in Fig. 1). As a result, excess electron energy is emitted in the form of a high harmonic of the driving laser frequency. Here it is assumed that the electron is initially localized at j 0 atomic site, as was mentioned above.
To account for contribution of all electrons in a lattice, we multiply the dipole moment from the shift j → j + n, q → q + n in a * j (t)a j (t ). Since these two phases cancel each other, the contribution from each electron in the system to the harmonic emission will be exactly the same.

A. Comparison of Bloch-Bloch and Wannier-Bloch models
We begin by comparing the emitted HHG spectrum from the Wannier-Bloch approach, Eq. (22), with the spectrum obtained from the delocalized Bloch-Bloch approach. The Bloch-Bloch model was implemented following [18], with the cosine band structure approximated by a Taylor expansion up to the fourth order and integration in momentum space replaced by a saddle-point approximation. Integration over the ionization time was done numerically using a Gaussian quadrature routine and the Fourier integral was performed as an FFT. Our approach exhibits good agreement with the Bloch-Bloch model, reproducing the plateau, the cutoff and the standard odd harmonic structure. The two spectra differ mainly in the low-order harmonics region, suggesting that localization (or recombination with the parent atom) may play greater importance in the production of low-order harmonics, as is confirmed in Fig. 2(b). Overall, this comparison confirms that the Wannier-Bloch picture reproduces the essential features of the Bloch-Bloch model for the emitted harmonics. To calculate the contribution of a given ∆j, we apply FFT on the part of the dipole d(t) composed only of terms in Eq. (22) for which |j − j | is equal ∆j. The relative length of displacement between the ionization and recombination sites is given by D s = ∆j a, which we call the length of Electron-recombination at a Relative Atom-center (ERA).
As Fig. 2(b) shows, in certain parts of the spectrum, even ∆j = 10 paths contribute considerably to the total harmonic emission. This is in clear contract to HHG in atomic gas, where the electron has to recombine with its parent atom, corresponding to ∆j = 0 contributions only. In the next section, we attempt to understand why even relatively distant atomic sites can contribute significantly to the total emission spectrum in solids.

B. Wannier-Bloch picture
In order to further investigate the contribution of different ERA to the HHG spectrum, we calculate harmonic emission for a set of different laser-parameters. This also allows us to establish under what conditions the Wannier-Bloch approach may be a more adequate description relative to the Bloch-Bloch one [18]. Due to computational constraints in the HHG spectra calculations using the band-structure of Sec. IIIA, here we focus on the narrower valence band case. For the latter we are able to scan a wider range of parameters and analyze the band structure influence on the different ERA contributions to the HHG.
In order to compute the HHG spectra, we fix the optical axis of ZnO (with polarization of the laser in the Γ−A direction) [18] and use: a = 9.83 a.u., I c = 0.02175 a.u., I v = −0.00295 a.u. and E g = 0.1213 a.u. Also, following the formulation in [18], d vc (k) = i.e. for the case when the electron recombines to the same atomic-parent site from where it was previously excited to the conduction band. The longer the electron recombination displacement D s is, the lower is its contribution to the HHG spectra. Both panels of Fig. 3 show that, while the harmonics of odd orders decay very fast with ∆j, there is a relatively large signal between the 8-th and 15-th harmonic which is preserved also for larger values of ∆j. This signal corresponds to the energy gap between valence and conduction bands, which spreads from about ∼ 8ω 0 (band gap for k = 0) to 14.5ω 0 (maximum energy gap for k = π a ). The greater contribution from large ∆j processes near the maximum and minimum energy gaps can be understood as resulting from the high density of possible interband transitions involving opposite band edges. Those transitions occur between states with narrowly defined momenta (near the band extrema) which require broad spatial coherence as reflected in large recombination lengths. Figure 3 shows the typical features of HHG spectra: namely odd harmonics are present, the signal is strongest for the low-order harmonics (1st, 3rd), and a plateau region and a cutoff can be easily distinguished. As would be expected for interband emission, the cutoff is located near the harmonic equivalent of the maximum energy difference between the conduction and valence bands. The region of the plateau exhibits pronounced interference structures -there is no clear recognition of even-odd harmonic symmetry. This behavior is also typical of the harmonic spectra in atoms in the limit of short pulses [28,29].
The spectra in Fig 3 can be compared with the experimental results shown in Fig. 3 of Ref. [5], specifically for the crystal angle 0 • . In this experimental data, up to the 13-th harmonic is distinguishable, which is close to the cut-off value obtained in our calculations.
As predicted by our model, a signal near the bandgap energy is observed in Ref. [5] (where, however, it is attributed to fluorescence effects). In contrast to our results, in the experiment both odd and even harmonics appear in the spectrum. This is likely an effect of a symmetry breaking in the 3D-ZnO lattice, which we do not take into account here. Laser intensity was set to I 0 = 5 × 10 11 W/cm 2 , laser pulse length to 10 laser periods and ϕ = 0 rad.
To investigate how harmonic emission scales as a function of wavelength, we calculate the harmonic spectrum for different values of the laser wavelength, λ 0 . Note that the saddle point method was not used to solve the momentum integral in Eq. (22). The full integration is a reasonable choice while the energy variation of the dipole-radiation phase is comparable to the driving field frequency [30]. The argument is based on the fact that the dipole-radiation phase takes into account electron propagation in the conduction band.
Nevertheless, for larger energy band-structure, the saddle point method is suitable for calculating the momentum integral in Eq. (22).
The results are shown in Fig. 4. Plot (a) shows HHG spectra for wavelengths in a range of 800-3500 nm. The frequency axis for each wavelength is scaled in ω 0 = 2πc/λ 0 units.
It is observed that the cut-off moves to lower harmonics while the wavelength decreases.
However, one may expect that the cut-off stays constant in terms of photon energy because of a well-defined maximum energy difference between the conduction and valence bands of about 6 eV. This effect is shown on Fig. 4 (b), where the spectra of Fig. 4 (a) are replotted as a function of photon energy. Hence, the cut-off of the harmonic spectrum is in good agreement with the value calculated from the band structure. From the so-called "action phase" in Eq. (22), i.e. ϕ(p, t, t ), we can infer that the maximum harmonic energy produced in a solid lattice should be limited by the band dispersion relation (this result is consistent with prior findings [18]).
Previously, it was found that HHG in solids scales linearly with the electric field strength [30]. Thereby, to further test the Wannier-Bloch model, we calculate HHG spectra for different values of the laser electric field. Figure 5 shows HHG spectrum for electric field amplitudes, E 0 , in the range between 0.4-1.0E 0 , where E 0 = 3.779 × 10 −3 a.u. corresponds to laser intensity I 0 = 5 × 10 11 W/cm 2 . As expected, decreasing the laser intensity results in shifting the cutoff to lower order harmonics. This can be seen more clearly in Fig. 6 where the spectral intensity of the signal is plotted for all odd harmonics at each laser intensity.
The fast decrease of the harmonic yield for decreasing electric field strength at low E 0 is a typical low-field behavior in transitions induced by electric fields such as e.g. interband Zener transitions. The decrease in cutoff with decreasing electric field strength was demonstrated experimentally in [30]. Although, it is difficult to determine the exact position of the cutoff due to the few harmonics present for the lowest intensity cases, we trace in For providing additional key information on localization of HHG in solids, our results agree in many aspects with prior experimental and theoretical observations [5,8,18,19]. In particular, (i) the HHG cutoff shows a dependence on the maximum energy difference between the valence and conduction bands as well as on the laser wavelength and peak intensity, (ii) for long laser pulses and few-cycle laser fields the model depicts the full odd and a continuum spectra, respectively, and (iii) we find a direct link between the emitted harmonic spectrum shape and the band-structure.

IV. CONCLUDING REMARKS
By using localized atomic sites in the valence band and delocalized functions in the conduction band, our model has the closest parallels to harmonic generation from atomic gas. As such, it allows one to access contributions of individual lattice sites, and hence assess the degree of localization of HHG in solids -something that has previously been inaccessible.
In particular, we can describe a process in which an electron initially localized at the j -th atom in the valence band has a certain probability to be excited to the conduction band, where it is accelerated to a high energy before recombining either to the parent atom, at j -th site, or (with different probability) to any other j-th atom in the lattice.
Different displacements of the electron recombination atomic-sites, i.e. ∆j = |j − j |, give different contributions to the harmonic spectrum. The approach developed here allows to extract all of these contributions. In particular, the main contribution was found to be given by ∆j = 0, or electron recombining at the same atomic site it was excited from. Especially for the case of narrow bands in the band structure, lower ∆j contribute by far the most to the harmonic spectrum, signifying substantial localization in the HHG process. On the other hand, we found enhanced contribution of high ∆j in case of wider valence and conduction bands. This enhanced delocalization is likely due to increased mobility of the electrons in a lattice. In all cases, distant neighbor contributions were highest near the band-gap energy.
This suggests that harmonic yield near the band gap energy should decline less (relative to other harmonics) with increasing ellipticity of laser light, since elliptical polarization tends to suppress local contributions.
Note that our results by means of the Wannier-Bloch approach demonstrate a different interpretation than those presented by the Bloch-Bloch picture in ref. [18]. While our model recreates the conventional atomic picture, the Bloch-Bloch model is based on the electronhole pair recombination. Thereby, it is clear that both approaches, while predicting similar total HHG spectra, provide different interpretations of the HHG process.
Our results suggest that it should be possible to control the localization of the HHG process by varying experimental parameters. Hence, by quantifying site specific contributions, our work paves the way to controlling HHG efficiency, creation of attosecond pulses and imaging of the electronic wave function in a crystal lattice [9].