## Abstract

We analyze high-order harmonic generation (HHG) in a disordered semiconductor within the context of the Anderson model of disorder. Employing the theoretical methods pioneered for the study of disordered metals, we show that disorder is a source of ultrafast dephasing of the HHG signal in semiconductors. Furthermore, it is shown that the dephasing effect induced by disorder on HHG spectra depends on both strength and correlation length of the disorder and very weakly on the frequency and intensity of the laser. Our results suggest that HHG has the potential to be a new spectroscopic tool for the analysis of disordered solids.

© 2020 Optical Society of America

## 1. INTRODUCTION

In a typical solid high-order harmonic generation (HHG) experiment, a wide bandgap semiconductor [1–5], or insulator [6–8], is made to interact with an intense laser field (approximately $1\;{\rm TW}/{{\rm cm}^2}$) with mid-infrared wavelength (2–5 µm). Driven by the external field, the crystal emits characteristic electromagnetic spectra containing many well-resolved peaks in correspondence to the harmonics of the laser frequency. It has already be demonstrated that solid HHG has many useful applications, such as the optical reconstruction of the band structure [9,10] or the investigation of ultrafast electronic processes in solids [11]. Furthermore, the potential of HHG for the study of ferromagnetic materials [12], Mott insulators [13], and other exotic materials [14] has been carefully investigated. It is therefore not surprising that solid HHG has attracted much attention in recent years. Nevertheless, as yet, a complete physical understanding of solid HHG has not been obtained.

Actually, it was shown in [4,15] that simulations based on the semiconductor Bloch equations (SBEs) [16] can reproduce the main features of HHG experiments. Nevertheless, to do so, it is required to introduce in the SBE a semi-phenomenological dephasing time of the order of 1 fs. These ultrashort dephasing times are problematic because a physical mechanism that can produce them has not clearly been identified [17]. It has been conclusively proven that Coulomb interactions cannot be responsible for the ultrashort dephasing times in HHG. In fact, simulations based either on the time-dependent density functional approach (TDDFT) [18–20] or the Hartree–Fock approximation [21] show that, even in the presence of Coulomb interactions, the simulated HHG spectra do not show clean harmonic peaks resembling those observed in experiments. Furthermore, these simulations have confirmed the validity of the single-particle picture in the analysis of HHG from solids. For this reason, HHG has been almost universally modeled [22–31] as the interaction of a set of independent electrons with a spatially uniform and time-dependent electric field and a time-independent and spatially periodic crystalline field.

A few papers have of late analyzed how the lack of perfect translation symmetry in the crystal [32–41] may affect the HHG process in solids. These studies show beyond doubt that the HHG process is strongly sensitive to the presence of impurities or dopants in the crystal. In particular, in Refs. [32,33], we have pointed out that disordered crystals can emit HHG spectra that strongly resemble those observed in experiments. In fact, a model disordered crystal possesses HHG radiation composed of well-defined odd harmonics of the fundamental laser frequency. In the past work, we have used only a simple heuristic argument to explain how disorder may affect the HHG process in solids. In this paper, we extend the results of our previous work and, adopting the tools developed for the analysis of disordered metals [42–46], show that disorder is a source of ultrafast dephasing for HHG in solids. The calculations reveal that the dephasing effect induced by disorder depends on both strength and correlation length of the disorder and, very weakly, on the intensity and frequency of the laser field. This result suggests that HHG could be used in the future as a source of detailed information about disordered solids.

## 2. TWO-BAND DISORDERED MODEL

We consider the interaction of a semiconductor with a strong short laser pulse. As reported in the introduction, in the study of solid HHG, it is acceptable to simulate the semiconductor as an ensemble of noninteracting electrons characterized by a one-particle Hamiltonian ${\hat H_0}$ and a one-particle density matrix $\rho$ [18]. To further simplify the problem, we only consider a one-dimensional two-band tight-binding Hamiltonian. The generalization to a more complex situation is straightforward.

The single-particle Hamiltonian of the electron in the solid is ${\hat H_0} = {\hat H_I} + {\hat H_D}$, where ${\hat H_I}$ is the Hamiltonian of an ideal crystal that possesses a perfect crystalline structure, and ${\hat H_D}$ is a small correction that models the effects of impurities and defects in the semiconductor. Specifically, we assume that the ideal crystal Hamiltonian ${\hat H_I}$ has an energy spectrum made of only one valence and one conduction band,

In the presence of a linearly polarized spatially uniform laser field $F(t)$, the crystal Hamiltonian can be written as $\hat H(t) = {\hat H_0} + \textit{eF}(t)x$ or in Wannier representation as

The density matrix of the solid satisfies the Liouville equation, $i\hbar \frac{{d\rho (t)}}{{dt}} = [\hat H(t),\rho (t)]$, with the formal solution, $\rho (t) = \hat U\rho (0){\hat U^\dagger }$, and the evolution operator $\hat U$ solves the Schrödinger equation $i\hbar \frac{{d\hat U}}{{dt}} = \hat H(t)\hat U$. For $t \gt 0$, the density matrix of the solid has the form $\rho (t) = \sum\nolimits_{m = - N}^N {c_0}|m(t)\rangle \langle m(t)|$, where the vector $|m(t)\rangle$ is the solution of the one-particle Schrödinger equation,

