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,

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

 

Fig. 1. (a) Schematic pictorial description of the generic Feynman trajectories contributing to the dipole moment. (b) Generic Feynman trajectories satisfying the stationary-phase conditions Eq. (15). The curves are not solutions of the Newton equations of motion but only a pictorial description of the Green functions ${\bar G_c}$ (blue asterisks), ${G_v}$ (red circles), $G_v^*$ (solid black line). The vertical dotted line represents the time instant $t^\prime $ when the conduction-band electron is created.

Download Full Size | PDF

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.

 

Fig. 2. HHG spectrum generated by a disordered semiconductor with ${\Gamma _1} = 0.25\hbar {\omega _L}$, interacting with a laser field with $I = 2.0\;{\rm TW}/{{\rm cm}^2}$ and $\lambda = 2.5 \;{\unicode{x00B5}{\rm m}}$. Blue line is the total spectrum; red line is the intraband contribution.

Download Full Size | PDF

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}$.

 

Fig. 3. As in Fig. 2, but with ${\Gamma _1} = 0.7\hbar {\omega _L}$.

Download Full Size | PDF

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.

 

Fig. 4. Detailed view of the HHG spectrum as a function of the disorder level in the solids. The topmost curve (black solid line with crosses) is the HHG spectrum generated by a crystal without disorder ${\Gamma _1} = 0.0$, the middle curve (dashed red) is the spectrum emitted by a solid with disorder strength ${\Gamma _1} = 0.25\hbar {\omega _L}$, and the bottom curve (blue line) is the spectrum radiated by the solid when ${\Gamma _1} = 0.7\hbar {\omega _L}$. Intensity and frequency of the laser are as in Fig. 2. The three curves have been shifted by arbitrary factors to improve visibility.

Download Full Size | PDF

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.