High order harmonic generation in solids: a review on recent numerical methods

ABSTRACT

Interaction of intense lasers with solid materials offers an alternative way to achieve high-order harmonic generation (HHG). Since the underlying mechanisms of the harmonic emission remain uncertain, a large number of theoretical studies have been proposed to explain the experimental findings in this fast-growing field. Here, we review the primary numerical methods and present a brief perspective of HHG in solids.

Graphical Abstract

KEYWORDS:

• Ultrafast laser
• strong field
• high-order harmonic generation
• solid materials
• numerical simulations
• laser-crystal interaction

1. Introduction

Beyond the investigations on high-order harmonic generation (HHG) in gaseous media, recent experimental efforts have demonstrated that such a nonperturbative nonlinear HHG phenomenon can be observed from solid materials in strong laser fields [1–10], which offers new opportunities to obtain efficient extreme ultraviolet sources [7,11–13] and study ultrafast electronic dynamics in the condensed phase [14–16].

Generally speaking, two primary sources contribute to solid HHG: interband polarization and intraband current. Analogous to the three-step model in atomic or molecular gases, the involved interband and intraband processes in a crystal solid can be described as three stages [17]: (i) electron tunneling excitation from the valence band to the conduction band, (ii) electron (hole) acceleration on the conduction (valence) band, which will lead to an intraband current, and (iii) electron–hole recombination gives rise to a time-dependent interband polarization. However, solid and atomic HHG exhibit a different cutoff law, which is determined by the maximum kinetic energy an electron can obtain while being accelerated by the electric field. In experiment and theory [1,4,7,18,19], the cutoff is found as linear dependence with laser wavelength and field strength in solid HHG instead of a quadratic scale in atomic case. In addition, theoretical studies predict that the maximum photon energy (i.e. cutoff energy) in HHG spectrum is limited by the largest transition energy among the different energy bands in solids [17,20,21]. In general, the cutoff energy presents a linear dependence with laser wavelength and field strength dominantly in intraband HHG and is limited to the largest energy difference among energy bands in interband HHG, due to the different generation mechanisms in intraband and interband HHG.

In 2011, Ghimire et al. [1,2] performed pioneering experiments using ZnO crystal in mid-infrared (MIR) laser field and considered intraband Bloch oscillations with the electrons and holes reflected at an edge of the Brillouin zone as the mechanism of HHG. Subsequently, solid HHG experiments were also accomplished in longer wavelengths toward THz fields [3,4], and intraband HHG is found to be dominant in long wavelength regime. However in the experiment, Vampa et al. [5] claimed that the interband polarization created by electron–hole recollision is the main source of solid HHG for MIR driving pulses by adding a weak second-harmonic beam to study the emission process. Moreover, Hohenleutner et al. [6] measured HHG in a bulk solid directly in the time domain and found electrons with tunneling interference among different bands play an important role in ultrafast harmonic radiation. Besides, Luu et al. [7] applied a sub-cycle synthesized field with bandwidths that cover the ultraviolet and near-infrared spectral regions to a silicon dioxide crystal and showed that coherent extreme ultraviolet radiation is emitted with photon energies about 40 eV. In their work, anharmonic electron motion is used to explain the generation of high-frequency component when electrons are driven with oscillation in lattice by the field. They also observed that the total photon yield of fused silica is higher than the noble gases.

In order to find more potential solid media for HHG, some special materials like low dimensional [22–31] such as graphene [22–25], monolayer molybdenum disulfide (MoS₂) [22,26], hexagonal boron nitride (h-BN) [27,28] or strongly correlated systems [16,32,33] have attracted considerable attention. For example, Yoshikawa et al. [22] observed that HHG signal in graphene is dramatically enhanced under elliptically polarized excitation. Compared to MoS₂ bulk, Liu et al. [26] found that the harmonic efficiency per layer exhibits significant enhancement in layered structures. Tancogne-Dejean et al. [34] also predicted that monolayer and few layer materials are promising for efficient HHG because of their strong inhomogeneity of the electron–nuclei potential. For strongly correlated systems such as the Mott insulator [16,32,33], Silva et al. [16] found that high-harmonic generation can be used to time-resolve ultrafast non-equilibrium many-body dynamics. Murakami et al. [33] demonstrated that the Mott insulator shows a stronger high-harmonic intensity than a semiconductor model with the same dispersion, even if the semiconductor bands are broadened by impurity scattering to mimic the incoherent scattering in the Mott insulator.

Nowadays, the complete waveform characterization of sub-cycle driver pulse with a precisely defined carrier-envelope phase (CEP) [35–38] makes it possible to coherently control the sub-cycle dynamics in solids. Shirai et al. [38] reported the HHG in a Si membrane driven by CEP-controlled sub-cycle MIR pulses generated through two-color filamentation. They found that the high-harmonic spectrum reaching the ultraviolet region shows clear evidence of the CEP dependence and the complex CEP dependence can be explained by the interference between different orders of harmonics. In addition, plasmonic high-harmonic generation [39,40] drew attention as a means of producing coherent extreme ultraviolet radiation in solids. To date, orientation-dependent high harmonic spectrum in solids has also become the focus both in experiment [41–45] and theory [34,46–49]. You et al. [42] found that the harmonic signal is strongly enhanced or diminished when the electron wavepacket is directed into or away from neighboring atomic sites. In fact, the discovery of HHG in solids provides an all-optical technique to reconstruct the Berry curvature [45] and band structure [50,51].