with initial condition, $|m(t = 0)\rangle = |{w_\textit{mv}}\rangle$. We can expand $|m(t)\rangle$ as $|m(t)\rangle = \sum\nolimits_{n = - N}^N [b_c^m(n,t)|{w_\textit{nc}}\rangle + b_v^m(n,t)|{w_\textit{nv}}\rangle ]$, where the time-dependent coefficients $b_c^m(n,t),b_v^m(n,t)$ are solutions of the coupled differential equations,Once the one-particle density matrix is known, it is possible to calculate any one-particle observable. In particular, we can write the current as

In absence of disorder, all the matrix elements $\langle m(t)|\hat p|m(t)\rangle$ are identical. When disorder is present, the $\langle m(t)|\hat p|m(t)\rangle$ are in general different, and a sum over all the possible initial conditions in Eq. (5) is needed to evaluate the current. For each different initial condition, the ket $|m(t)\rangle$ describes an electron that, at $t = 0$, is localized on the lattice site m. This particle evolves under the action of the laser field, the crystalline field, and a given disorder potential (characterized by a given realization of the ${\alpha _n}$, for $n = - N, \ldots ,N$). We note that summing over the $2 N+ {1}$ initial conditions for a given configuration of disorder (a single realization of the ${\alpha _n}$ for $n = - N, \ldots ,N$ and $2 N+ {1}$ different initial conditions $|m(t = 0)\rangle = |{w_\textit{mv}}\rangle$, $m = - N \ldots N$) is equivalent to summing over $2 N+ {1}$ different disorder configurations for a single initial condition (a fixed m, for example $m = 0$, and $2 N+ {1}$ disorder configurations characterized by $2 N+ {1}$ different sets of $\alpha _n^m$ that satisfy the condition $\alpha _n^m = {\alpha _{n - m}}$ for $m = - N\ldots N$). Since the number of sites $2 N+ {1}$ is a macroscopic number in a crystal, we make the assumption that the sum over $2 N+ {1}$ initial conditions is equivalent to that over all the possible configurations of disorder.

## 3. KELDYSH APPROXIMATION AND LONG-WAVELENGTH LIMIT

In Keldysh approximation [50] and letting the number of sites $2N + 1$ in the crystal go to infinity, we can write the equation of motion [Eq. (3)] of the electron in the two-band model in quasi-momentum space as

The long-wavelength approximation is based on the following observation. Since the lasers employed in typical solid HHG experiments have only moderate intensities (${\simeq} 1\;{\rm TW/cm}^2$ [1]), only electrons with small quasi-momentum can tunnel from the valence to the conduction band, thus contributing to the HHG process [22]. For this reason, we can approximate the $\cos$ functions in Eq. (6) as quadratic polynomials ($\cos (ap) \simeq 1 - {a^2}{p^2}/2$) and formally consider the quasi-momentum $p$ as defined on the whole real axis ($p \in ( - \infty , + \infty )$). We then define two continuous functions ${\phi _c}(x,t),{\phi _v}(x,t)$ that satisfy the Scrhödinger-like equations of motion,

As explained earlier, to simulate the effect of disorder in HHG experiments, the equations of motion Eq. (7) are solved for the initial same conditions ${\phi _c}(x,t = 0) = 0,{\phi _v}(x,t = 0) = \delta (x)$ and different configurations of disorder [different realization of the function D(x)]. Within the long-wavelength approximation, we then approximate the interband dipole moment for a single realization of disorder $d(t) = {d_\textit{cv}}\sum\nolimits_n (b_\textit{nc}^*(t){b_\textit{nv}}(t) + {\rm c.c})$ as $d(t) = ({d_\textit{cv}}/a)\int {\rm d}x(\phi _c^*(x,t){\phi _v}(x,t) + {\rm c.c})$. The total dipole moment of the system is the sum over all the configurations of disorder that we write formally as $\bar d(t) = ({d_\textit{cv}}/a)$$\int {\rm d}x(\bar \phi _c^*(x,t){\phi _v}(x,t) + {\rm c.c})$. ${\bar \phi _c}(x,t)$ is the average of the function ${\phi _c}(x,t)$ over the different configurations of disorder.

The equation for ${\phi _c}$ can be solved using the Green function approach [47,52],

