Disorder Effects on the Origin of High-Order Harmonic Generation in Solids

We consider noninteracting electrons coupled to laser fields, and study perturbatively the effects of the lattice potential involving disorder on the harmonic components of the electric current, which are sources of high-order harmonic generation (HHG). By using the Floquet-Keldysh Green functions, we show that each harmonic component consists of the coherent and the incoherent parts, which arise respectively from the coherent and the incoherent scatterings by the local ion potentials. As the disorder increases, the coherent part decreases, the incoherent one increases, and the total harmonic component of the current first decreases rapidly and then approaches a nonzero value. Our results highlight the importance of the periodicity of crystals, which builds up the Bloch states extending over the solid. This is markedly different from the traditional HHG in atomic gases, where the positions of individual atoms are irrelevant.


I. INTRODUCTION
High-order harmonic generation (HHG), which was traditionally studied in atomic and molecular gases [1], has recently been extended to solids [2][3][4][5][6][7][8][9][10][11]. Owing to the state-of-the-art optics technology, HHG is important not only as the foundation of the attosecond physics [12,13] but also as a probe to study electron dynamics in solids under strong laser electric fields. HHG from solids, unlike those from gases, should reflect band structures and crystal symmetries [14], and active studies are ongoing to elucidate their principles and applications.
Mechanisms and characteristics of HHG in solids are often theoretically studied by effective two-band models or, equivalently, tight-binding models. For semiconductors, two-band models [15][16][17][18][19][20][21][22] revealed that both the interband transitions and the intraband dynamics are important sources of HHG although their relative importances seem intricate. The intraband contribution is enhanced when the band dispersion has more anharmonicity [2], whereas the interband contribution remains significant even without intraband dynamics in the valence band [23,24]. The mechanism of HHG in semiconductors is still under active debate. The effective two-band models are also useful to study HHG in, for example, graphene [25,26], Mott insulators [27,28], charge-density-wave materials [29,30], superconductors [31,32], and topological insulators [33,34], where both the inter-and the intra-band contributions play important roles. We emphasize that both contributions originate from the periodic lattice potential, which is presupposed and expressed by the band dispersions and the interband coupling in those effective models.
To discuss the very origin of HHG in solids, different kinds of models are used, where the periodic lattice potential appears explicitly in the Hamiltonian [35][36][37][38][39][40][41][42][43][44][45][46][47]. In these models, it is manifest that there is no harmonic generation in the absence of the periodic lattice potential no matter how strong the laser field becomes (see e.g. Ref. [46]). Once the lattice potential is introduced, the quadratic energy dispersion is folded to form the energy bands, and the laser field causes both the anharmonic intraband dynamics and the interband transitions [41][42][43]. HHG has been studied in these models and the results seem consistent with those of the effective models. However, it has not been well studied how HHG changes when the lattice potential is not perfectly periodic. Given that the periodic lattice potential is the origin of HHG in solids, how important is the perfect periodicity of the lattice potential?
In this paper, we address this question and investigate the origin of HHG in solids by considering the highharmonic current (HHC), which gives rise to HHG, under strong laser fields on a weak lattice potential involving disorder. By invoking the Floquet-Keldysh Green functions, we show that the HHC consists of the coherent and incoherent parts, and rapidly decreases to a nonzero value as the disorder increases. On the basis of these results, we highlight the difference between the mechanisms of the HHG in solids and gases: the positions of atoms are relevant in solids while irrelevant in gases. These results show the importance of the periodicity of the lattice potential in HHG from solids, and reinforce the fact that its origin is the coherent dynamics of the Bloch state. This paper is organized as follows. In Sec. II, we formulate the problem of calculating HHG from disordered solids by using the Floquet-Keldysh formalism. In Sec. III, we analyze the disorder effects on HHG by both analytical and numerical approaches, uncovering the coherent and incoherent nature of HHG in solids. In Sec. IV, we discuss the essential difference between HHG in solids and gases by using our results. Finally in Sec. V, we summarize our study with concluding remarks.

II. FORMULATION OF THE PROBLEM
In this section, we formulate the problem of calculating HHC driven by a strong ac electric field in the presence of the disordered lattice potential. We begin by solving the problem in the absence of the potential, and show that arXiv:1905.05205v3 [cond-mat.dis-nn] 10 Dec 2019 there is no HHC. Then we introduce our model of the lattice potential [see Eq. (6)], and derive the formula for the HHC in terms of the Floquet-Keldysh Green function.

