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

TL;DR

This paper investigates how disorder in solid lattices affects high-order harmonic generation, revealing the interplay between coherent and incoherent scattering and emphasizing the importance of crystal periodicity.

Contribution

It introduces a perturbative analysis of disorder effects on HHG in solids using Floquet-Keldysh Green functions, highlighting the role of lattice periodicity.

Findings

Coherent harmonic components decrease with disorder

Incoherent components increase as disorder grows

Total harmonic signal approaches a nonzero value despite disorder

Abstract

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.

Figures

Equations

$${\hat{H}}_0(t)=\frac{{\hat{\mathbf{p}}}^2}{2m}-\frac{e}{m}\hat{\mathbf{p}}\cdot \mathbf{A}(t),$$

$${\psi }_{\mathbf{k}}(\mathbf{r},t)$$

$${\left[H_0^{\text{F}}(\mathbf{k})\right]}_{mn}$$

$$=({ϵ}_{\mathbf{k}}+n\Omega ){\delta }_{mn}+{\alpha }_{\mathbf{k}}({\delta }_{m,n+1}+{\delta }_{m,n-1}),$$

$$|{\psi }_0^M(\mathbf{k})⟩$$

$${ϵ}_{\mathbf{k}}^M$$

$$V(\mathbf{r})$$

$$g_{mn}^{R/A}(\mathbf{k},\omega )$$

$$g_{mn}^K(\mathbf{k},\omega )$$

$${\mathbf{j}}_n(\mathbf{k})=-i{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{v}_{\mathbf{k}}{g}_{n0}^{<}(\mathbf{k},\omega ),$$

$$G_{mn}^{(2)}(\mathbf{k},\omega )$$

$$=\sum _{p,q}\int \frac{d^d{k}^{\prime }}{(2\pi {)}^d}{g}_{mp}(\mathbf{k},\omega )V_{\mathbf{k}-{\mathbf{k}}^{\prime }}{g}_{pq}({\mathbf{k}}^{\prime },\omega )V_{{\mathbf{k}}^{\prime }-\mathbf{k}}{g}_{qn}(\mathbf{k},\omega ),$$

$$G_{mn}^{(2)}(\mathbf{k},\omega )$$

$$g_{mn}(\mathbf{k},\omega )$$

$$⟨{\mathbf{j}}_n⟩=-i\int \frac{d^{d}k}{(2\pi {)}^d}\int \frac{d\omega }{2\pi }{v}_{\mathbf{k}}⟨G_{n0}^{<(2)}(\mathbf{k},\omega )⟩.$$

$$\sum _{a,b}⟨e^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{r}}_a}{e}^{i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{r}}_b}⟩$$

$$=\sum _{\mathbf{B}}\delta (\mathbf{k}-{\mathbf{k}}^{\prime }+\mathbf{B})e^{-{\sigma }^2|\mathbf{B}{|}^2}+(1-e^{-{\sigma }^2|\mathbf{k}-{\mathbf{k}}^{\prime }{|}^2}),$$

$$⟨{\mathbf{j}}_n⟩={\mathbf{j}}_n^{\text{C}}+{\mathbf{j}}_n^{\text{IC}},$$

$${\mathbf{j}}_n^{\text{C}}$$

$${\mathbf{j}}_n^{\text{IC}}$$

$$f_n(\mathbf{k},{\mathbf{k}}^{\prime })$$

$$\mathbf{j}(t)$$

$$=-\frac{ie}{2m}\int d^{d}r{⟨{\hat{\psi }}^{†}(\mathbf{r},t)[\nabla \hat{\psi }(\mathbf{r},t)]-[\nabla {\hat{\psi }}^{†}(\mathbf{r},t)]\hat{\psi }(\mathbf{r},t)⟩}_0,$$

$$\mathbf{j}(t)=\int \frac{d^{d}k}{(2\pi {)}^d}{v}_{\mathbf{k}}{⟨{\hat{\psi }}_{\mathbf{k}}^{†}(t){\hat{\psi }}_{\mathbf{k}}(t)⟩}_0,$$

$$G^{<}(\mathbf{k};t,t^{\prime })$$