In order to describe the role of disorder in HHG, we need to find the average propagator ${\bar G_c}$, which is obtained by integrating over all possible configurations of disorder [42–47] ${\bar G_c}(x,x^\prime ,t,t^\prime ) \;\;=\;\; \Theta (t - t^\prime )\int [dD]\int_{x(t^\prime ) = x^\prime }^{x(t) = x} [{\rm d}x]{e^{\frac{i}{\hbar }S\big( {[x]_{tt^\prime }^{xx^\prime }} \big)}}$. The functional integration of the disorder configuration is represented in the formula by the formal expression $\int [{\rm d}D]$. In particular, we consider a disorder potential of the form $D(x) = {W_0}\sum\nolimits_{i = 1}^{{N_D}} v(x - {r_i})$, where $v(x)$ is a dimensionless function, ${W_0}$ characterizes the strength of the disorder, and ${r_i}$ are random variables. This is the model of disorder introduced by Edwards and Gulyaev [45,47]. The function D describes the potential generated by ${N_D}$ impurities randomly positioned in space. The random variables ${r_i}$ are assumed to be described by a uniform probability distribution. It is convenient to choose $v(x) = {e^{ - {x^2}/{L^2}}}$. With the above choices, the functional integral over disorder configurations, equivalent to an integral over the variables ${r_i}$, can be exactly calculated, and the result is for weak and dense disorder (small ${W_0}$, large ${N_D}$) [42–46] (and from now always for $t \gt t^\prime $),

The average propagator can be rewritten as

The functional integral above cannot usually be performed in closed form, and it is usually estimated in the so-called first cumulant approximation [42,52–54] ${ \langle }{e^{\Phi (x,x^\prime ,t,t^\prime )}} \rangle \simeq {e^{\bar \Phi (x,x^\prime ,t,t^\prime )}}$, where

#### A. Dipole Moment

The interband part of the dipole moment $\bar d(t)$ averaged over the disorder configurations is written in term of the Green functions as

The Fourier spectrum of the HHG radiation $|P(\omega )|^2$ is related to Fourier transform of $\bar d(t)$ via the relation $P(\omega ) = \int {e^{i\omega t}}\dot {\bar d}(t){\rm d}t$. Using the expression in Eq. (13) for $\bar d(t)$, the Pauli–Van Vleck form of the propagators [42] (${\bar G_c} = A{e^{\frac{i}{\hbar }(S_c^\textit{cl} - i\sigma )}},{G_v} = A{e^{\frac{i}{\hbar }(S_v^\textit{cl} + { \epsilon _G}t)}}$) and adopting the stationary-phase approximation in calculating the integrals in Eq. (13), we find that the Feynman trajectories that significantly contribute to the HHG process must satisfy the conditions [15,57–62]

Since $\bar \Phi$ depends only on the difference $x - x^\prime $, we can combine the last two conditions into a single relation,

In terms of the trajectories for the valence-band electron ${e_v}$, valence-band hole ${\bar e_v}$, and conduction-band electron ${e_c}$ introduced in the discussion following Eq. (13), Eq. (15) stipulates that the Feynman trajectories contributing to HHG are those describing a valence-band electron ${e_v}$ and a valence-band hole ${\bar e_v}$ propagating along the same trajectory from the initial position $x = 0$ at time $t = 0$ to the position $x = x^\prime $ at $t = t^\prime $, where the valence-band electron is destroyed and a conduction-band electron ${e_c}$ is created. Later on, in the time interval $(t^\prime ,t)$, ${e_c}$ and ${\bar e_v}$ propagate along their respective trajectories until they meet each other again at the position $x$ at time $t$. Since in the time interval $(0,t^\prime )$ the valence-band electron and valence-band hole move along the same path, the Feynman trajectories relevant for HHG are those describing an electron–hole pair created at time $t^\prime $ at $x^\prime $ and recolliding at time $t$ at $x$.

The first two conditions in Eq. (14) enforce energy conservation at the moment of the creation $t^\prime $ and destruction $t$ of the particle–antiparticle system. In particular, the sum of Eqs. (14b) and (14a) tells us that to emit a photon of energy $\hbar \omega$, the electron–hole pair must acquire from the external field, on average, a kinetic energy equal to $\hbar \omega - (\frac{{\partial {\sigma _I}}}{{\partial t}} + \frac{{\partial {\sigma _I}}}{{\partial t^\prime }})$. The term $\frac{{\partial {\sigma _I}}}{{\partial t}} + \frac{{\partial {\sigma _I}}}{{\partial t^\prime }}$ accounts for the fact that the electron, on average, will dissipate energy in inelastic collisions while traveling in the disordered medium. Equation (14b) mandates that at $t^\prime $ energy of the valence-band electron [63], ${E_v}(t^\prime ) = - \frac{{\partial S_v^\textit{cl}}}{{\partial t^\prime }}([x]_{t^\prime 0}^{x^\prime 0})$ differs from the energy of the electron ${E_c}(t^\prime ) = \frac{{\partial S_c^\textit{cl}}}{{\partial t^\prime }}([x]_{tt^\prime }^{xx^\prime })$ by ${\tilde \epsilon _G} = { \epsilon _G} - \frac{{\partial {\sigma _I}}}{{\partial t^\prime }}$ or, more explicitly, $\Delta E(t^\prime ) = \frac{{\partial S_v^\textit{cl}}}{{\partial t^\prime }}([x]_{t^\prime 0}^{x^\prime 0}) + \frac{{\partial S_c^\textit{cl}}}{{\partial t^\prime }}([x]_{tt^\prime }^{xx^\prime }) = - {\tilde \epsilon _G}$. We observe that for the trajectories contributing to HHG, those satisfying Eq. (15), the energy of the valence-band hole at ${t_1}$ is the opposite of that of the valence-band electron ${E_v}(t^\prime )$. The quantity $\Delta E(t^\prime )$ can therefore be considered as the total energy of the electron–hole pair at time ${t_1}$. We can write the generic electron–hole trajectories ${x_{{e_c}}},{x_{{{\bar e}_v}}}$ for the pair created at $x^\prime $ at $t^\prime $ and recolliding at $x$ at $t$ as