A. Electron dynamics without potential
We begin by analyzing the electron dynamics in the absence of the lattice potential. Let us consider noninteracting electrons coupled to a homogeneous ac electric field at frequency Ω in d dimensions. We ignore the spin degree of freedom since it merely doubles our results. We represent the ac field by the vector potential A(t) = A 0 cos Ωt in the velocity gauge. The Hamiltonian is given bŷ where m and e (< 0) are the mass and the electric charge of the electron, respectively, andp is the momentum operator.
The solutions of the time-dependent Schrödinger equation, i∂ t ψ(t) =Ĥ 0 (t)ψ(t), are known as the Volkov states [48][49][50] ( = 1 throughout this paper). Their wave functions are characterized by the momentum k and given by where k = k 2 /2m, α k = eA 0 · k/(mΩ), and J n (z) denotes the n-th Bessel function of the first kind. The Volkov states carry no HHC no matter how strong the ac electric field is. In fact, their paramagnetic and diamagnetic currents are given as j para (k, t) = −(ie/m) drψ k (r, t) * ∇ψ k (r, t) = ek/m and j dia (k, t) = −(e 2 /m) drψ k (r, t) * A(t)ψ k (r, t) = (e 2 /m)A(t) in appropriate wave-function normalizations. These results explicitly show that the Volkov state carries no HHC, or the Fourier components at nΩ with |n| ≥ 2.
We note that j dia (k, t) = (e 2 /m)A(t) is satisfied not only by the Volkov state but also by any states. Thus we ignore this part and focus on j para (k, t) in the following. Unlike j dia (k, t), j para (k, t) may involve harmonics when a lattice potential exists and k is no longer a good quantum number. Our problem is thus to study the Fourier components of j para (k, t) in the presence of the potential. We remark that, when the lattice potential is periodic, the energy bands are formed and the paramagnetic current has both diagonal and off-diagonal matrix elements in the band basis. In other words, the paramagnetic current involves both contributions from the interband transition and the intraband dynamics of electrons.
For later use, we express the Volkov states in terms of the Floquet theory [51]. The Floquet Hamiltonian H F 0 (k) corresponding to Eq. (1) reads where m, n ∈ Z are the Floquet indices, and H 0 (k, t) is the Fourier component ofĤ 0 (t). Its eigenvectors |ψ M 0 (k) and eigenvalues M k are labeled by M ∈ Z and given by where |φ n (k) is the Floquet state corresponding to the wave function ∝ e ik·r e −inΩt . We note that |ψ M 0 (k) is the Floquet representation of the Volkov state (2).

B. Model of lattice potential
To investigate the effects of the lattice potential on the HHC, we introduce the following lattice potential, where U denotes the strength of the potential, u(r) is a dimensionless function localized at r ∼ 0, and r a denotes the position of the scattering center, i.e. the ion in solids [52]. Without loss of generality, we assume dr u(r) = 0 since a constant shift of the total energy changes no physical consequences. This assumption implies that the Fourier component V k = U a u k e −ik·ra vanishes at k = 0, and thus the diagonal matrix elements of the lattice potential vanish in the momentum basis: k|V |k = V k−k = 0, where |k is the momentum eigenstate with eigenvalue k.
Suppose that r a 's form an approximate Bravais lattice such as the simple square lattice (see Fig. 1). Let R a denote the position of each lattice point a of the Bravais lattice, which is characterized by a set of d integers {c i a ∈ Z | i = 1, 2, · · · , d} as R a = d i=1 c i a a i with the primitive vectors a i . Then we introduce a small deviation δr a and define r a as r a = R a + δr a . Unlike R a 's, r a 's do not have the exact discrete translational symmetry. We note that our model describes the amorphous solids where the ions do not array periodically rather than the doped semiconductors [53][54][55] where the impurity potentials are periodically added and the electrons are trapped to them.
We assume that each of δr a 's is an independent Gaussian random variable. Its probability density function is given by P (δr a ) = (2πσ 2 ) −1/2 exp −δr 2 a /(2σ 2 ) , where the standard deviation σ quantifies the randomness of the lattice. Our assumption of the independence of the variables means, for instance, δr a δr b = δr a δr b for a = b, where · · · denotes the average over the random variables. In the following, we analyze the HHC for a given set of {r a } a and take its average over the deviations {δr a } a .
Here we show the band structure of our model. In the absence of the lattice potential V (r), the energy dispersion is parabolic as shown in Fig. 2, where the parabolic dispersion is folded into the first Brillouin zone for convenience. In the presence of V (r) without disorder (when there is disorder, the band picture is no longer valid since the discrete translational symmetry is broken), the energy gaps open at the degeneracy points of the folded dispersion, and the energy bands are formed as schematically shown in the inset of Fig. 2. We note that the energy gap is O(U ) when the lattice potential is a perturbation as we assume later. This means that our model is closer to narrow-gap semiconductors, metals, or Dirac materials rather than semiconductors.