$$={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{\int }_{-\infty }^{\infty }\frac{d{\omega }^{\prime }}{2\pi }{\tilde{G}}^{<}(\mathbf{k};\omega ,{\omega }^{\prime })e^{i\omega t}{e}^{-i{\omega }^{\prime }{t}^{\prime }}$$

$$=\sum _{m,n}{\int }_0^{\Omega }\frac{d\omega }{2\pi }{\int }_0^{\Omega }\frac{d{\omega }^{\prime }}{2\pi }{\tilde{G}}_{mn}^{<}(\mathbf{k};\omega ,{\omega }^{\prime })e^{i(\omega -m\Omega )t}{e}^{-i({\omega }^{\prime }-n\Omega )t^{\prime }},$$

$$G^{<}(\mathbf{k};t,t)$$

$${\mathbf{j}}_n$$

$$=-i\int \frac{d^{d}k}{(2\pi {)}^d}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{v}_{\mathbf{k}}{G}_{n,0}^{<}(\mathbf{k},\omega ).$$

$$g_{mn}^{R/A}(\mathbf{k},\omega )$$

$$=\sum _M\frac{J_{M-m}({\alpha }_{\mathbf{k}})J_{M-n}({\alpha }_{\mathbf{k}})}{\omega -{ϵ}_{\mathbf{k}}^M\pm i\eta },$$

$$g_{mn}^K(\mathbf{k},\omega )$$

$$G_{mn}^{(2)}(\mathbf{k},\omega )$$

$$=\sum _{p,q}\int \frac{d^d{k}^{\prime }}{(2\pi {)}^d}{g}_{mp}(\mathbf{k},\omega )V_{\mathbf{k}-{\mathbf{k}}^{\prime }}{g}_{pq}({\mathbf{k}}^{\prime },\omega )V_{{\mathbf{k}}^{\prime }-\mathbf{k}}{g}_{qn}(\mathbf{k},\omega ),$$

$$G_{mn}^{(2)}(\mathbf{k},\omega )$$

Full text

I Introduction

High-order harmonic generation (HHG), which was traditionally studied in atomic and molecular gases Brabec and Krausz (2000), has recently been extended to solids Ghimire et al. (2011); Schubert et al. (2014); Hohenleutner et al. (2015); Luu et al. (2015); Ndabashimiye et al. (2016); Vampa et al. (2015a); You et al. (2017a); Yoshikawa et al. (2017); Higuchi et al. (2017); Kaneshima et al. (2018). Owing to the state-of-the-art optics technology, HHG is important not only as the foundation of the attosecond physics Ghimire et al. (2014); Calegari et al. (2016) 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 Vampa et al. (2015b), 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 Golde et al. (2006, 2008, 2011); Vampa et al. (2014, 2015c); Higuchi et al. (2014); Tamaya et al. (2016a, b) 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 Ghimire et al. (2011), whereas the interband contribution remains significant even without intraband dynamics in the valence band Osika et al. (2017); Catoire et al. (2018). 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 Mikhailov and Ziegler (2008); Mikhailov (2016), Mott insulators Murakami et al. (2018); Murakami and Werner (2018), charge-density-wave materials Nag et al. (2018); Ikeda et al. (2018), superconductors Kawakami et al. (2018); Yonemitsu (2018), and topological insulators Bauer and Hansen (2018); Jürß and Bauer (2019), 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 Plaja and Roso-Franco (1992); Faisal and Kamiński (1997); Faisal et al. (2005); Yan (2008); Korbman et al. (2013); Park (2014); Wu et al. (2015); Du and Bian (2017); Ikemachi et al. (2017); Jia et al. (2017); Ikemachi et al. (2018); Ikeda (2018); Navarrete et al. (2019). 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. Ikeda (2018)). 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 Wu et al. (2015); Du and Bian (2017); Ikemachi et al. (2017). 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 high-harmonic 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 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.

II.1 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 $\Omega$ 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 $\bm{A}(t)=\bm{A}_{0}\cos\Omega t$ in the velocity gauge. The Hamiltonian is given by

[TABLE]

where $m$ and $e$ $(<0)$ are the mass and the electric charge of the electron, respectively, and $\hat{\bm{p}}$ is the momentum operator.