Finally, Eq. (14d) requires that, at time $t^\prime $ when the electron–hole pair is created, the momentum of the conduction-band electron ${p^c}([x]_{tt^\prime }^{xx^\prime },\tau = t^\prime ) = - \frac{{\partial S_c^\textit{cl}}}{{\partial x^\prime }}([x]_{tt^\prime }^{xx^\prime })$ differs from the momentum of the valence-band electron ${p^v}([x]_{t^\prime 0}^{x^\prime 0},\tau = t^\prime ) = \frac{{\partial S_v^\textit{cl}}}{{\partial x^\prime }}([x]_{t^\prime 0}^{x^\prime 0})$ by the amount $\frac{{\partial {\sigma _I}}}{{\partial x^\prime }}([x]_{tt^\prime }^{xx^\prime })$,

Equation (17) defines the distance $(x - x^\prime )$ covered by the electron–hole pair before annihilation as an implicit function of the recollision time $(t - t^\prime )$ and of the parameters of the laser field and disorder. It is not possible to write explicitly the solutions of Eq. (17), but it is possible to see that there are solutions with small travel distance $(x - x^\prime \lt L)$, when the correlation length is large ($L \gg {a_0}$) and laser strength and frequency have the values commonly used in HHG experiments. In fact, usually in HHG experiments, the ratio ${F_0}/\omega _L^2$ is of the order of ${a_0}$ (with ${F_0}$ the laser peak strength and ${\omega _L}$ its frequency); therefore, for $(x - x^\prime ) \ll L$, we can set ${e^{ - A}} \simeq 1$ [Eq. (10)] and write

In summary, the Feynman trajectories that contribute to the HHG process in our model disordered semiconductor are almost identical to those of a defect-free solid [15,58]. Nevertheless, in the HHG process in disordered solids, the contribution of each trajectory to the HHG signal is weighted by a factor ${e^{{\sigma _R}(x \simeq x^\prime ,t,t^\prime )}}$. This factor, calculated in Eq. (11), strongly suppresses the contribution of trajectories with long travel times (large $t - t^\prime $) [32]. It has been observed in Refs. [15,32,58] that to obtain a clean HHG spectrum it is required to eliminate the contribution of all the trajectories with travel time longer than a quarter of the laser period ${T_L}/4 = \pi /2{\omega _L}$. In our model, this condition is met if the parameters that characterize the disorder, $\rho ,{W_0},L$, are such that ${-}{\sigma _R}(x = x^\prime ,t - t^\prime = {T_L}/4) \ge 1$. When only short trajectories contribute to HHG, the current $J(t)$ satisfies the relation [64] $J(t) \simeq - J(t - {T_L}/2)$ for $t$ in the generic ${m}$th half-cycle, $t \in [m{T_L},(m + 1/2){T_L}]$. In fact, the current is generated by electron–hole quantum trajectories with short recombination times. These short trajectories in adjacent half-cycles are almost identical because the laser envelope does not change significantly over such a short time interval. As a result, the Fourier transform of the current vanishes for $\omega = 2n{\omega _L}$ (with $n$ any integer), and the spectrum $|J(\omega )|^2$ is made up of odd harmonics of the laser frequency.

## 4. NUMERICAL RESULTS

We pause at this point to stress that the main finding of this paper is that the HHG spectrum generated by a disordered crystal is made up only of peaks in correspondence of, or very close to, the harmonics of the laser frequency ${\omega _L}$, when a sufficient amount of disorder is present. This result is important because all the reported experimental HHG spectra are indeed built of a series of sharp peaks at frequencies that are integer multiples (or very close to integer multiples), i.e., harmonics, of the laser frequency. This fact is clearly shown in the experimental papers Refs. [1–10]