C. Floquet-Keldysh formalism
We analyze the HHC by using the Floquet-Keldysh formalism, which is a method of combining the Floquet theory and the non-equilibrium Green function (see e.g., Refs. [56,57]). The non-equilibrium Green function has three components, the retarded, the advanced and the Keldysh Green functions. The retarded and the advanced Green functions have information of the energy spectrum, whereas the Keldysh Green function has that of the energy spectrum and the distribution function.
Without the lattice potential, the analytical solutions of the Green functions denoted by g are available since all the eigenstates of the Floquet Hamiltonian H F 0 (k) are obtained. The retarded and advanced Green functions are given by [58,59] (see also Ref. [46]) where m, n ∈ Z are the Floquet indices. We note that the band degrees of freedom are included in k; for instance, |k| < π/a lat (π/a lat < |k| < 2π/a lat ) corresponds to the lowest band (the second-lowest band) in Fig. 2.
Here, we assume a finite relaxation time 1/η (> 0), which stems from the electron correlation and so on. Because of the finite relaxation time, the Keldysh component has a Lorentzian form, where n k is a distribution function. In Sec. III B of the numerical calculation part, we assume that only the lowest band is occupied by the electrons and the other bands are empty (see Fig. 2): n k = 1 for |k| = k ≤ k F and n k = 0 for k > k F , where k F = π/a lat is the Fermi momentum and a lat is the lattice constant.
Once we have the non-equilibrium Green functions, the Fourier components of the paramagnetic current j para (k, t) at frequency nΩ is obtained as where v k = ek/m is the electron velocity, and g < = (g A − g R + g K )/2 is the lesser Green function [See Appendix A for the derivation of Eq. (9)]. From Eq. (7)-(9), we have j n (k) = δ n0 n k ek/m, which is consistent with j para (k, t) = ek/m obtained above directly by the wave function (2). We note that, in the absence of the lattice potential, the net paramagnetic current j n = d d k/(2π) d j n (k) vanishes because j n (k) = −j n (−k) if n k = n −k (Note that the diamagnetic current j dia exists but involves only the fundamental frequency Ω as noted above).

III. ANALYSIS OF HIGH-ORDER HARMONIC CURRENTS
In this section, we analyze the HHC in the presence of the disordered lattice potential on the basis of the The solid line with (n, k) indicates the n-photon state with momentum k, and the dot denotes the potential scattering. After random averaging, jn(k) consists of the two diagrams [see Eqs. (15)- (18)]. The first (second) diagram on the righthand side corresponds to the two scatterings occurring at different sites (the same site).
Floquet-Keldysh formalism. We perform the leadingorder calculation of the Green functions and derive the expressions for the HHC [see Eqs. (15)-(18)]. We then numerically evaluate those expressions and discuss the disorder effects on the HHC.