The solutions of the time-dependent Schrödinger equation, $i\partial_{t}\psi(t)=\hat{H}_{0}(t)\psi(t)$, are known as the Volkov states Keldysh (1965); Faisal (1973); Reiss (1980) ($\hbar=1$ throughout this paper). Their wave functions are characterized by the momentum $\bm{k}$ and given by

[TABLE]

where $\epsilon_{\bm{k}}=\bm{k}^{2}/2m$, $\alpha_{\bm{k}}=e\bm{A}_{0}\cdot\bm{k}/(m\Omega)$, 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 $\bm{j}_{\text{para}}(\bm{k},t)=-(ie/m)\int d\bm{r}\psi_{\bm{k}}(\bm{r},t)^{\ast}\nabla\psi_{\bm{k}}(\bm{r},t)=e\bm{k}/m$ and $\bm{j}_{\text{dia}}(\bm{k},t)=-(e^{2}/m)\int d\bm{r}\psi_{\bm{k}}(\bm{r},t)^{\ast}\bm{A}(t)\psi_{\bm{k}}(\bm{r},t)=(e^{2}/m)\bm{A}(t)$ in appropriate wave-function normalizations. These results explicitly show that the Volkov state carries no HHC, or the Fourier components at $n\Omega$ with $|n|\geq 2$.