We further stress that *ab initio* numerical
calculations based on model Hamiltonian describing a perfect crystal,
i.e., a crystal without impurities or defects, are not able to reproduce
the experimental HHG spectra. As a support to this statement, we mention,
just as a few examples among many others, the results in Ref. [21], where the spectra are calculated
using the Hartree–Fock approach, the numerical computations in
Refs. [17,20] with spectra obtained within the Kohn–Sham density
functional scheme, or the findings in Refs. [15,29] based on
the use of the SBE with or without dephasing time. All these simulations
show HHG spectra with both harmonic peaks and a large number of
nonharmonic peaks at variance with the experiments. The only simulations,
other than those based on the disordered Hamiltonian used in our paper and
in our previous works [32,33], that produce HHG spectra that
resemble the experimental ones are those based on the SBE with ultrashort
dephasing time (dephasing time of the order of 1 fs). Nevertheless,
simulations based on the SBE with dephasing time cannot considered as
*ab initio* simulations. In fact, the SBE with
a dephasing time cannot be derived from a microscopic many-body
Schrödinger equation (see, for example, [16,65,66]). The SBEs with a dephasing time are
semi-phenomenological equations where the effects such as the lattice
vibrations (phonons), the Coulomb interaction beyond the Hartree–Fock
approximation, or disorder are approximately treated with the introduction
of an arbitrary parameter (the dephasing time) whose value is fixed by
trying to match the simulated value of some observable quantity to its
experimental counterpart as closely as possible. In the case of HHG, the
value of dephasing time is chosen to reproduce the emitted spectra. The
SBEs with dephasing time are certainly very useful; nevertheless, they
cannot reveal which physical phenomenon is responsible for the fact that
the experimental HHG spectra are composed of harmonic peaks only.

We now compare the predictions of the long-wavelength model against the numerical solutions of the equation of motions [Eq. (3)]. In the numerical calculations, the external field has the form

In Fig. 2, we plot the spectrum emitted by a disordered semiconductor interacting with a field of intensity $I = 2.0\;{\rm TW}/{{\rm cm}^2}$ and wavelength $\lambda = 2.5\;\unicode{x00B5}{\rm m}$. The strength of the disorder potential is ${\Gamma _1} = 0.25\hbar {\omega _L}$. We observe that, in this case, the HHG spectrum presents both harmonic peaks and nonharmonic peaks. In detail, we note that in the HHG spectrum between the bandgap frequency (approximately $8{\omega _L}$) and the cutoff (approximately at $15{\omega _L}$), there are three odd harmonic peaks at $11{\omega _L}$, $13{\omega _L}$, $15{\omega _L}$; one even harmonic peak at $12\hbar {\omega _L}$; and 5 nonharmonic peaks. This result is in agreement with the long-wavelength model since ${-}{\sigma _I}(x = x^\prime ,T/4) \ll 1$. In fact, such a small value of ${\sigma _I}$ means that the disorder is too weak to effectively suppress the contributions of trajectories with long recollision times to the dipole moment. As a result, the HHG spectrum above the bandgap frequency contains a large number of nonharmonic peaks.

In Fig. 3, we present the spectrum radiated by a disordered semiconductor with ${\Gamma _1} = 0.70\hbar {\omega _L}$, with all other parameters of the semiconductor and of the laser field being the same as in the previous calculation. We note that, for the given value of disorder, the long-wavelength model predicts that ${-}{\sigma _I}(x = x^\prime ,T/4) \simeq 1$. According to the discussion in the previous sections, this means that the disordered solid should radiate an HHG spectrum composed only of harmonic peaks. In fact, as we have shown already in Refs. [32,33], the HHG spectrum of the disordered solid is for all frequencies above the bandgap (${\simeq} 8{\omega _L}$) much weaker than that generated by the perfect solid, and it is made almost entirely of odd harmonics of the laser frequency. In particular, we note that, in the region [$8{\omega _L}$-$15{\omega _L}$], there are now four harmonic peaks ($9{\omega _L}$, $11{\omega _L}$, $13{\omega _L}$, $15{\omega _L}$), and only one nonharmonic peak close to $8{\omega _L}$.

To further illustrate the role of disorder in the ultrafast optical response of a crystal, we plot in Fig. 4 HHG spectra, in the range between the bandgap frequency and the cutoff, as a function of the amount disorder in the crystal. In particular, the results presented in this figure support our claim that the role of disorder in HHG is well described through the long-wavelength model. In fact, it is clearly shown in this figure that, only when the condition ${\sigma _R}(x = x^\prime ,t - t^\prime = \pi /2{\omega _L}) \ge 1$ is met (for disorder strength ${\Gamma _1} = 0.7\hbar {\omega _L}$), the HHG spectra are made up, almost entirely, of odd harmonics of the laser frequency. It is useful to stress that, even in experimental HHG spectra (for example [1,3]), some nonharmonic peaks are visible close to the bandgap frequency.

In conclusion, we have investigated HHG from an imperfect semiconductor, and we have observed that the HHG spectra generated resemble those observed in experiments when sufficiently strong disorder is present in the solid. Further, a simple theory has been presented that explains how the impurities in the crystal can generate the ultrafast dephasing effect that is responsible for the clean and well-resolved HHG spectra in the numerical simulations. This dephasing effect depends critically on the strength and correlation length of the disorder. Our results therefore show that (1) disorder is a source of ultrafast dephasing of the HHG signals, and (2) it may be responsible for the well-resolved experimental HHG spectra. Furthermore, our results suggest that it may be possible to utilize HHG as a novel experimental tool for the analysis of disordered solids.