Accompanied by the rapid experimental progresses in solid HHG, many theoretical works [4,6,14,17–21,29,30,34,46–50,52–93] have been intensively developed. In this review, we focus on the primary and widely used numerical methods based on time-dependent Schrödinger equation (TDSE) [18,46,52–65], semiconductor Bloch equations (SBE) [4,6,11,14,17,20,21,29,47,48,50,66–78] and time-dependent density-functional theory (TDDFT) [30–32,34,81–84]. Keep in mind that it is merely a subset of rapidly emerging theoretical studies about solid HHG. Last but not the least, we summarize current status and some challenges of present numerical simulations which are yet to be improved or solved.

2. Numerical method by solving the time-dependent Schrödinger equation

As early as in 1992, numerical computation on solution of TDSE was performed on the basis of Bloch functions by Plaja et al. [52] to treat HHG in a crystalline solid. In their model, the electron wave function in the crystal is obtained from an appropriately chosen pseudopotential for matching the case of pure silicon, and harmonic generation proves to be more feasible when laser frequencies are below gap resonance. Then, Bachau et al. [94] investigated the irradiation of wide band-gap insulators with a femtosecond Ti:Sapphire laser pulse by using a quantum mechanical approach based on the solution of TDSE and analyzed in detail the heating mechanism due to the sequence of direct interbranch transitions in the conduction band. Moreover, Korbman et al. [95] studied the polarization response of a model dielectric resembling SiO₂ by solving TDSE under external field and they point out that interband transitions create coherent superpositions of states that manifest themselves as quantum beats in the polarization. The transitions between different conduction bands play an important role in the generation of ionization [96] and HHG [53] in periodic solids.

For the insights into the HHG process, TDSE was solved by two numerical approaches based on different state bases (Bloch state basis and Houston state basis) for an electron interacting with a periodic potential [18]. We briefly review these two state bases of TDSE model below.

2.1. Bloch state basis

Usually, linearly polarized laser field is considered to propagate through the crystal along the optical axis in the simulation. When laser-solid interaction is simplified in one dimension, TDSE in velocity-gauge reads as:

(1) $ih\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left|\psi (t)⟩=\left(\frac{{\hat{p}}^2}{2m}+V(\hat{x})+\frac{e}{m}A(t)\hat{p}\right)\right|\psi (t)⟩,$(1)

where the vector potential A(t) is employed by the dipole approximation A(x, t) ≈ A(t) because the wavelengths we are interested in are much larger than the lattice constant, and is related to the electric field by:

(2) $A(t)=-\underset{-\mathrm{\infty }}{\overset{t}{\int }}E(t^{\prime })d{t}^{\prime }.$(2)

$\hat{p}$ is the momentum operator and V(x) is the periodic lattice potential. According to Bloch’s theorem, each Bloch state can be written as a product of a plane wave and a periodic function:

(3) $⟨x|{\varphi }_k^n⟩=e^{ikx}{U}_k^n(x),$(3)

where U_kⁿ (x) satisfy $U_k^n(x+a_0)=U_k^n(x).$ For each value of the quasimomentum k, the electron wave function is represented in the basis of Bloch states:

(4) $|{\psi }_k(t)⟩=\sum _n{C}_k^n(t)|{\varphi }_k^n⟩,$(4)

where Bloch states $|{\varphi }_k^n⟩$ are evaluated by solving the single-electron stationary Schrödinger equation:

(5) $\left[\frac{{\hat{p}}^2}{2m}+V(\hat{x})\right]|{\varphi }_k^n⟩=E_n(k)|{\varphi }_k^n⟩,$(5)

where n is the band index and the eigenvalues E_n(k) represent the dispersion relations of the bands. Finally, the time-dependent laser-induced macroscopic current j(t) can be expressed by the sum of the current in each of the different k channel:

(6) $j(t)=-\frac{e}{m}\int \left\{\mathrm{R}\mathrm{e}\left[⟨{\psi }_k|\hat{p}|{\psi }_k⟩\right]+eA(t)\right\}dk.$(6)

The harmonic spectra can be obtained by the Fourier transformation of the macroscopic current.

2.2. Houston state basis

The Houston state [97] is the best as an adiabatic basis in which the lattice momentum that would be ћk₀ in the absence of a field, which has the time dependence [18]:

(7) $\mathrm{h}k(t)=\mathrm{h}{k}_0+eA(t).$(7)

In the Houston approach, TDSE reads as:

(8) $i\mathrm{h}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left|\psi (t)⟩=\left[\frac{{(\hat{p}+eA)}^2}{2m}+V(\hat{x})\right]\right|\psi (t)⟩.$(8)

The Houston states are related to the Bloch states with lattice momentum ћk(t) by [98]:

(9) $\left|{\tilde{\varphi }}_{k_0}^n(t)\right.⟩=e^{-ieA\hat{x}/\mathrm{\hslash }}\left|{\varphi }_{k(t)}^n\right.⟩,$(9)

where $\hat{x}$ is the position operator. The time-dependent wave function with initial lattice momentum ћk₀ is expressed in Houston states:

(10) $|\psi (t)⟩=\sum _n{a}_{k_0}^n(t)|{\tilde{\varphi }}_{k_0}^n(t)⟩.$(10)