A. Analytical calculations of HHC
In the presence of the lattice potential (6), it is difficult to obtain the exact solutions to the Green functions, which are denoted by G. We invoke the perturbation theory in terms of the potential amplitude U to approximately obtain G, and calculate the leading-order correction of the HHC induced by the potential by Eq. (9) with g < replaced by G < .
The HHC does not appear up to the first order of U . In fact, at the zeroth order, G corresponds to g and there is no HHC as mentioned above. Besides, the O(U ) corrections to the Green functions induce no HHC. This is because V k vanishes at k = 0 and the corrections are made only on the off-diagonal components G kk (k = k ), while the HHC originates from the diagonal elements [see Eq. (9)].
The HHC first appears at the second order. According to the second-order Feynmann rules for the nonequilibrium Green function [60], the O(U 2 ) corrections to the full Green functions are given by where both G (2) mn (k, ω) and g mn (k, ω) are the 2 × 2 matrices, The technical details to derive and evaluate Eq. (10) are shown in Appendix B.
To obtain the results for the HHC, we take the average of the HHC over the random variables δr a 's, We note that the random variables appear only in the scattering vertices |V k−k | 2 in Eq. (13) in the following form: where B denotes the reciprocal vector. To obtain Eq. (14), we have used e ik·δra = e −σ 2 k 2 /2 and a e ik·Ra = B δ(k + B). The first and the second terms on the right-hand side of Eq. (14) correspond respectively to the terms with a = b and a = b on the left-hand side. In other words, these terms imply the double potential scattering at two different sites and the same site.
Corresponding to this decomposition, we decompose the HHC into two parts, where j C n (j IC n ) stems from the double scattering at two different sites (the same site), and has the coherent (incoherent) nature as discussed below [61]. Each current is given by where f n (k, k ) is These relations are the main results of the present work. We note j n = 0 for even n's when the system is inversion symmetric after random averaging. One can prove this by noticing that contributions from ±k cancel out each other.
We remark the similarity and the difference between j C n and j IC n . The similarity is that the n-th component derives from the Feynman diagram in which the initial and the final momenta are the same but the Floquet indices ("photon numbers") differ by n (see Fig. 3). The difference is that j IC n is contributed by any internal momentum whereas j C n by the discrete internal momenta shifted by reciprocal lattice vectors. This difference reflects the coherent (incoherent) nature of the potential scatterings in j C n (j IC n ). The physical meaning of the scattering becomes more obvious by introducing the band picture (see Fig. 2). The scattering process in j C n with the reciprocal lattice vector B can be interpreted as the interband transition without changing the lattice momentum, or k modulo B. On the other hand, the incoherent current j IC n is contributed by the scattering with any k and thus involves both the inter-and intra-band transitions.
When the disorder is very small, or σ 0, the n-th harmonic current j n is dominated by j C n . In the limit of σ → 0, j IC n vanishes and j C n reduces to the result for the perfectly periodic lattice. In this limit, j C n is the largest because the phase factor e i(k−k )·(ra−r b ) for the double scattering at any pair of sites a and b becomes unity. Namely, the scatterings at different sites are all coherent. As the disorder increases, j C n exponentially decays in σ 2 . This is because the vertex phase factors e i(k−k )·ra fluctuate and the scatterings at different pairs of sites work destructively.
On the other hand, j IC n becomes dominant when the disorder is very large. In the limit of σ → ∞, j C n vanishes and j IC n converges to a nonzero value. In this limit, the fluctuations δr a are so large that the phase factor e i(k−k )·(ra−r b ) is nonvanishing only for a = b. In other words, j IC n consists of the incoherent sum of the contributions from each scattering center. This local nature of j IC n manifests as the presence of any momentum k in Eq. (17). We note that j IC n vanishes at σ = 0, where the lattice potential is perfectly periodic. Thus j IC n is specific to disordered systems.