## APPENDIX A: DERIVATION OF $\bar \Phi$

The starting point for the calculation of $\bar \Phi$ is Eq. (9) in the main text. Since we have chosen $v(x)$ a Gaussian function, the correlation function is a Gaussian $\gamma (x(\tau ^\prime ) - x(\tau )) = \sqrt \pi L{e^{ - {{(x(\tau ^\prime ) - x(\tau ))}^2}/2{L^2}}}$, where $L$ is the parameter that determines the length-scale of the correlation of the disorder. The first cumulant can now be written as

*cl*is used to remind that the actions $S_k^\textit{cl},S_c^\textit{cl}$ are calculated for the trajectory satisfying the classical Newton equation of motions with the boundary conditions $x(\tau = t^\prime ) = x^\prime ,x(\tau = t) = x$. The first cumulant approximation can therefore be written as

## Funding

Ministry of Science and Technology, Taiwan; National Taiwan University (108L893201 and 108L104048).

## Disclosures

The authors declare no conflict of interests.

## REFERENCES

**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. **7**,
138–141 (2011). [CrossRef]

**2. **G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, “Linking high harmonics from
gases and solids,” Nature **522**, 462–464
(2015). [CrossRef]

**3. **S. Gholam-Mirzaei, J. Beetar, and M. Chini, “High harmonic generation in
ZnO with a high-power mid-IR OPA,” Appl. Phys.
Lett. **110**, 061101
(2017). [CrossRef]

**4. **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 terahertz
high-harmonic generation by dynamical Bloch
oscillations,” Nat. Photonics **8**, 119–123
(2014). [CrossRef]

**5. **P. Xia, C. Kim, F. Lu, T. Kanai, H. Akiyama, J. Itatani, and N. Ishii, “Nonlinear propagation effects
in high-harmonic generation in reflection and transmission from
gallium arsenide,” Opt. Express **26**, 29393–29400
(2018). [CrossRef]

**6. **T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, “Extreme ultraviolet
high-harmonic spectroscopy of solids,”
Nature **521**,
498–502 (2015). [CrossRef]

**7. **M. Garg, H. Y. Kim, and E. Goulielmakis, “Ultimate waveform
reproducibility of extreme ultraviolet pulse by high-harmonic
generation in quartz,” Nat. Photonics **12**, 291–296
(2018). [CrossRef]

**8. **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 **534**, 520–524
(2016). [CrossRef]

**9. **G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, D. D. Klug, and P. B. Corkum, “All-optical reconstruction of
crystal band structure,” Phys. Rev.
Lett. **115**, 193603
(2015). [CrossRef]

**10. **A. A. Lanin, E. A. Stepanov, A. B. Fedotov, and A. M. Zheltykov, “Mapping the electron band
structure by intraband high-harmonic generation in
solids,” Optica **4**,
516–519 (2017). [CrossRef]

**11. **M. Hohenleutner, F. Langer, O. Schubert, M. Knorr, U. Huttner, S. W. Koch, and R. Huber, “Real-time observation of
interfering crystal electrons in high-harmonic
generation,” Nature **523**, 572–575
(2015). [CrossRef]

**12. **G. P. Zhang, M. S. Si, M. Murakami, Y. H. Baiand, and T. F. George, “Generating high-order optical
and spin harmonics from ferromagnetic monolayers,”
Nat. Commun. **9**, 3031
(2018). [CrossRef]

**13. **Y. Murakami, M. Eckstein, and P. Werner, “High-harmonic generation in
Mott insulators,” Phys. Rev. Lett. **121**, 057405
(2018). [CrossRef]

**14. **D. Bauer and K. K. Hansen, “High-harmonic generation in
solid with and without topological edge states,”
Phys. Rev. Lett. **120**,
177401 (2018). [CrossRef]

**15. **G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, “Theoretical analysis of
high-harmonic generations in solids,” Phys.
Rev. Lett. **113**, 073901
(2014). [CrossRef]

**16. **H. Haug and S. W. Koch, *Quantum Theory of the Optical and
Electronic Properties of Semiconductors*
(World Scientific,
2004).

**17. **I. Floss, C. Lemell, G. Wachter, M. Pickem, J. Burgdörfer, X.-M. Tong, K. Yabana, and S. A. Sato, “Ab initio simulation of
high-order harmonic generation in solids,”
Phys. Rev. A **97**,
011401 (2018). [CrossRef]

**18. **N. Tancogne-Dejean, O. D. Mücke, F. X. Kärtner, and A. Rubio, “Impact of electronic band
structure in high-harmonic generation spectra of
solids,” Phys. Rev. Lett. **118**, 087403
(2017). [CrossRef]

**19. **T. Otobe, “High-harmonic generation in
α-quartz by electron-hole recombination,”
Phys. Rev. B **94**,
235512 (2016). [CrossRef]