Finally, the total current can be calculated:

(11) $j_{tot}=-\frac{e}{m}Re\left[\sum _{n{n}^{\prime }}{a}_{n{k}_0}^{\ast }{a}_{n^{\prime }{k}_0}⟨{\varphi }_{{ }_{k(t)}}^n|\hat{p}|{\varphi }_{{ }_{k(t)}}^{n^{\prime }}⟩\right].$(11)

The advantage of the Houston basis is that the electron dynamics naturally separates into an intraband and interband contributions:

(12) $j_{intra}=-\frac{e}{m}\sum _n{\left|a_{n{k}_0}\right|}^2⟨{\varphi }_{{ }_{k(t)}}^n|\hat{p}|{\varphi }_{{ }_{k(t)}}^n⟩,$(12)

(13) $j_{inter}=-\frac{e}{m}Re\left[\sum _{n{n}^{\prime }(n\ne n^{\prime })}{a}_{n{k}_0}^{\ast }{a}_{n^{\prime }{k}_0}⟨{\varphi }_{{ }_{k(t)}}^n|\hat{p}|{\varphi }_{{ }_{k(t)}}^{n^{\prime }}⟩\right].$(13)

2.3. Application examples of TDSE method

Based on these two state bases of TDSE model, Wu et al. [18] showed that high harmonic spectrum in solid exhibits multiple plateaus and they declared that the primary plateau is due to transitions between the valence band and the lowest conduction band, whereas the secondary plateau and more generally higher frequencies in the spectrum are due to contributions from higher-lying bands. In Ref. [18], the authors compared HHG spectra calculated in the three-band Houston basis and the 51-band Bloch basis at different field strengths 4.5 × 10¹¹ W/cm² shown in Figure 1(a) and 1.0 × 10¹² W/cm² in Figure 1(b). The initial condition used in both Houston basis and Bloch basis is a single k state (k = 0) in the valence band where the band gap is the smallest. It can be found that the Houston spectrum matches that of the Bloch spectrum very well in lower laser intensity as shown in Figure 1(a). However, since the Houston calculation only includes three bands, in Figure 1(b) it cannot reproduce the high-frequency part of the spectrum due to higher-lying bands from the Bloch calculation in higher laser intensity.

Figure 1. HHG from a three-band Houston basis calculation and a 51-band Bloch basis calculation, for intensities of a 3.2-μm driving field: (a) 4.5 × 10¹¹ W/cm² and (b) 1.0 × 10¹² W/cm². The black thin line is calculated from the conduction band only in the Houston basis. Adapted from Ref. [18].

Soon afterwards, Ndabashimiye et al. [54] experimentally compared HHG from rare gases in the solid and gas phase and confirmed the multiple-plateau structure of HHG in the rare-gas solids. In detail, owing to the weak van der Waals interaction, rare-gas solids are a near-ideal medium as the role of high density and periodicity to study HHG process in the solid state. In their work [54], HHG from solid argon and krypton has been studied using different laser intensities and wavelengths. Figure 2(a) shows the high-harmonic spectra measured in solid argon as a function of laser intensity at the 1333 nm pump laser. It is clear to see that the second plateau appears when the laser intensity increases. Similar phenomenon was observed for solid krypton. Experimental results are successfully reproduced by solving TDSE in Houston state basis [54]. A four-band model was performed to simulate the HHG in solid argon. Figure 2(b) shows the calculated harmonic spectrum for comparison in their work. The energies for the relevant bands at Γ point are illustrated in the inset. TDSE simulation indicates that the harmonic strength and the cutoff energy of each plateau depend on the energy separation and the coupling strength between different pairs of bands. For example, the first plateau originates from the coupling between levels 1 and 2, and the second plateau from the coupling between 2 and 3. Meanwhile, the cutoff energy of the second plateau (dashed white curve) is determined by the energy difference between the field-dressed levels 1 and 3 [54]. More details and explanations of the TDSE method they used can be found in Ref. [55].

Figure 2. The measured HHG spectrum of solid argon as function of laser intensity using 1333 nm drive laser. (b) Simulated HHG spectrum by solving four-band TDSE model. The energies for the relevant bands at Γ point are illustrated in the inset. The dashed rectangle shows the corresponding range of the experimental data in (a). The dashed white curve indicates the predicted cutoff energy of the second plateau based on the energy difference between field-dressed levels 1 and 3. Adapted from Ref. [54].

In addition to Bloch basis and Houston basis, TDSE can be numerically solved by using B-spline basis sets in coordinate space and the harmonic spectra obtained agree well with the results simulated by TDSE in k-space using Bloch states [56]. Recently, to simulate the orientation-dependent high harmonic spectrum in solids, two-dimensional (2D) TDSE [46,57,58] was solved by exploring HHG from a 2D square lattice which is described by the 2D potential V(x, y). By comparing with the 1D case, Jin et al. [58] found that the generation dynamics have a significant difference due to the existence of many crossing points in the 2D band structure. In particular, the higher conduction bands can be excited step by step via these crossing points and the total contribution of the harmonic is given by the mixing of transitions between different clusters of conduction bands to the valence band [58].

3. Numerical method by solving semiconductor Bloch equations