B. Numerical evaluation and quantitative analysis
Now we numerically evaluate Eqs. (16) and (17) for a choice of the potential u(r) and perform quantitative analyses. We work in one dimension for simplicity. We adopt u(x) = e −32x 2 /a 2 lat cos(16x/a lat ) + e −32x 2 /a 2 lat −2 , which is localized around x = 0 and satisfies dx u(x) = 0. Here a lat = 5Å is a typical lattice constant and we set the parameters as k F /a lat = π, Ω = 0.27 eV, and τ = 1/η = 48 fs [16]. Since Ω is above the band gap of O(U ), our setup is closer to narrow-gap semiconductors, metals, or Dirac materials rather than semiconductors. Thus, in applying our ideal model calculations directly to experimental situations, one should be careful because a very intense laser may cause the damage of sample solids, which is neglected in our calculation.
The crossover between the two limits σ = 0 and σ → ∞ is shown in Fig. 4(a), where |j C n | and |j IC n | evaluated by Eqs. (15)-(18) are plotted against σ for n = 5, 17, 29, and 35. As discussed above, |j C n | exponentially decays whereas |j IC n | firstly increases at σ/a lat ∼ 0 and approaches a nonzero value. The coherent part is larger than the incoherent one for (σ/a lat ) 2 0.1, and vice versa for (σ/a lat ) 2 0.1. This leads to the crossover of the HHC between the regimes where the coherent and the incoherent parts are dominant. The harmonic current | j n |, which is the sum of j C n and j IC n , is plotted against σ for n = 5, 17, 29, and 35 in Fig. 4(b). From the figure, we confirm that | j n | decreases exponentially at σ ∼ 0 (see also Ref. [62]), and approaches a nonzero value as σ increases. This is the crossover from the coherent scattering regime to the incoherent one. Also, we emphasize that the harmonic currents j n at σ = 0 is much larger than those at σ → ∞. This result shows the importance of the periodicity of the lattice potential to obtain large HHG.
We remark that the initial decrease of the HHC is more rapid for the higher order n. In fact, at small σ in Fig. 4(b), we observe | j n | ∝ e −(2πσ/a lat ) 2 for n = 5 and 17 whereas | j n | ∝ e −(4πσ/a lat ) 2 for n = 29 and 35. This implies that the low-(high-)order coherent current is dominantly contributed from the potential scattering with B = 2π/a lat (4π/a lat ) in Eq. (16). In terms of the band picture, the scattering with B = 2π/a lat (4π/a lat ) corresponds to the scattering from the lowest band to the second-lowest band (the third-lowest band).
The disorder has nontrivial effects on the structure of the harmonic spectrum through this difference. In Fig. 5(b) (E 0 = A 0 Ω = 10.6 MV/cm), we observe two plateaus in the absence of disorder σ = 0. As σ increases, the second plateau decreases more rapidly than the first one, and the only first plateau remains in the limit of σ → ∞, where j n = j IC n . This difference between the decay speeds of the two plateaus is consistent with the observation in Fig. 4(b), which shows that the higher-order harmonic currents decay more rapidly than the lower-order ones. Therefore, as discussed above, these imply that the second plateau is made up by the coherent current due to B = 4π/a lat whereas the first plateau is by both the coherent part due to B = 2π/a lat and the incoherent one.
We show the results for a weaker or a stronger laser field. Figure 5(a) shows the result for a weaker laser field (E 0 = 5.3 MV/cm), where the second plateau vanishes even for σ = 0. This implies that the scattering with large B is hard to occur for a weak laser because it needs the multiphoton absorption. For a stronger laser field (E 0 = 18.0 MV/cm), we obtain qualitatively similar result to that for E 0 = 10.6 MV/cm. As shown in Fig. 5(c), the first and the second plateaus merge to form a wide plateau at σ = 0. However, the second plateau decreases to vanish as σ → ∞, and only the first plateau remains in the limit. We note that the first plateau is wider for the stronger field E 0 .

IV. DISCUSSIONS
Although our setup is different from that for semiconductors as mentioned above, it is intriguing to point out similarities of our results with the HHG experiments on amorphous silica [63,64]. In this experiment [63], they have observed high harmonics with order n < 34 by the peak field of E 0 ∼ 200 MV/cm, while their efficiency is lower than those in crystalline quartz. Our results have two features observed in this experiment. First, we have shown that, as increases, the HHC decreases but approaches a nonzero value rather than vanishes. Second, our first plateau extends over n > 20 at E 0 > 18.0 MV/cm even in the presence of the disorder. These two characteristics are in line with the surprising observation that such a high order has been detected in amorphous solids.
Let us finally discuss the difference between the mechanisms of the HHG in gases and solids. The HHG in an atomic gas is explained by the celebrated three-step model [65]. According to this theory, the tunneling ionization occurs at an atom, an electron propagates in the laser field, and the electron and the ion recombine to produce radiation. These processes occur at individual atoms and their positions are irrelevant as shown in Fig. 6(a).
Our results highlight the essential difference of the HHG in solids from that in gases [66]. In solids, the HHG is induced by the potential scatterings by the ions, whose positions are of crucial importance. The HHG becomes the largest when they are aligned periodically and the scatterings at different positions are coherent. To put this in the momentum space, the electrons in the periodic potential are in the Bloch states whose wave functions extend over the entire crystal and are compatible with the periodicity [67]. When driven by a strong laser field, the Bloch states produce the HHC and hence the HHG as schematically depicted in Fig. 6(b). When the scattering centers fluctuate and the disorder sets in, the HHC rapidly decreases. As schematically illustrated in Fig. 6(c), the disorder disturbs the coherence of the Bloch states, and the resultant HHC, or the HHG, is suppressed.