**20. **K. K. Hansen, T. Deffge, and D. Bauer, “High-order harmonic generation
in solid slabs beyond the single-active-electron
approximation,” Phys. Rev. A **96**, 053418 (2017). [CrossRef]

**21. **T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto, M. Kuwata-Gonokami, and K. L. Ishikawa, “Time-dependent Hartree-Fock
study of electron-hole interactions effects on high-order harmonic
generation from periodic crystals,” Phys. Rev.
A **98**, 023415
(2018). [CrossRef]

**22. **M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B. Gaarde, “High-harmonic generation from
Bloch electrons in solids,” Phys. Rev.
A **91**, 043839
(2015). [CrossRef]

**23. **T. Higuchi, M. I. Stockman, and P. Hommelhoff, “Strong-field perspective on
high-harmonic radiation from bulk solids,”
Phys. Rev. Lett. **113**,
213901 (2014). [CrossRef]

**24. **N. Moiseyev, “Selection rules for harmonic
generation in solids,” Phys. Rev. A **91**, 053811 (2015). [CrossRef]

**25. **P. G. Hawkins, M. Y. Ivanov, and V. S. Jakovlev, “ Effect of multiple conduction
bands on high-harmonic emission from dielectrics,”
Phys. Rev. A **91**,
013405 (2015). [CrossRef]

**26. **T.-Y. Du and X.-B. Bian, “Quasi-classical analysis of
the dynamics of the high-order harmonic generation from
solids,” Opt. Express **25**, 151–158
(2017). [CrossRef]

**27. **J.-B. Li, X. Zhang, S.-J. Yue, H.-M. Wu, B.-T. Hu, and H.-C. Du, “Enhancement of the second
plateau in solid high-order harmonic spectra by the two-color
fields,” Opt. Express **25**, 18603–18613
(2017). [CrossRef]

**28. **X. Liu, X. Zhu, X. Zhang, D. Wang, P. Lan, and P. Lu, “Wavelength scaling of the
cutoff energy in the solid high harmonic generation,”
Opt. Express **25**,
29216–29224 (2017). [CrossRef]

**29. **S. Jiang, J. Chen, H. Wei, C. Yu, R. Lu, and C. D. Lin, “Role of the transition dipole
amplitude and phase on the generation of odd and even high-order
harmonics in crystals,” Phys. Rev.
Lett. **120**, 253201
(2018). [CrossRef]

**30. **Y. W. Kim, T.-J. Shao, H. Kim, S. Han, S. Kim, M. Ciappina, X.-B. Bian, and S.-W. Kim, “Spectral interference in high
harmonic generation from solids,” ACS
Photon. **6**,
851–857 (2019). [CrossRef]

**31. **N. Tsatrafyllis, S. Kühn, M. Dumergue, P. Foldi, S. Kahaly, E. Cormier, I. A. Gonoskov, B. Kiss, K. Varju, S. Varro, and P. Tzallas, “Quantum optical signatures in
a strong laser pulse after interaction with
semiconductors,” Phys. Rev. Lett. **122**, 193602
(2019). [CrossRef]

**32. **G. Orlando, C.-M. Wang, T.-S. Ho, and S.-I. Chu, “High-order harmonic generation
in disordered semiconductors,” J. Opt. Soc.
Am. B **35**,
680–688 (2018). [CrossRef]

**33. **G. Orlando, T.-S. Ho, and S.-I. Chu, “Macroscopic effects in
high-order harmonic generation in disordered
semiconductors,” J. Opt. Soc. Am. B **36**, 1873–1880
(2019). [CrossRef]

**34. **S. Almalki, A. M. Parks, G. Bart, P. B. Corkum, T. Brabec, and C. R. McDonald, “High harmonic generation
tomography of impurities in solids: conceptual
analysis,” Phys. Rev. B **98**, 144307
(2018). [CrossRef]

**35. **T. Huang, X. Zhu, L. Li, X. Liu, P. Lan, and P. Lu, “High-order harmonic generation
of a doped semiconductor,” Phys. Rev.
A **96**, 043425
(2017). [CrossRef]

**36. **C. Yu, K. K. Hansen, and L. B. Madsen, “Enhanced high-order harmonic
generation in donor-doped band-gap materials,”
Phys. Rev. A **99**,
013435 (2019). [CrossRef]

**37. **C. Yu, K. K. Hansen, and L. B. Madsen, “High-order harmonic generation
in imperfect crystals,” Phys. Rev. A **99**, 063408 (2019). [CrossRef]

**38. **J. Ma, C. Zhang, H. Cui, Z. Ma, and X. Miao, “Theoretical investigation of
the electron dynamics in high-order harmonic generation process from
the doped periodic potential,” Chem. Phys.
Lett. **744**, 137207
(2020). [CrossRef]

**39. **K. Chinzei and T. N. Ikeda, “Disorder effects on the origin
of high-order harmonic generation in solids,”
Phys. Rev. Res. **2**,
013033 (2020). [CrossRef]