The dynamics of strong lasers interacting with semiconductors can be theoretically described by the SBE [98–100]. Recently, SBE model has been successfully used to investigate the properties of harmonic emission and electric current in solids. In this section, we introduce the SBE method as well as its improvement for solid HHG study.

3.1. Two-band model

The coupled interband and intraband dynamics of semiconductors under strong laser field described by the SBEs within a two-band model read [20,66–72]:

(14) $i\frac{\mathrm{\partial }}{\mathrm{\partial }t}{p}_k=\left({\epsilon }_c(k)-{\epsilon }_v(k)-i\frac{1}{T_2}\right)p_k+\left(f_k^c-f_k^v\right)d(k)E(t)+iE(t){\mathrm{\nabla }}_k{p}_k,$(14)

(15) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}{f}_k^{c(v)}=-2\mathrm{I}\mathrm{m}\left[d(k)E(t)p_k^{\ast }\right]+E(t){\mathrm{\nabla }}_k{f}_k^{c(v)}.$(15)

In Equations (14) and (15), f_k^c and f_k^v are the populations of the lowest conduction band and the highest valence band, respectively. Electrons are initially filled in valence band. p_k is microscopic interband polarization between conduction band and valence band. d(k) is the dipole element, and E(t) is the temporal laser field. T₂ is the dephasing time. ε_c(k) and ε_v(k) are k-dependent energy bands of the lowest conduction band and the highest valence band, respectively, that are usually calculated from first-principles calculations.

Because of the motion of the carriers within the bands under laser field, the intraband electric current J(t) and macroscopic interband polarization P(t) are given by:

(16) $J(t)=\sum _{\lambda }{\int }_{BZ}e{\upsilon }_k^{\lambda }{f}_k^{\lambda }(t)dk,$(16)

(17) $P(t)={\int }_{BZ}\left[d(k)p_k(t)+c.c.\right]dk,$(17)

where υ_k^λ is the group velocity generated from the derivative of the bands and λ is the band index.

Based on two-band SBEs, Golde et al. [20,67,68] investigated the emitted radiation of a semiconductor nanostructure after excitation with an extremely intense ultrashort laser pulse [67,68] and THz field [20]. They found the nontrivial coupling of interband and intraband dynamics leads to efficient HHG. In particular, the intraband acceleration significantly modifies the dynamics of the interband polarization which results in a strong enhancement of HHG. Yu et al. [70] paid attention to the importance of the dipole moment between valence and conduction bands, and found that the shape of k-space dependent dipole moment plays an important role in harmonic emission, because it determines the interband transition process [101] in solid HHG. They improved the SBE model and performed the simulations utilizing real dipole moment calculated from state-of-the-art first-principles code.

Later works [71,72] also established the important roles of the transition dipole and its phase in HHG from solids. They found that the strength [71] and orientation dependence [72] of even harmonics are attributed to the variation of transition dipole phase from excitation to recombination. In Ref [72], the authors have taken into account the transition dipole phase calculated by tight-binding model and employ the familiar SBEs method within the simplest 1D two-band model to calculate the angle-dependent HHG from ZnO. The simulation results are shown in Figure 3(b) and compared to Figure 3(a) from the experiment of Ref [43]. One can see that the agreement between the two-band model and the experiment is quite good. Even harmonics appear at angles away from θ = 90° as expected.

Figure 3. (a) Orientation-dependent HHG spectra, redrawn from the experimental data in Ref [43]. (b) Spectra calculated by considering the transition dipole phase from the SBEs method within the simplest 1D two-band model. Adapted from Ref. [72].

Similar to the two-band SBE model as described above, a three-dimensional (3D) single active electron model called two-band density matrix equations (DME) [17,29,50,73–76] is suggested and numerically solved to study the HHG in bulk crystals when exposed to intense lasers:

(18) $\dot{\pi }(\mathbf{K},t)=-\frac{\pi (\mathbf{K},t)}{T_2}-i\mathrm{\Omega }(\mathbf{K},t)w(\mathbf{K},t)e^{-iS(\mathbf{K},t)},$(18)

(19) ${\dot{n}}_m(\mathbf{K},t)=i{s}_m{\mathrm{\Omega }}^{\ast }(\mathbf{K},t)\pi (\mathbf{K},t)e^{iS(\mathbf{K},t)}+c.c.,$(19)

where n_m is the valence (m = v) and conduction (m = c) band population, and w = n_v – n_c is the population difference. Initially all electrons are in the valence band; s_m = −1, 1 for m = v, c, respectively. $S(\mathbf{K},t)={\int }_{-\mathrm{\infty }}^t{\epsilon }_g\left[\mathbf{K}+\mathbf{A}(t^{\prime })\right]d{t}^{\prime }$ is the classical action, $\mathrm{\Omega }(\mathbf{K},t)=\mathbf{F}(t)\mathbf{d}\left[\mathbf{K}+\mathbf{A}(t)\right]$ is the Rabi frequency, T₂ is the dephasing time, and d(k) is the transition dipole moment. In essence, the main difference between two-band DME and SBE models is just the frame transformation [17]. In DME method, HHG in solids is also determined by intraband j_ra and interband j_er contribution. Their expressions for the current averaged over the unit cell are given by:

(20) ${\mathbf{j}}_{ra}(t)=\sum _{m=c,v}{\int }_{\overline{BZ}}{\mathbf{v}}_m\left[\mathbf{K}+\mathbf{A}(t)\right]n_m(\mathbf{K},t)d^3\mathbf{K},$(20)