V. CONCLUSIONS
We have studied the disorder effects on the HHG, or the HHC, in solids by considering the lattice potential (6), where the scattering centers fluctuate around a Bravais lattice. We have noted that, if there is no potential, the electrons are in the Volkov states and the HHC does not exist no matter how strong the laser field is. Then we have turned on a weak potential and analyzed the induced HHC by means of the perturbation theory. In other words, we have focused on the origin of the HHG in solids.
Our main results (15)- (18) have been obtained by the Floquet-Keldysh formalism, stating that the n-th harmonic current j n averaged over the fluctuations of the potential is the sum of the coherent part j C n and the incoherent one j IC n . Here j C n and j IC n are due to the double scattering at two different sites a = b and the same site a = b, respectively. In the absence of the disorder (σ = 0), where the lattice is perfectly periodic, j C n is the largest and j IC n vanishes. As σ increases, j C n decreases exponentially and j IC n increases to be dominant. Namely, the HHC exhibits a crossover from the coherent to the incoherent ones as the disorder increases. We have numerically shown that the total HHC j n is the largest in the limit of no disorder. Besides, we have shown that the disorder significantly deforms the structure of the harmonic spectrum. As shown in Fig. 5, the first plateau is robust against the disorder whereas the second plateau is strongly suppressed by the disorder.
Our results shed light on the difference of the mechanisms of the HHG in gases and solids. In atomic gases, the HHG occurs at individual atoms and their positions are not relevant. In solids, however, the positions of the scattering centers are crucially important, and the HHG is the largest when they form a periodic lattice. This implies in the momentum space that the periodic potential makes up the Bloch states with coherence, which produce the largest radiation, and the disorder breaks the coherence to suppress the HHG.
Finally, we comment on the applicability of our results to experiments. In this study, we have considered a simplified model of the disordered lattice potential and analyzed the situation in which the photon energy is larger than the gap. This situation is different from those of typical HHG experiments in semiconductors and rather close to those in (semi)metals and narrow-gap semiconductors. Thus one should be careful to compare the experimental results with our theoretical ones. To study real disordered materials, one may need a more detailed Hamiltonian and possibly the interaction between electrons [68]. We leave this as an important open issue.
We assume G < (k; t, t ) = G < (k; t + T, t + T ), meaning that the Green function is periodic with respect to the center of time (t + t )/2. This assumption implies that we consider a nonequilibrium steady state represented by a mixed state of the Floquet eigenstates. From this assumption, we haveG < mn (k; ω, ω ) = 2πδ(ω − ω )G < mn (k, ω), and then the Green function at the same time t = t is given by Replacing the dummy variables of summation (m, n) with (l, n) (l = m−n) and using the relation G < m+1,n+1 (k, ω) = G < mn (k, ω − Ω), we obtain As a result, the Fourier component of j(t) at frequency nΩ is obtained as This is the formula by which we calculate the HHC from the non-equilibrium Green function.
Let us concretely calculate the second-order correction of the Green function from Eq. (B1)-(B3). From Eq. (B3), the retarded, the advanced and the Keldysh components are where we have defined the following abbreviations: Therefore, the HHC at frequency nΩ is where we have used G < = (G A − G R + G K )/2. We perform the ω-integration in Eq. (B10) by calculating the residues. Since the integrand vanishes fast enough for |ω| → ∞, we replace the integral path with the closed semicircle in the upper or the lower halves of the ω-plane. First, the ω-integration of AAA and RRR terms vanish because their integrands have poles only in the upper or the lower halves of the ω-plane (see Figs. 8(a) and (d)). Next, we perform the ω-integrations of RRK, KAA and RKA terms. From Eqs. (B1) and (B2), we calculate where we have calculated the residues owing to the second equality (see Fig. 8(b)). In the same way, we obtain .
(B13) Figure 8 shows the poles of the each term in Eq. (B10). Summing up all these terms, we obtain the following expression f n (k, k ) = U 2 |u k−k | 2 (n k − n k ) nΩ − 4iη nΩ − 2iη .
where we have used a J a (x)J b−a (y) = J b (x + y) and V k = U a u k e −ik·ra . Taking the average over the random variables δr a 's [see Eq. (14) in the main text] and substituting Eq. (B14) and (B15) to Eq. (B10), we obtain Eqs. (15)- (18) in the main text.
Appendix C: Relaxation-time insensitivity of the results Figure 9 shows the harmonic spectra for two relaxation times τ = 48 fs and 24 fs. The panel (a) is the same as Fig. 3(a) in the main text, and the panel (b) is calculated with τ changed and all the other parameters fixed. Evidently, the harmonic spectra are almost the same. Thus we confirm that the physical consequences are not sensitive to the relaxation time.