**40. **A. Pattanayak, M. S. Mrudul, and G. Dixit, “Influence of vacancy defects
in solid high-order harmonic generation,”
Phys. Rev. A **101**,
013404 (2020). [CrossRef]

**41. **M. S. Mrudul, N. Tancogne-Dejean, A. Rubio, and G. Dixit, “High-harmonic generation by
spin-polarized defects in solids,” Comput.
Mater. **6**, 10
(2020). [CrossRef]

**42. **D. C. Khandekar and S. V. Lawande, “Feynman path integrals: some
exact results and applications,” Phys.
Rep. **137**,
115–229 (1986). [CrossRef]

**43. **V. Bezak, “Path integral theory of an
electron gas in a random potential,” Proc. R.
Soc. London A **315**,
339–354 (1970). [CrossRef]

**44. **A. V. Chaplik, “The Green’s function and
mobility of an electron in a random potential,”
Sov. Phys. JETP **26**,
797–800 (1968).

**45. **S. F. Edwards and V. B. Gulyaev, “The density of states of an
highly impure semiconductor,” Proc. Phys.
Soc. **83**,
495–496
(1964). [CrossRef]

**46. **T. Lukes, “On the electronic structure of
disordered systems,” Philos. Mag. **12**(118),
719–724 (1965). [CrossRef]

**47. **E. Akkermans and G. Montambaux, *Mesoscopic Physics of Electrons and
Phonons* (Cambridge
University, 2007).

**48. **J. M. Ziman, *Model of Disorder*
(Cambridge University,
1979).

**49. **E. I. Blount, “Formalisms of band
theory,” Solid State Phys. **13**, 305–373
(1962). [CrossRef]

**50. **L. V. Keldysh, “Ionization in the field of a
strong electromagnetic wave,” Sov. Phys.
JETP **20**,
1307–1314 (1965).

**51. **D. Bennhardt, P. Thomas, A. Weller, M. Lindberg, and S. W. Koch, “Influence of Coulomb
interactions on the photon echo in disordered
semiconductors,” Phys. Rev. B **43**, 8934–8945
(1991). [CrossRef]

**52. **R. P. Feynman and A. R. Hibbs, *Quantum Mechanics and Path
Integrals* (McGraw-Hill,
1965).

**53. **V. Samathiakanit, “An average propagator of a
disordered system,” J. Phys. A **6**, 632–639
(1973). [CrossRef]

**54. **V. Sayakanit, “Mobility of an electron in a
Gaussian random potential,” Phys.
Lett. **59**,
461–463 (1977). [CrossRef]

**55. **D. C. Khandekar, V. A. Singh, K. V. Bhagwat, and S. V. Lawande, “Functional integral approach
to positionally disordered systems,” Phys.
Rev. B **33**,
5482–5486 (1986). [CrossRef]

**56. **I. S. Gradshteyn and I. M. Ryzhik, *Tables of Integrals, Series, and
Products* (Academic,
2015).

**57. **J. Y. Yan, “Theory of excitonic high-order
sideband generation in semiconductors under a strong terahertz
field,” Phys. Rev. B **78**, 075204 (2008). [CrossRef]

**58. **G. Vampa, C. R. McDonald, G. Orlando, P. B. Corkum, and T. Brabec, “Semiclassical analysis of high
harmonic generation in bulk crystals,” Phys.
Rev. B **91**, 064302
(2015). [CrossRef]

**59. **G. Vampa and T. Brabec, “Merge of high harmonic
generation from gases and solids and its implications for attosecond
science,” J. Phys. B **50**, 083001 (2017). [CrossRef]

**60. **N. B. Delone and V. P. Krainov, “Energy and angular electron
spectra for the tunnel ionization of atoms by strong low-frequency
radiation,” J. Opt. Soc. Am. B **8**, 1207–1211
(1991). [CrossRef]

**61. **P. B. Corkum, “Plasma perspective on strong
field multiphoton ionization,” Phys. Rev.
Lett. **71**,
1994–1997 (1993). [CrossRef]

**62. **M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, “Theory of high-harmonic
generation by low-frequency laser fields,”
Phys. Rev. A **49**,
2117–2132 (1994). [CrossRef]

**63. **L. D. Landau and E. M. Lifshitz, *Mechanics*
(Pergamon,
1960).

**64. **M. Protopapas, D. G. Lappas, C. H. Keitel, and P. Knight, “Recollisions, bremsstrahlung,
and attosecond pulses from intense laser fields,”
Phys. Rev. A **53**,
R2933–R2936 (1995). [CrossRef]

**65. **V. Turkowski and C. A. Ullrich, “Time-dependent
density-functional theory for ultrafast interband
excitations,” Phys. Rev. B **77**, 075204 (2008). [CrossRef]

**66. **S. Y. Kruchinin, “Non-Markovian pure dephasing
in a dielectric excited by a few-cycle laser pulse,”
Phys. Rev. A **100**,
043839 (2019). [CrossRef]