(21) ${\mathbf{j}}_{er}(t)=\frac{d}{dt}{\int }_{\overline{BZ}}\mathbf{d}\left[\mathbf{K}+\mathbf{A}(t)\right]\pi (\mathbf{K},t)e^{iS(\mathbf{K},t)}{d}^3\mathbf{K}+c.c.,$(21)

where v_m is the band velocity defined by ${\mathbf{v}}_m(\mathbf{k})={\mathrm{\nabla }}_{\mathbf{k}}{E}_m(\mathbf{k})$.

Based on two-band DME, Vampa et al. [17] elucidated two main mechanisms that interband polarization and intraband current contribute to the HHG in ZnO (wurtzite structure), and they emphasized that interband current is the dominant mechanism for HHG. In detail, Figure 4 shows harmonic spectrum from interband (blue line) and intraband (red line) currents [17]. The orientation of the reciprocal lattice of ZnO is chosen so that x‖Γ–M, y‖Γ–K, and z‖Γ–A (optical axis). The used driving laser is a temporal Gaussian envelope with a FWHM equal to 10 cycles. The field strength is F₀ = 0.003 a.u. and frequency is ω₀ = 0.014 a.u. Clearly, HHG can be divided into two regions, i.e. perturbative regime and nonperturbative regime (plateau region) by the energy of minimum band gap. The HHG in perturbative regime takes place below the dashed line in Figure 4, the two contributions become comparable with intraband HHG being slightly stronger. Whereas in the nonperturbative regime, interband HHG is always dominant by at least 2 orders of magnitude. Moreover, one can see that spectra in Figure 4(a) are very noisy and do not show a clear odd harmonic structure in plateau region, because of the interference of multiple returns with full dephasing time T₂ = ∞. Then, the odd harmonic spectrum starts to manifest when reducing the dephasing times to half cycle in Figure 4(b) as higher returns are suppressed. For even shorter dephasing time in Figure 4(c), a clean odd harmonic spectrum emerges.

Figure 4. Simulated harmonic spectrum of ZnO crystal from interband (blue line) and intraband (red line) currents; (a), (b), (c) show dephasing times T₂ = ∞, T₀/2, T₀/4, respectively. T₀ = 2π/ω₀ = 10.9 fs. The dashed black vertical lines represent the (minimum) band gap at the center of the Brillouin zone. Adapted from Ref. [17].

Moreover, Vampa et al. found that interband and intraband HHG exhibit different wavelength dependence, and interband emission is dominant for the MIR laser, whereas intraband emission dominates the far-infrared range [73]. Then, interband Bloch oscillations were identified as the mechanism for the higher plateaus in solid HHG from DME model [74]. Liu et al. [75] theoretically investigated HHG in ZnO by DME model with elliptically polarized MIR light, and they demonstrated that the sensitivity for the ellipticity dependence of the interband harmonics above the band gap is determined by both the harmonic order and the intensity of driving light. They also pointed out HHG from elliptically polarized driving field carries the significant signature of the band structure in solids [75]. Interestingly, the effect of quantum confinement on HHG in semiconductors was investigated by systematically varying the width of a model quantum nanowire in DME model [29], indicating that low-dimension systems present a potential pathway to increase yield and photon energy of HHG in solids.

3.2. Multi-band model

For better describing the time evolution of polarizations and carrier occupations of solids under intense field, the multi-band SBE [4,6,11,21,47,48,78] read:

(22) $\begin{array}{rl}& i\frac{\mathrm{\partial }}{\mathrm{\partial }t}{p}_{\mathbf{k}}^{h_i{e}_j}=\left({\epsilon }_{\mathbf{k}}^{e_j}+{\epsilon }_{\mathbf{k}}^{h_i}-i\frac{1}{T_2}\right)p_{\mathbf{k}}^{h_i{e}_j}-\left(1-f_{\mathbf{k}}^{e_j}-f_{\mathbf{k}}^{h_i}\right){\mathbf{d}}_{\mathbf{k}}^{e_j{h}_i}\mathbf{E}(t)+i\mathbf{E}(t){\mathrm{\nabla }}_{\mathbf{k}}{p}_{\mathbf{k}}^{h_i{e}_j}\\ & \mathbf{E}(t)\sum _{e_{\lambda }\ne e_j}\left({\mathbf{d}}_{\mathbf{k}}^{e_{\lambda }{h}_i}{p}_{\mathbf{k}}^{e_{\lambda }{e}_j}-{\mathbf{d}}_{\mathbf{k}}^{e_j{e}_{\lambda }}{p}_{\mathbf{k}}^{h_i{e}_{\lambda }}\right)+\mathbf{E}(t)\sum _{h_{\lambda }\ne h_i}\left({\mathbf{d}}_{\mathbf{k}}^{h_{\lambda }{h}_i}{p}_{\mathbf{k}}^{h_{\lambda }{e}_j}-{\mathbf{d}}_{\mathbf{k}}^{e_j{h}_{\lambda }}{p}_{\mathbf{k}}^{h_i{h}_{\lambda }}\right),\end{array}$(22)