We note that $\bm{j}_{\text{dia}}(\bm{k},t)=(e^{2}/m)\bm{A}(t)$ is satisfied not only by the Volkov state but also by any states. Thus we ignore this part and focus on $\bm{j}_{\text{para}}(\bm{k},t)$ in the following. Unlike $\bm{j}_{\text{dia}}(\bm{k},t)$, $\bm{j}_{\text{para}}(\bm{k},t)$ may involve harmonics when a lattice potential exists and $\bm{k}$ is no longer a good quantum number. Our problem is thus to study the Fourier components of $\bm{j}_{\text{para}}(\bm{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 Shirley (1965). The Floquet Hamiltonian $H^{\text{F}}_{0}(\bm{k})$ corresponding to Eq. (1) reads

[TABLE]

where $m,n\in\mathbb{Z}$ are the Floquet indices, and $H_{0}(\bm{k},t)$ is the Fourier component of $\hat{H}_{0}(t)$. Its eigenvectors $\ket{\psi^{M}_{0}(\bm{k})}$ and eigenvalues $\epsilon^{M}_{\bm{k}}$ are labeled by $M\in\mathbb{Z}$ and given by

[TABLE]

where $\ket{\phi^{n}(\bm{k})}$ is the Floquet state corresponding to the wave function $\propto e^{i\bm{k}\cdot\bm{r}}e^{-in\Omega t}$. We note that $\ket{\psi^{M}_{0}(\bm{k})}$ is the Floquet representation of the Volkov state (2).

II.2 Model of lattice potential

To investigate the effects of the lattice potential on the HHC, we introduce the following lattice potential,

[TABLE]

where $U$ denotes the strength of the potential, $u(\bm{r})$ is a dimensionless function localized at $\bm{r}\sim 0$, and $\bm{r}_{a}$ denotes the position of the scattering center, i.e. the ion in solids 111 We neglect the changes of the potential amplitude $U$ and the local ion potential $u(\bm{r})$ by the disorder. This simplification has also been made in theoretical studies of the solid-state HHG (see e.g. Ref. Yu et al. (2019a)) . Without loss of generality, we assume $\int d\bm{r}\,u(\bm{r})=0$ since a constant shift of the total energy changes no physical consequences. This assumption implies that the Fourier component $V_{\bm{k}}=U\sum_{a}u_{\bm{k}}e^{-i\bm{k}\cdot\bm{r}_{a}}$ vanishes at $\bm{k}=0$, and thus the diagonal matrix elements of the lattice potential vanish in the momentum basis: $\braket{\bm{k}}{\hat{V}}{\bm{k}}=V_{\bm{k}-\bm{k}}=0$, where $\ket{\bm{k}}$ is the momentum eigenstate with eigenvalue $\bm{k}$.

Suppose that $\bm{r}_{a}$’s form an approximate Bravais lattice such as the simple square lattice (see Fig. 1). Let $\bm{R}_{a}$ denote the position of each lattice point $a$ of the Bravais lattice, which is characterized by a set of $d$ integers $\{c_{a}^{i}\in\mathbb{Z}\mid i=1,2,\cdots,d\}$ as $\bm{R}_{a}=\sum_{i=1}^{d}c_{a}^{i}\bm{a}_{i}$ with the primitive vectors $\bm{a}_{i}$. Then we introduce a small deviation $\delta\bm{r}_{a}$ and define $\bm{r}_{a}$ as $\bm{r}_{a}=\bm{R}_{a}+\delta\bm{r}_{a}$. Unlike $\bm{R}_{a}$’s, $\bm{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 Huang et al. (2017); Almalki et al. (2018); Yu et al. (2019a) where the impurity potentials are periodically added and the electrons are trapped to them.

We assume that each of $\delta\bm{r}_{a}$’s is an independent Gaussian random variable. Its probability density function is given by $P(\delta\bm{r}_{a})=(2\pi\sigma^{2})^{-1/2}\exp\left[-\delta r_{a}^{2}/(2\sigma^{2})\right]$, where the standard deviation $\sigma$ quantifies the randomness of the lattice. Our assumption of the independence of the variables means, for instance, $\langle\delta\bm{r}_{a}\delta\bm{r}_{b}\rangle=\langle\delta\bm{r}_{a}\rangle\langle\delta\bm{r}_{b}\rangle$ for $a\neq b$, where $\langle\cdots\rangle$ denotes the average over the random variables. In the following, we analyze the HHC for a given set of $\{\bm{r}_{a}\}_{a}$ and take its average over the deviations $\{\delta\bm{r}_{a}\}_{a}$.

Here we show the band structure of our model. In the absence of the lattice potential $V(\bm{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(\bm{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.

II.3 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. Oka and Aoki (2009); Morimoto and Nagaosa (2016)). 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^{\text{F}}_{0}(\bm{k})$ are obtained. The retarded and advanced Green functions are given by Faisal and Genieser (1989); Faisal (1991) (see also Ref. Ikeda (2018))

[TABLE]

where $m,n\in\mathbb{Z}$ are the Floquet indices. We note that the band degrees of freedom are included in $\bm{k}$; for instance, $|\bm{k}|<\pi/a_{\text{lat}}$ ($\pi/a_{\text{lat}}<|\bm{k}|<2\pi/a_{\text{lat}}$) corresponds to the lowest band (the second-lowest band) in Fig. 2. Here, we assume a finite relaxation time $1/\eta$ $(>0)$, which stems from the electron correlation and so on. Because of the finite relaxation time, the Keldysh component has a Lorentzian form,

[TABLE]

where $n_{\bm{k}}$ is a distribution function. In Sec. III.2 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_{\bm{k}}=1$ for $|\bm{k}|=k\leq k_{F}$ and $n_{\bm{k}}=0$ for $k>k_{F}$, where $k_{F}=\pi/a_{\mathrm{lat}}$ is the Fermi momentum and $a_{\mathrm{lat}}$ is the lattice constant.

Once we have the non-equilibrium Green functions, the Fourier components of the paramagnetic current $\bm{j}_{\text{para}}(\bm{k},t)$ at frequency $n\Omega$ is obtained as

[TABLE]

where $v_{\bm{k}}=e\bm{k}/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 $\bm{j}_{n}(\bm{k})=\delta_{n0}n_{\bm{k}}e\bm{k}/m$, which is consistent with $\bm{j}_{\text{para}}(\bm{k},t)=e\bm{k}/m$ obtained above directly by the wave function (2). We note that, in the absence of the lattice potential, the net paramagnetic current $\bm{j}_{n}=\int d^{d}k/(2\pi)^{d}\bm{j}_{n}(\bm{k})$ vanishes because $\bm{j}_{n}(\bm{k})=-\bm{j}_{n}(-\bm{k})$ if $n_{\bm{k}}=n_{-\bm{k}}$ (Note that the diamagnetic current $\bm{j}_{\text{dia}}$ exists but involves only the fundamental frequency $\Omega$ 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 Floquet-Keldysh formalism. We perform the leading-order 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.

III.1 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_{\bm{k}}$ vanishes at $\bm{k}=0$ and the corrections are made only on the off-diagonal components $G_{\bm{k}\bm{k}^{\prime}}$ ($\bm{k}\neq\bm{k}^{\prime}$), 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 non-equilibrium Green function Rammer and Smith (1986), the $O(U^{2})$ corrections to the full Green functions are given by

[TABLE]

where both $G^{(2)}_{mn}(k,\omega)$ and $g_{mn}(k,\omega)$ are the $2\times 2$ matrices,

[TABLE]

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 $\delta\bm{r}_{a}$’s,

[TABLE]

We note that the random variables appear only in the scattering vertices $|V_{\bm{k}-\bm{k}^{\prime}}|^{2}$ in Eq. (13) in the following form:

[TABLE]

where $\bm{B}$ denotes the reciprocal vector. To obtain Eq. (14), we have used $\braket{e^{i\bm{k}\cdot\delta\bm{r}_{a}}}=e^{-\sigma^{2}\bm{k}^{2}/2}$ and $\sum_{a}e^{i\bm{k}\cdot\bm{R}_{a}}=\sum_{\bm{B}}\delta(\bm{k}+\bm{B})$. The first and the second terms on the right-hand side of Eq. (14) correspond respectively to the terms with $a\neq 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,

[TABLE]

where $\bm{j}_{n}^{\text{C}}$ ($\bm{j}_{n}^{\text{IC}}$) stems from the double scattering at two different sites (the same site), and has the coherent (incoherent) nature as discussed below 222 The coherent (incoherent) current $\bm{j}_{n}^{\text{C}}$ ($\bm{j}_{n}^{\text{IC}}$) is contributed by both the intra- and the inter-band dynamics in the three-step model in the momentum space Vampa et al. (2014). In the real space, the HHG is produced by the electrons recolliding with the various ion sites. As discussed in Ref. Ghimire and Reis (2019), these two interpretations describe the same phenomenon and thus are consistent with each other (see also Discussions section in the main text). . Each current is given by

[TABLE]

where $f_{n}(\bm{k},\bm{k}^{\prime})$ is

[TABLE]

These relations are the main results of the present work. We note $\langle\bm{j}_{n}\rangle=0$ for even $n$’s when the system is inversion symmetric after random averaging. One can prove this by noticing that contributions from $\pm\bm{k}$ cancel out each other.

We remark the similarity and the difference between $\bm{j}_{n}^{\text{C}}$ and $\bm{j}_{n}^{\text{IC}}$. 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 $\bm{j}_{n}^{\text{IC}}$ is contributed by any internal momentum whereas $\bm{j}_{n}^{\text{C}}$ by the discrete internal momenta shifted by reciprocal lattice vectors. This difference reflects the coherent (incoherent) nature of the potential scatterings in $\bm{j}_{n}^{\text{C}}$ ($\bm{j}_{n}^{\text{IC}}$). The physical meaning of the scattering becomes more obvious by introducing the band picture (see Fig. 2). The scattering process in $\bm{j}_{n}^{\text{C}}$ with the reciprocal lattice vector $\bm{B}$ can be interpreted as the interband transition without changing the lattice momentum, or $\bm{k}$ modulo $\bm{B}$. On the other hand, the incoherent current $\bm{j}_{n}^{\text{IC}}$ is contributed by the scattering with any $\bm{k}^{\prime}$ and thus involves both the inter- and intra-band transitions.

When the disorder is very small, or $\sigma\simeq 0$, the $n$-th harmonic current $\langle\bm{j}_{n}\rangle$ is dominated by $\bm{j}_{n}^{\text{C}}$. In the limit of $\sigma\to 0$, $\bm{j}_{n}^{\text{IC}}$ vanishes and $\bm{j}_{n}^{\text{C}}$ reduces to the result for the perfectly periodic lattice. In this limit, $\bm{j}_{n}^{\text{C}}$ is the largest because the phase factor $e^{i(\bm{k}-\bm{k}^{\prime})\cdot(\bm{r}_{a}-\bm{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, $\bm{j}_{n}^{\text{C}}$ exponentially decays in $\sigma^{2}$. This is because the vertex phase factors $e^{i(\bm{k}-\bm{k}^{\prime})\cdot\bm{r}_{a}}$ fluctuate and the scatterings at different pairs of sites work destructively.

On the other hand, $\bm{j}_{n}^{\text{IC}}$ becomes dominant when the disorder is very large. In the limit of $\sigma\rightarrow\infty$, $\bm{j}_{n}^{\text{C}}$ vanishes and $\bm{j}_{n}^{\text{IC}}$ converges to a nonzero value. In this limit, the fluctuations $\delta\bm{r}_{a}$ are so large that the phase factor $e^{i(\bm{k}-\bm{k}^{\prime})\cdot(\bm{r}_{a}-\bm{r}_{b})}$ is nonvanishing only for $a=b$. In other words, $\bm{j}_{n}^{\text{IC}}$ consists of the incoherent sum of the contributions from each scattering center. This local nature of $\bm{j}_{n}^{\text{IC}}$ manifests as the presence of any momentum $\bm{k}^{\prime}$ in Eq. (17). We note that $\bm{j}_{n}^{\text{IC}}$ vanishes at $\sigma=0$, where the lattice potential is perfectly periodic. Thus $\bm{j}_{n}^{\text{IC}}$ is specific to disordered systems.

III.2 Numerical evaluation and quantitative analysis

Now we numerically evaluate Eqs. (16) and (17) for a choice of the potential $u(\bm{r})$ and perform quantitative analyses. We work in one dimension for simplicity. We adopt $u(x)=e^{-32x^{2}/a_{\text{lat}}^{2}}\cos(16x/a_{\text{lat}})+e^{-32x^{2}/a_{\text{lat}}^{2}-2}$, which is localized around $x=0$ and satisfies $\int dx\,u(x)=0$. Here $a_{\text{lat}}=5\,\text{\AA}$ is a typical lattice constant and we set the parameters as $k_{F}/a_{\text{lat}}=\pi$, $\hbar\Omega=0.27$ eV, and $\tau=1/\eta=48$ fs Golde et al. (2008). Since $\hbar\Omega$ 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 $\sigma=0$ and $\sigma\to\infty$ is shown in Fig. 4(a), where $|\bm{j}_{n}^{\text{C}}|$ and $|\bm{j}_{n}^{\text{IC}}|$ evaluated by Eqs. (15)–(18) are plotted against $\sigma$ for $n=5,17,29$, and $35$. As discussed above, $|\bm{j}_{n}^{\text{C}}|$ exponentially decays whereas $|\bm{j}_{n}^{\text{IC}}|$ firstly increases at $\sigma/a_{\text{lat}}\sim 0$ and approaches a nonzero value. The coherent part is larger than the incoherent one for $(\sigma/a_{\text{lat}})^{2}\lesssim 0.1$, and vice versa for $(\sigma/a_{\text{lat}})^{2}\gtrsim 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 $|\langle\bm{j}_{n}\rangle|$, which is the sum of $\bm{j}_{n}^{\text{C}}$ and $\bm{j}_{n}^{\text{IC}}$, is plotted against $\sigma$ for $n=5,17,29$, and $35$ in Fig. 4(b). From the figure, we confirm that $|\langle\bm{j}_{n}\rangle|$ decreases exponentially at $\sigma\sim 0$ (see also Ref. Orlando et al. (2018)), and approaches a nonzero value as $\sigma$ increases. This is the crossover from the coherent scattering regime to the incoherent one. Also, we emphasize that the harmonic currents $\langle\bm{j}_{n}\rangle$ at $\sigma=0$ is much larger than those at $\sigma\rightarrow\infty$. 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 $\sigma$ in Fig. 4(b), we observe $|\langle\bm{j}_{n}\rangle|\propto\mathrm{e}^{-(2\pi\sigma/a_{\text{lat}})^{2}}$ for $n=5$ and 17 whereas $|\langle\bm{j}_{n}\rangle|\propto\mathrm{e}^{-(4\pi\sigma/a_{\text{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\pi/a_{\text{lat}}$ $(4\pi/a_{\text{lat}})$ in Eq. (16). In terms of the band picture, the scattering with $B=2\pi/a_{\text{lat}}$ ($4\pi/a_{\text{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}\Omega=10.6$ MV/cm), we observe two plateaus in the absence of disorder $\sigma=0$. As $\sigma$ increases, the second plateau decreases more rapidly than the first one, and the only first plateau remains in the limit of $\sigma\to\infty$, where $\bm{j}_{n}=\bm{j}_{n}^{\text{IC}}$. 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\pi/a_{\text{lat}}$ whereas the first plateau is by both the coherent part due to $B=2\pi/a_{\text{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 $\sigma=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 $\sigma=0$. However, the second plateau decreases to vanish as $\sigma\to\infty$, 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 You et al. (2017b); Jürgens et al. (2019). In this experiment You et al. (2017b), they have observed high harmonics with order $n<34$ by the peak field of $E_{0}\sim 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 B. Corkum (1993). 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 333 In the limit of $\sigma\rightarrow\infty$, our model corresponds to strongly disordered solids rather than gases. This is because the HHG in our model derives from multiscattering processes at several ion sites while that in gases does from individual atoms. . 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 Ghimire and Reis (2019). 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 $\langle\bm{j}_{n}\rangle$ averaged over the fluctuations of the potential is the sum of the coherent part $\bm{j}_{n}^{\text{C}}$ and the incoherent one $\bm{j}_{n}^{\text{IC}}$. Here $\bm{j}_{n}^{\text{C}}$ and $\bm{j}_{n}^{\text{IC}}$ are due to the double scattering at two different sites $a\neq b$ and the same site $a=b$, respectively. In the absence of the disorder ($\sigma=0$), where the lattice is perfectly periodic, $\bm{j}_{n}^{\text{C}}$ is the largest and $\bm{j}_{n}^{\text{IC}}$ vanishes. As $\sigma$ increases, $\bm{j}_{n}^{\text{C}}$ decreases exponentially and $\bm{j}_{n}^{\text{IC}}$ 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 $\langle\bm{j}_{n}\rangle$ 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 Yu et al. (2019b). We leave this as an important open issue.

acknowledgements

This work was supported by JSPS KAKENHI Grant No. JP18K13495. K. C. acknowledges financial support provided by the Advanced Leading Graduate Course for Photon Science at the University of Tokyo.

Appendix A Derivation of formula for HHC

We derive the formula for the HHC described by the non-equilibrium Green function (see Eq. (4) in the main text). Let $\hat{\psi}(\bm{r},t)$ denote the second-quantized field operator in the Heisenberg picture, and then the electric current is given by

[TABLE]

where $\braket{\cdots}_{0}$ denotes the expectation value in the quantum state, not random averaging. Substituting the Fourier expansion of $\hat{\psi}(\bm{r},t)$ and performing the integration of $\bm{r}$, we obtain

[TABLE]

where $\hat{\psi}_{\bm{k}}(t)$ denotes the Fourier component of $\hat{\psi}(\bm{r},t)$. Here we define the lesser Green function as $G^{<}(\bm{k};t,t^{\prime})\equiv i\left\langle\hat{\psi}^{\dagger}_{\bm{k}}(t^{\prime})\hat{\psi}_{\bm{k}}(t)\right\rangle_{0}$. The Fourier expansion of $G^{<}(\bm{k};t,t^{\prime})$ with respect to $t$ and $t^{\prime}$ is

[TABLE]

where we have used $\int_{-\infty}^{\infty}d\omega e^{i\omega t}=\sum_{m}\int_{0}^{\Omega}e^{i(\omega-m\Omega)t}$ and $\tilde{G}^{<}_{mn}(\bm{k};\omega,\omega^{\prime})\equiv\tilde{G}^{<}(\bm{k};\omega-m\Omega,\omega^{\prime}-n\Omega)$ ($0\leq\omega,\omega^{\prime}\leq\Omega$). For convenience, we extend the domain of $\tilde{G}^{<}_{mn}(\bm{k};\omega,\omega^{\prime})$ by $\tilde{G}^{<}_{m+1,n+1}(\bm{k};\omega,\omega^{\prime})=\tilde{G}^{<}_{mn}(\bm{k};\omega-\Omega,\omega^{\prime}-\Omega)$. We note that $m$ and $n$ represent the Floquet indices.

We assume $G^{<}(\bm{k};t,t^{\prime})=G^{<}(\bm{k};t+T,t^{\prime}+T)$, meaning that the Green function is periodic with respect to the center of time $(t+t^{\prime})/2$. This assumption implies that we consider a nonequilibrium steady state represented by a mixed state of the Floquet eigenstates. From this assumption, we have $\tilde{G}^{<}_{mn}(\bm{k};\omega,\omega^{\prime})=2\pi\delta(\omega-\omega^{\prime})G^{<}_{mn}(\bm{k},\omega)$, and then the Green function at the same time $t=t^{\prime}$ is given by

[TABLE]

Replacing the dummy variables of summation $(m,n)$ with $(l,n)$ ($l=m-n$) and using the relation $G^{<}_{m+1,n+1}(\bm{k},\omega)=G^{<}_{mn}(\bm{k},\omega-\Omega)$, we obtain

[TABLE]

As a result, the Fourier component of $\bm{j}(t)$ at frequency $n\Omega$ is obtained as

[TABLE]

This is the formula by which we calculate the HHC from the non-equilibrium Green function.

Appendix B Calculation of HHC

Here, we concretely calculate the $O(U^{2})$ correction of the HHC and derive Eqs. (15)-(18). For convenience, we repeat the calculation procedure in the main text. In the absence of the lattice potential, we obtain the exact bare Green function as follows,

[TABLE]

where we have used $\sum_{M}\ket{\psi^{M}_{0}(\bm{k})}\bra{\psi^{M}_{0}(\bm{k})}=1$, $\braket{\phi^{m}(\bm{k})}{\psi^{M}_{0}(\bm{k})}=J_{M-m}(\alpha_{\bm{k}})$, and $H^{\text{F}}_{0}(\bm{k})\ket{\psi^{M}_{0}(\bm{k})}=\epsilon_{\bm{k}}^{M}\ket{\psi^{M}_{0}(\bm{k})}$. In the presence of the potential, the exact full Green function is not available. However, the lowest-order correction is obtained by the second-order Feynmann rules for the nonequilibrium Green function Rammer and Smith (1986),

[TABLE]

where both $G_{mn}^{(2)}(\bm{k},\omega)$ and $g_{mn}(\bm{k},\omega)$ are the $2\times 2$ matriices,

[TABLE]

Figure 7 shows the Feynmann diagram of $G^{(2)}$. Once we obtain the Green function, we can calculate the high harmonics from Eq. (25) with $G^{<}=(G^{A}-G^{R}+G^{K})/2$.

Let us concretely calculate the second-order correction of the Green function from Eq. (26)-(28). From Eq. (28), the retarded, the advanced and the Keldysh components are

[TABLE]

where we have defined the following abbreviations:

[TABLE]

Therefore, the HHC at frequency $n\Omega$ is

[TABLE]

where we have used $G^{<}=(G^{A}-G^{R}+G^{K})/2$.

We perform the $\omega$-integration in Eq. (35) by calculating the residues. Since the integrand vanishes fast enough for $|\omega|\rightarrow\infty$, we replace the integral path with the closed semicircle in the upper or the lower halves of the $\omega$-plane. First, the $\omega$-integration of $AAA$ and $RRR$ terms vanish because their integrands have poles only in the upper or the lower halves of the $\omega$-plane (see Figs. 8(a) and (d)). Next, we perform the $\omega$-integrations of $RRK,KAA$ and $RKA$ terms. From Eqs. (26) and (27), we calculate

[TABLE]

where we have calculated the residues owing to the second equality (see Fig. 8(b)). In the same way, we obtain

[TABLE]

Figure 8 shows the poles of the each term in Eq. (35). Summing up all these terms, we obtain the following expression

[TABLE]

where we have used $\sum_{a}J_{a}(x)J_{b-a}(y)=J_{b}(x+y)$ and $V_{\bm{k}}=U\sum_{a}u_{\bm{k}}e^{-i\bm{k}\cdot\bm{r}_{a}}$. Taking the average over the random variables $\delta\bm{r}_{a}$’s [see Eq. (14) in the main text] and substituting Eq. (39) and (40) to Eq. (35), 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 $\tau=48$ $\mathrm{fs}$ and $24$ $\mathrm{fs}$. The panel (a) is the same as Fig. 3(a) in the main text, and the panel (b) is calculated with $\tau$ 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.