(23) $\begin{array}{rl}& i\frac{\mathrm{\partial }}{\mathrm{\partial }t}{p}_{\mathbf{k}}^{e_i{e}_j}=\left({\epsilon }_{\mathbf{k}}^{e_j}-{\epsilon }_{\mathbf{k}}^{e_i}-i\frac{1}{T_2}\right)p_{\mathbf{k}}^{e_i{e}_j}+\left(f_{\mathbf{k}}^{e_j}-f_{\mathbf{k}}^{e_i}\right){\mathbf{d}}_{\mathbf{k}}^{e_j{e}_i}\mathbf{E}(t)+i\mathbf{E}(t){\mathrm{\nabla }}_{\mathbf{k}}{p}_{\mathbf{k}}^{e_i{e}_j}\\ & +\mathbf{E}(t)\sum _{e_{\lambda }\ne e_j}{\mathbf{d}}_{\mathbf{k}}^{e_{\lambda }{e}_i}{p}_{\mathbf{k}}^{e_{\lambda }{e}_j}-\mathbf{E}(t)\sum _{e_{\lambda }\ne e_i}{\mathbf{d}}_{\mathbf{k}}^{e_j{e}_{\lambda }}{p}_{\mathbf{k}}^{e_i{e}_{\lambda }}\\ & +\mathbf{E}(t)\sum _{h_{\lambda }}\left({\mathbf{d}}_{\mathbf{k}}^{h_{\lambda }{e}_i}{p}_{\mathbf{k}}^{h_{\lambda }{e}_j}-{\mathbf{d}}_{\mathbf{k}}^{e_j{h}_{\lambda }}{(p_{\mathbf{k}}^{h_{\lambda }{e}_i})}^{\ast }\right),\end{array}$(23)

(24) $\begin{array}{rl}& i\frac{\mathrm{\partial }}{\mathrm{\partial }t}{p}_{\mathbf{k}}^{h_i{h}_j}=\left({\epsilon }_{\mathbf{k}}^{h_i}-{\epsilon }_{\mathbf{k}}^{h_j}-i\frac{1}{T_2}\right)p_{\mathbf{k}}^{h_i{h}_j}+\left(f_{\mathbf{k}}^{h_i}-f_{\mathbf{k}}^{h_j}\right){\mathbf{d}}_{\mathbf{k}}^{h_j{h}_i}\mathbf{E}(t)+i\mathbf{E}(t){\mathrm{\nabla }}_{\mathbf{k}}{p}_{\mathbf{k}}^{h_i{h}_j}\\ & +\mathbf{E}(t)\sum _{h_{\lambda }\ne h_j}{\mathbf{d}}_{\mathbf{k}}^{h_{\lambda }{h}_i}{p}_{\mathbf{k}}^{h_{\lambda }{h}_j}-\mathbf{E}(t)\sum _{h_{\lambda }\ne h_i}{\mathbf{d}}_{\mathbf{k}}^{h_j{h}_{\lambda }}{p}_{\mathbf{k}}^{h_i{h}_{\lambda }}\\ & +\mathbf{E}(t)\sum _{e_{\lambda }}\left({\mathbf{d}}_{\mathbf{k}}^{e_{\lambda }{h}_i}{(p_{\mathbf{k}}^{h_j{h}_{\lambda }})}^{\ast }-{\mathbf{d}}_{\mathbf{k}}^{h_j{e}_{\lambda }}{p}_{\mathbf{k}}^{h_i{e}_{\lambda }}\right),\end{array}$(24)

(25) $\begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}{f}_{\mathbf{k}}^{e_i}=-2\mathrm{I}\mathrm{m}\left[\sum _{e_{\lambda }\ne e_i}{\mathbf{d}}_{\mathbf{k}}^{e_i{e}_{\lambda }}\mathbf{E}(t){\left(p_{\mathbf{k}}^{e_{\lambda }{e}_i}\right)}^{\ast }+\sum _{h_{\lambda }}{\mathbf{d}}_{\mathbf{k}}^{e_i{h}_{\lambda }}\mathbf{E}(t){\left(p_{\mathbf{k}}^{h_{\lambda }{e}_i}\right)}^{\ast }\right]\\ +\mathbf{E}(t){\mathrm{\nabla }}_{\mathbf{k}}{f}_{\mathbf{k}}^{e_i}-\frac{1}{2T_1}\left(f_{\mathbf{k}}^{e_i}-f_{-\mathbf{k}}^{e_i}\right),\end{array}$(25)

(26) $\begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}{f}_{\mathbf{k}}^{h_i}=-2\mathrm{I}\mathrm{m}\left[\sum _{h_{\lambda }\ne h_i}{\mathbf{d}}_{\mathbf{k}}^{h_{\lambda }{h}_i}\mathbf{E}(t){\left(p_{\mathbf{k}}^{h_i{h}_{\lambda }}\right)}^{\ast }+\sum _{e_{\lambda }}{\mathbf{d}}_{\mathbf{k}}^{e_{\lambda }{h}_i}\mathbf{E}(t){\left(p_{\mathbf{k}}^{h_i{e}_{\lambda }}\right)}^{\ast }\right]\\ +\mathbf{E}(t){\mathrm{\nabla }}_{\mathbf{k}}{f}_{\mathbf{k}}^{h_i}-\frac{1}{2T_1}\left(f_{\mathbf{k}}^{h_i}-f_{-\mathbf{k}}^{h_i}\right).\end{array}$(26)

Here, f_k^e and f_k^h are the populations of electrons and holes, respectively. We assume that electrons are initially filled in valence bands. p_k^λλ′ is microscopic interband polarizations between conduction bands and valence bands. λ = e, h is the index, which specifies either an electron or a hole. T₁ represents the phenomenological damping of the antisymmetric part of the carriers. T₂ is the dephasing time. ε_k^e = ε_k^c and ε_k^h = – ε_k^v are k-dependent energy of the corresponding carriers in conduction or valence bands. Similar to the two-band model, energy bands (i, j) and transition dipole elements d_k^λλ′ are calculated from first-principles calculations.

Due to the motion of the carriers within the bands under laser field, the total time-dependent macroscopic interband polarization P(t) and intraband electric current J(t) are given by:

(27) $P(t)=\sum _{\lambda ,{\lambda }^{\prime },\mathbf{k}}\left[{\mathbf{d}}_{\mathbf{k}}^{\lambda {\lambda }^{\prime }}{p}_{\mathbf{k}}^{\lambda {\lambda }^{\prime }}(t)+c.c.\right],$(27)

(28) $J(t)=\sum _{\lambda ,\mathbf{k}}\left[-2{\upsilon }_{\mathbf{k}}^{\lambda }{f}_{\mathbf{k}}^{\lambda }(t)\right],$(28)

where υ_k^λ is the group velocity which can be calculated from the derivative of the bands. Finally, the total intensity of high-harmonic spectrum can be obtained by:

(29) $S(\omega )\propto {\left|\omega \mathbf{P}(\omega )+i\mathbf{J}(\omega )\right|}^2.$(29)

In multi-band SBE model, the quantum interference from different excitation path among each pair of bands is thought to be very important and responsible for even harmonic radiation [4,47,48,78]. In Figure 5, taking a three-band system as an example, different excitation path from valence band h₁ to the conduction band e can be distinguished. As discussed in Ref [78]., a direct excitation from the valence band h₁ to the conduction band e (blue arrow in Figure 5(b) consists of a series of terms containing only odd powers of E(t) result in resonances at odd multiples of the driving field frequency. Thus, only odd harmonic orders are generated in the process of direct excitation. However, if E(t) is strong enough to create vacancies in h₂, transitions between both valence bands become possible. As a result, an additional indirect excitation path, h₁ to h₂ to e opens up shown in Figure 5(c). In contrast with the direct excitation path, this indirect excitation path contains an additional factor and evolves into a series of even powers in E(t). Then, even harmonic orders will be obtained by Fourier transform of the even powers of E(t). Moreover, quantum interference of the direct and indirect excitation paths shown in Figure 5(d) was used to control the harmonic emission in the time domain directly [6,21].

Figure 5. (a) Simplified band structure schematics consisting of two valence bands h₁ and h₂ (red lines) and one conduction band e (blue line). The allowed transitions are labeled and marked by black arrows. (b) Schematic of the direct excitation path. The transition h₁ to e is marked by the blue arrow. (c) Indirect excitation path consisting of the transition h₁ to h₂ (red arrow) and h₂ to e (three-dimensional red arrow). (d) Full quantum interference (QI): combination of the direct excitation path and the indirect excitation path. Adapted from Ref. [6].

Similar results can be found in the orientation-dependent high harmonic spectrum by multi-band SBE model [48]. In Ref [48], three valence bands and five conduction bands are used in SBE simulation. Figure 6(a) is the simulated crystal orientation-dependent high harmonic spectra from h-BN. Significantly, the generation of even harmonics can be found when laser field is polarized around high symmetric direction Γ–M. However, even harmonics disappear clearly if one selectively turns off the indirect excitation pathways as shown in Figure 6(b). Therefore, quantum interference of crystal electrons in the direct and indirect excitation pathways from strong coupling between the valence bands or conduction bands can be regarded as the physical mechanism of even harmonics generation.

Figure 6. (a) Simulated orientation-dependent high harmonic spectra from h-BN by using three valence bands and five conduction bands. The peak intensity of laser is 5.0 × 10¹² W/cm², wavelength is 800 nm and pulse duration is 10 optical cycles. (b) Identical to (a) but by selectively turning off the indirect excitation pathways. (c) High harmonic spectra along the high-symmetry directions Γ–M and Γ–K in the Brillouin zone. Adapted from Ref [48].

Furthermore, two main plateaus were found in the orientation-dependent high harmonic spectra from h-BN, which exhibit different anisotropic distributions as shown in Figure 6(a). Figure 6(c) shows the details of double plateau structure in Γ–M and Γ–K directions. The appearance of two plateaus indicates that strong interband couplings are involved in multiple bands [65], and the couplings of higher conduction band with other valence bands dominate the harmonic emission of second plateau.

4. Numerical method based on time-dependent density-functional theory

By employing an ab initio approach based on TDDFT [102,103], researchers have studied the microscopic origin of HHG [30–32,34,81–84] and electrical currents [15,79,80] in solids. The time-dependent Kohn-Sham (KS) equation reads

(30) $i{\mathrm{\partial }}_t{\psi }_i(\mathbf{r},t)=H(\mathbf{r},t){\psi }_i(\mathbf{r},t),$(30)

where the index i runs over the occupied KS orbitals ψ_i with the Hamiltonian:

(31) $H(\mathbf{r},t)=\frac{1}{2}{\left[-i\mathrm{\nabla }+\mathbf{A}(t)\right]}^2+{\hat{V}}_{ion}+\int d{\mathbf{r}}^{\prime }\frac{n({\mathbf{r}}^{\prime },t)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}+{\hat{V}}_{XC}(\mathbf{r},t),$(31)

where A(t) is vector potential, ${\hat{V}}_{ion}$ is periodic lattice potential, and ${\hat{V}}_{XC}$ is the exchange and correlation potential which plays a key role in describing relaxation and dissipation in an interacting many-electron system [104,105]. The valence electron density is given as $n(\mathbf{r},t)=\sum _i{\left|{\psi }_i(\mathbf{r},t)\right|}^2$. After solving KS equation, time- and space-dependent microscopic current density is received as

(32) $\mathbf{j}(\mathbf{r},t)=\left|e\right|\sum _i\frac{1}{2}\left[{\psi }_i^{\ast }(\mathbf{r},t)(-i\mathrm{\nabla }+\mathbf{A}(t)){\psi }_i(\mathbf{r},t)+c.c.\right].$(32)

Then, the macroscopic current density J(t) along the laser polarization direction F₀ is given by the average of j(r, t) over the unit cell with volume Ω:

(33) $J(t)=\frac{1}{\mathrm{\Omega }}{\int }_{\mathrm{\Omega }}d\mathbf{r}\mathbf{j}(\mathbf{r},t).{\mathbf{F}}_{\mathbf{0}}/\left|F_0\right|.$(33)

Taking advantage of TDDFT model, Tancogne-Dejean et al. [34] found that the HHG spectrum of bulk silicon does not change if consider either the full evolution of the Hartree and the exchange-correlation parts of the KS Hamiltonian or the time evolution in a static ground-state potential in Figure 7(a), and the discovery justifies the independent particle approximation assumed in most previously published HHG models. Here, the pulse duration of laser is 25 fs, with a sin² envelope of the vector potential. The peak intensity inside matter is taken to be I₀ = 10¹¹ W/cm², and the wavelength is 3000 nm [34]. Besides, they investigated the effect of the laser polarization on the HHG emission. In Figure 7(b), the TDDFT simulations clearly display an anisotropic emission of harmonics while rotating the polarization around the [001] crystallographic axis.

Figure 7. (a) Simulated HHG spectra from bulk silicon, for polarization along Γ-X, computed within the LDA (LDA, black line) and within the LDA, but freezing the coulomb and exchange-correlation terms to their ground-state value (LDA-FreezeHXC, red line). (b) The anisotropic HHG spectra calculated from TDDFT model by rotating the laser polarization around the [001] crystallographic direction. Adapted from Ref. [34].

Moreover, on basis of TDDFT method, HHG was studied in one-dimensional structures of sizes from single nucleus up to hundreds of nuclei, and the dependence of system size on the observed HHG cutoffs is found [31]. Bauer et al. [30] studied the two topological phases of a finite, one-dimensional, periodic structure, and astonishingly, harmonic yield is found to differ by up to 14 orders of magnitude for the two topological phases. Tancogne-Dejean et al. [32] demonstrated the importance of a dynamically modulated Hubbard U with the ACBN0 functional in the description of the high-harmonic generation of NiO. Combining TDDFT model and multi-band SBEs model, ab initio multiscale simulations of solid-state HHG were also performed by Floss et al. [82], they found strong propagation and field inhomogeneity effects on the harmonic spectra in the dense medium. The propagation effects can be considered by combining the microscopic SBE solution for the current density J(t) of individual cells with the solution of Maxwell’s equations:

(34) $\frac{1}{c^2}{\mathrm{\partial }}_t^2\vec{A}(\vec{r},t)-{\mathrm{\nabla }}^2\vec{A}(\vec{r},t)=\frac{4\pi }{c}\vec{J}(\vec{r},t)$(34)

where the cells are placed on a mesoscopic grid along the propagation direction.

5. Conclusion and outlook

In conclusion, we introduce recent progress of HHG in solid materials, and summarize the development of three primary numerical methods including TDSE, SBE, and TDDFT models. Based on Bloch state basis and Houston state basis, TDSE model is able to investigate the electronic dynamics in periodic potential. However, it is difficult to simulate the HHG of real solid materials as the method starting from the model potential. The SBE could be improved from two-band model to multi-band model to study the real systems and capture the main features of HHG in solids if we use the accurate energy bands and transition dipole moments from first-principles calculations. The future efforts should be made to obtain the correct phase of transition dipole elements. One big challenge still exists because current first-principles codes always give random phases for transition dipole moments. TDDFT model seems as the optimal method to study the HHG in solid straightforwardly in real coordinate space, its problem may be huge time-consuming in computation, and sometimes it is hard to directly and intuitively analyze the physical process or mechanism within solid-state energy band picture. It should be noted that the propagation effects in the dense medium should be considered in all these methods as it does affect HHG in solid [82]. Also, how to properly incorporate the Berry curvature [45] in models and what role in solid HHG need more attention in the future theoretical simulations. Although several open questions still remain in this emerging field, it is hoped that this review can provide valuable reference and continuous endeavors from the relevant community will gradually lift the veil on solid HHG studies.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by National Natural Science Foundation of China (Grant No. 11704187) and the Natural Science Foundation of Jiangsu Province (Grant No. BK20170032).