Introduction to theory of high-harmonic generation in solids: tutorial

Lun Yue; Lun Yue; Mette B. Gaarde; Mette B. Gaarde

doi:10.1364/JOSAB.448602

1. INTRODUCTION

High-harmonic generation (HHG) is an extremely nonlinear optical process where a macroscopic system irradiated by intense laser light emits coherent radiation with frequencies many times that of the driving laser field. HHG in gases has facilitated the generation of attosecond light pulses, the probing and control of electrons on their natural timescales, and more generally laid the foundation for attosecond science [1,2]. The last decade has seen HHG extended to systems in the condensed phase and liquid phase [3–6]. Earlier measurements of HHG in reflection [7–10] and transmission geometries [11] have not explicitly established the non-perturbative nature of HHG. Since the initial observation of nonperturbative HHG by Ghimire et al. [12] in a bulk solid in 2011, HHG has been observed in semiconductors [13–16], dielectrics [17,18], rare-gas solids [19], monolayer materials [20–23], nanostructures [24–26], amorphous solids [27], doped systems [28], and topological insulators [29,30]. Solid-state HHG has exciting prospects for new compact attosecond light-source technologies [17,18,24–26], as well as novel ultrafast spectroscopy methods capable of probing band structures [31,32], Berry curvatures [21,33], and topological effects [29,30,34].

An intuitive semiclassical understanding of HHG in condensed phase systems as a three-step process is illustrated in Fig. 1, in both reciprocal space (top panels) and real space (bottom panels). In reciprocal space, the first step is the creation of an electron-hole pair via excitation of an electron from the valence band to the conduction band, usually by tunneling near the minimum band gap. Per the acceleration theorem [35,36], the time-dependent crystal momentum will follow the time-dependent vector potential, and the resulting carrier motion in the (nonparabolic) bands will lead to the emission of nonperturbative intraband harmonics [12,37]. At the same time, also interband harmonic radiation is emitted via recombination, with frequencies corresponding to the instantaneous band gap [14,38]. The interband radiation results from the coherence of the electron-hole pair, and is emitted when a stationary phase condition is satisfied for the phase that is accumulated during propagation in the bands. In real space, the three steps again consist of tunneling, which creates the electron-hole pair; propagation, which accelerates them apart in space and leads to intraband emission; and recollision, when the electron and hole reencounter each other in space and drive interband emission via recombination. The recollision corresponds to the stationary phase condition (from the reciprocal picture) being exactly satisfied. The real-space semiclassical model represents the generalization of the recollision model for gas-phase HHG [39,40] to the condensed phase. Although the picture illustrated in Fig. 1 is based on a semiclassical understanding of the HHG process, its interpretation has been supported by time-dependent strong-field light–matter simulations such as the single-particle time-dependent Schrödinger equation (TDSE) [41–44], time-dependent density-functional theory (TD-DFT) [45–51], SBEs [52–54], and density matrix approaches [38,55,56].

Fig. 1. Sketch of the recollision model for HHG. Upper panels: Tunneling, propagation, and recombination in reciprocal space. Lower panels: Equivalent physics in real space.

Download Full Size | PDF

Any calculation of HHG from a condensed-phase material interacting with a strong laser field involves at least three different types of calculations: Finding the initial state of the material, solving an equation of motion (EOM) to describe the time evolution, and calculating the observable current(s) that will make up the harmonic spectrum. For a periodic material, these calculations are most conveniently performed in reciprocal space, using Bloch states, so that both the structure and the dynamics of the system are described in terms of its band structure. As we will describe in more detail in the next couple of paragraphs, these calculations involve: (i) Calculating the band structure, including relevant matrix elements of the momentum or position operators, by solving the time-independent Schrödinger equation (TISE); (ii) Solving the EOM after choosing a structure as well as a laser gauge to calculate the dynamics, where the gauge choice will determine the exact form of the EOM; and (iii) Calculating the total current and perhaps different contributions to the current that have different physical meaning, such as the intraband and interband currents described above, or the anomalous current driven by the Berry curvature [21,33].

Starting with the structure calculation, there is an important, and sometimes subtle, complication that deserves special attention in the theoretical description of carrier dynamics in a solid. Since the different crystal momenta ${\textbf{k}}$ are treated independently in the structure calculation in reciprocal space, there is a ${\textbf{k}}$-dependent phase arbitrariness in the Bloch states, i.e., a gauge freedom, which we will henceforth refer to as the structure gauge. It is often favorable to pick a gauge (generally referred to as gauge fixing) where the Bloch functions are smooth, especially if one considers the dynamics in terms of electrons and holes explicitly moving along the bands. In the case of large carrier excursions in the Brillouin zone (BZ), as is often the case for carrier motion beyond the perturbative regime and the breakdown of the notion of effective masses, a BZ periodic gauge is also required. The construction of such a smooth and BZ periodic structure gauge is intricate, especially in the case of bands containing degeneracies [57–59]. As we will discuss in more detail in Section 2, the structure gauge is closely related to the Berry phase [60] that is accumulated under adiabatic motion, which in condensed matter theory has had profound implications such as the development of the modern theory of polarization [61] and the discovery of the quantization of adiabatic transport [62].

In the description of the time-dependent interaction between the crystal and the strong laser field there is another gauge freedom; namely, the choice of the laser gauge. While all physical observables in principle are gauge-invariant in a complete basis, for computational purposes the number of bands and number of ${\textbf{k}}$ points in the BZ should be truncated. This basis truncation critically depends on the choice of the laser gauge [56,63–66]. In addition, the choice of the laser gauge is linked to the fixing of the structure gauge [56], since the band-coupling terms are different in different laser gauges [67].

Finally, while the microscopic theory of HHG can be solved in the framework of the dipole approximation, the experimentally measurable signal is actually the macroscopic HHG response. Hence, a realistic treatment of HHG should involve the effect of beam propagation in the medium, as well as the propagation of the radiation from the near-field at the sample to the far-field at the detector. These effects requires the solution of the coupled microscopic response to Maxwell’s equations, which is very computationally demanding [55,68–72].

This tutorial provides a hands-on introduction to doing calculations of strong-field laser–matter interactions in the condensed phase, with an emphasis on the generation of high-order harmonics. Although the tutorial is primarily aimed at scientists specifically interested in performing calculations of ultrafast condensed-phase dynamics, we believe that many of the concepts and descriptions will be useful to anyone wishing to gain a deeper understanding of the current state and capabilities of solid-state HHG theory. Indeed, the decomposition of the current, the semiclassical recollision model, and the macroscopic effects can directly be used to interpret experiments. It is also worth noting that HHG in solids is an inherently interdisciplinary scientific field, attracting researchers from both the condensed-matter physics community and the strong-field/attoscience community, which originated with atomic, molecular, and optical physicists. We thus also aim to provide some unification of concepts across this diverse community of scientists.

We start the tutorial in Section 2 with an introduction to adiabatic states and related concepts such as the Berry phase, connections, and curvatures, all of which will be useful for reading the later sections. Section 3 deals with the time-independent problem of a crystalline solid, and provides methods to construct a smooth and BZ periodic structure gauge. The random gauge, the parallel transport (PT) gauge, the twisted PT (TPT) gauge and the Wannier gauge, as well as degenerate bands, will be covered. Section 4 treats the time-dependent microscopic problem of HHG, and presents the relevant EOMs in the velocity gauge (VG) and the length gauge (LG), as well as in the time-dependent adiabatic Houston basis. The advantages and drawbacks of the different methods are compared, and we present a concrete calculation example for HHG in a monolayer material. Section 5 is about the saddle-point method and the semiclassical solutions to the saddle point equations, i.e., the recollision model. Section 6 gives a brief introduction to the macroscopic propagation schemes for HHG and provides an example of a near-field to far-field propagation scheme and discusses the spatio-spectral properties of the far-field spectrum.

Atomic units, where the reduced Planck constant, the elementary charge, the Bohr radius and the electron mass are set to unity, are used throughout this work unless indicated otherwise.

2. ADIABATIC STATES, BERRY CONNECTIONS, AND CURVATURES

We start this tutorial by briefly introducing the adiabatic states and some general concepts that could be helpful to understand the rest of the tutorial. Concepts such as Berry phases, Berry connections, Berry curvatures, and gauge fixing will be discussed. For further reading, see [57,58,60,73–75].

A. Adiabatic States

When describing laser–matter interactions, the Hamiltonian $\hat H(t)$ is generally time-dependent, and it is often useful to describe the dynamics of the system using the adiabatic states. These are defined as the eigenstates of the instantaneous Hamiltonian at time $t$:

(1)$$\hat H(t)|n(t)\rangle = {\epsilon _n}(t)|n(t)\rangle ,$$

with $|n(t)\rangle$ an adiabatic state and ${\epsilon _n}(t)$ its energy. The TDSE reads

(2)$$i|\dot \Psi (t)\rangle = \hat H(t)|\Psi (t)\rangle ,$$

with $|\Psi (t)\rangle$ the quantum state, and the diacritic dot is used henceforth to denote full time derivatives. It is important to note that an instantaneous eigenstate satisfying Eq. (2) at a given $t$ is not uniquely defined, but rather carries an arbitrary phase factor; i.e., there is a gauge freedom. This gauge freedom and its consequences will recur many times throughout this document.

We consider the case where the adiabatic states are nondegenerate for all $t$. The degenerate case will be discussed for the Brillouin problem in Section 3.F. The wave function expanded into the adiabatic states reads

(3)$$|\Psi (t)\rangle = \sum\limits_n {c_{\!n}}(t){e^{i{\theta _n}(t)}}|n(t)\rangle ,$$

with ${\theta _n}(t) = - \int {\epsilon _n}(t^\prime){\rm d}t^\prime $ as the dynamical phase. We can get the EOM for the time-dependent coefficients ${c_{\!n}}(t)$ by inserting Eq. (3) into the TDSE of Eq. (2) and projecting onto $\langle n(t)|$,

(4)$${\dot c_n}(t) = - \langle n(t)|\dot n(t)\rangle {c_{\!n}}(t) - \sum\limits_{m \ne n} {c_m}(t)\langle n(t)|\dot m(t)\rangle {e^{i[{\theta _m}(t) - {\theta _n}(t)]}}.$$

The first term on the right-hand side of Eq. (4) depends only on the $n$th state, while the second term involves the nonadiabatic couplings $\langle n(t)|\dot m(t)\rangle$ responsible for the transfer of population between adiabatic states. Rewriting $\langle n(t)|\dot m(t)\rangle = - \langle n(t)|{\dot {\hat H}}(t)|m(t)\rangle /[{{\epsilon _n}(t) - {\epsilon _m}(t)}]$ for $m \ne n$, one can see that the nonadiabatic couplings are generally large when the adiabatic states are close to each other in energy.

In the adiabatic approximation, the terms in Eq. (4) involving the nonadiabatic couplings are neglected. For a system that starts in the $n$th adiabatic state, $\Psi ({t_0}) = |n({t_0})\rangle$, such that ${c_{\!n}}({t_0}) = {\delta _{\textit{mn}}}$, the time evolution proceeds as

(5)$$|\Psi (t)\rangle = {e^{i{\theta _n}(t)}}{e^{i{\gamma _n}(t)}}|n(t)\rangle ,$$

where

(6)$${\gamma _n}(t) = i\int_{{t_0}}^t {\rm d}t^\prime \langle n(t^\prime)|\dot n(t^\prime)\rangle$$

is the Berry phase or the geometric phase [57,60].

B. Berry Phase, Connection, and Curvature

The Berry phase ${\gamma _n}(t)$ in Eq. (6) emerged in the adiabatic state basis as the additional phase accumulated by a state beyond its dynamical phase ${\theta _n}$. The Berry phase is often termed the geometrical phase because of its geometric properties. These can be realized by assuming without loss of generality that the Hamiltonian depends on time through a set of parameters ${\textbf{R}}(t) = [{R_1}(t),{R_2}(t), \ldots]$, such that $\hat H(t) = \hat H[{\textbf{R}}(t)]$. The Berry phase can now be rewritten as an integral over a path ${\cal C}$ through this parameter space, from ${\textbf{R}}({t_0})$ to ${\textbf{R}}(t)$:

(7)$${\gamma _n}({\cal C}) = \int_{\cal C} {{\cal A}_n}({\textbf{R}}) \cdot {\rm d}{\textbf{R}},$$

and the purely real Berry connection is defined as

(8)$${{\cal A}_n}({\textbf{R}}) \equiv i\langle n({\textbf{R}})|{\partial _{\textbf{R}}}|n({\textbf{R}})\rangle ,$$

with ${\partial _{\textbf{R}}} = [{\partial _{{R_1}}},{\partial _{{R_2}}}, \ldots]$ the operator for the partial derivatives. Note that the Berry phase in this formulation only depends on the path and no longer on the time duration.

For a closed path, the Berry phase can be rewritten to a surface integral by the application of Stoke’s theorem, with ${\cal S}({\cal C})$ denoting a surface enclosed by ${\cal C}$. For convenience, we assume that the parameter space ${\textbf{R}}$ is 3D, such that ${\partial _{\textbf{R}}} = {\nabla _{\textbf{R}}}$, although derivations for higher dimensions proceed similarly. We have

(9)$${\gamma _n}({\cal C}) = \oint_{\cal C} {{\cal A}_n}({\textbf{R}}) \cdot {\rm d}{\textbf{R}} = \iint_{{\cal S}({\cal C})} \left[{{\nabla _{\textbf{R}}} \times {{\cal A}_n}({\textbf{R}})} \right] \cdot {\rm d}{\textbf{S}},$$

where the integrand in the surface integral is the Berry curvature

(10)$${\boldsymbol{\Omega} _n}({\textbf{R}}) \equiv {\nabla _{\textbf{R}}} \times {{\cal A}_n}({\textbf{R}}).$$

The Berry curvature is important because, like the Berry phase, it is independent under a gauge transform. This is in contrast to the Berry connection, which transforms as ${{\cal A}_n} \to {{\cal A}_n} + {\nabla _{\textbf{R}}}\chi$ under a gauge transform of the adiabatic states $|n\rangle \to |n\rangle {e^{- i\chi ({\textbf{R}})}}$, with $\chi$ an arbitrary real function. Note that the Berry connection and the Berry curvature transform similarly to the vector potential and magnetic field from classical electromagnetism, respectively.

The Berry curvature can be expressed in a form that is useful for numerical evaluations [76]. Using $\nabla \times (\psi \nabla \phi) = \nabla \psi \times \nabla \phi$, the resolution of identity, and $\langle m({\textbf{R}})|{\nabla _{\textbf{R}}}n({\textbf{R}})\rangle [{\epsilon _n}({\textbf{R}}) - {\epsilon _m}({\textbf{R}})] = \langle m({\textbf{R}})|{\nabla _{\textbf{R}}}\hat H({\textbf{R}})|n({\textbf{R}})\rangle$ for $m \ne n$ obtained by taking the gradient of Eq. (2), the Berry curvature can be written

(11)$$\begin{split}{\boldsymbol{\Omega} _n}({\textbf{R}})&= i\langle {\nabla _{\textbf{R}}}n({\textbf{R)}}| \times |{\nabla _{\textbf{R}}}n({\textbf{R)}}\rangle\\&= i\sum\limits_m \langle {\nabla _{\textbf{R}}}n({\textbf{R)}}|m({\textbf{R}})\rangle \times \langle m({\textbf{R}})|{\nabla _{\textbf{R}}}n({\textbf{R)}}\rangle\\&= i\sum\limits_{m \ne n} \frac{{\langle n({\textbf{R)}}|{\nabla _{\textbf{R}}}\hat H({\textbf{R}})|m({\textbf{R}})\rangle \times \langle m({\textbf{R}})|{\nabla _{\textbf{R}}}\hat H({\textbf{R}})|n({\textbf{R)}}\rangle}}{{{{\left[{{\epsilon _m}({\textbf{R}}) - {\epsilon _n}({\textbf{R}})} \right]}^2}}}.\end{split}$$

From this expression, the Berry curvature for the $n$th adiabatic state is intuitively seen to originate from the coupling to all the other adiabatic states, and the closer in energy a neighboring adiabatic state is, the larger the contribution. Equation (11) is useful for numerical evaluations since no derivatives are taken with respect to the adiabatic states.

As mentioned above, there is a gauge freedom in the adiabatic states, i.e., a phase-arbitrariness in the definition of the adiabatic states. One way to remove this arbitrariness of the phase is to enforce the parallel transport (PT) gauge condition on the eigenstates:

(12)$$0 = \langle n(t)|\dot n(t)\rangle = \dot{\textbf{R}} \cdot \langle n({\textbf{R}})|{\partial _{\textbf{R}}}|n({\textbf{R}})\rangle .$$

In the PT gauge, the “velocity” of the state is thus perpendicular to the state itself. With this gauge choice, the eigenstate $|n(t)\rangle$ is always single-valued as a function of $t$, and the integrand in Eq. (7) is zero along the path. However, if ${\cal C}$ is a closed path in parameter space (e.g., under periodic motion) such that ${\textbf{R}}({t_0}) = {\textbf{R}}(T)$, there is no guarantee that $|n[{\textbf{R}}({t_0})]\rangle$ is equal to $|n[{\textbf{R}}(T)]\rangle$ under PT. The phase difference between $|n[{\textbf{R}}({t_0})]\rangle$ and $|n[{\textbf{R}}(T)]\rangle$ is exactly the Berry phase ${\gamma _n}(T)$. If ${\gamma _n}(T) \ne 0$, the adiabatic state $|n\rangle$ in PT gauge is thus a multivalued function of ${\textbf{R}}$ (but still a single-valued function of $t$). Instead of the PT gauge, a sometimes more convenient gauge is the periodic gauge, in which one applies the additional requirement that the adiabatic state is single-valued along the path ${\cal C}$ in parameter space. Later, in Section 3 of this tutorial we will present methods to construct such a gauge in crystalline systems.

Finally, by neglecting adiabatic coupling in the second term on the right hand side of Eq. (4), the quantum evolution of $|\Psi (t)\rangle$ in Eq. (5) is correct only to the zeroth order in $\dot{\textbf{ R}}$; i.e., if it starts in the $n$th adiabatic state, it stays in that state, with no transition to other adiabatic states. The result can be extended to the first order by making the ansatz $|{\Psi ^{(1)}}\rangle = {e^{i{\theta _n}}}{e^{i{\gamma _n}}}[{|n\rangle + \dot{\textbf{R}}|\delta n\rangle}]$ in the TDSE and solving for $|\delta n\rangle$, resulting in the state [58,74]

(13)$$|{\Psi ^{(1)}}(t)\rangle = {e^{i{\theta _n}(t)}}{e^{i{\gamma _n}(t)}}\left[{|n(t)\rangle + i\sum\limits_{m \ne n} \frac{{\langle m(t)|\dot n(t)\rangle}}{{{\epsilon _m}(t) - {\epsilon _n}(t)}}|m(t)\rangle} \right].$$

The second term on the right-hand side of Eq. (13) results in the so-called anomalous current, which we will return to in Section 4. This expression, together with the Berry phase, has been crucial for the development of the theory of adiabatic charge transport [62] and the modern theory of polarization [57,61].

3. TIME-INDEPENDENT PROBLEM AND STRUCTURE GAUGES

In this section, we discuss the time-independent problem of a crystalline solid. While for a perfect crystal it can be argued that everything boils down to the TISE and Bloch’s theorem, the Bloch states can have an arbitrary crystal-momentum-dependent phase, representing a gauge degree of freedom. This section is devoted to the discussion of this gauge freedom and presents methods for gauge fixing. We start with Bloch’s theorem, the choice of crystal conventions, and how to calculate relevant matrix elements. We then discuss four different gauge choices in Sections 3.B–3.E and conclude with Section 3.F for gauge fixing in the nondegenerate case.

A. Bloch States and Coupling Matrix Elements

The single-particle eigenstates satisfying Bloch’s theorem are solutions to the crystalline TISE

(14)$${\hat H_0}\phi _n^{\textbf{k}}({\textbf{r}}) = E_n^{\textbf{k}}\phi _n^{\textbf{k}}({\textbf{r}}),$$

where the Hamiltonian ${\hat H_0} = \hat T + \hat V$ contains the kinetic energy operator $\hat T = \hat{{\textbf{p}}}^2/2$ and the periodic potential satisfying $V({\textbf{r}}) = V({\textbf{r}} + {\textbf{R}})$ with ${\textbf{R}}$ as a lattice vector. The band energies $E_n^{\textbf{k}}$ and Bloch states $\phi _n^{\textbf{k}}({\textbf{r}})$ in Eq. (14) depend on the band index $n$ and the crystal momentum ${\textbf{k}}$. According to Bloch’s theorem, the Bloch states can be written $\phi _n^{\textbf{k}}({\textbf{r}}) = {e^{i{\textbf{k}} \cdot {\textbf{r}}}}u_n^{\textbf{k}}({\textbf{r}})$, with $u_n^{\textbf{k}}({\textbf{r}})$ as a lattice cell periodic function. Although in this tutorial we focus on systems without electron–electron interactions, it has been shown that most of the conclusions can be extended to interacting systems, provided that the interacting system is adiabatically connected to a noninteracting one. Two systems are adiabatically connected if one can go to the other by following a suitable path in parameter space, i.e., in this case that the interacting system can be reduced to the noninteracting one by an adiabatic path that takes the finite interaction strength to zero [58,74]. For recent studies on electron-correlation effects in strong-field laser-solid interactions using TD-DFT, refer to [46,47,49–51].

Before proceeding, some details on the choice of lattice are warranted. We consider a $D$ dimensional crystal with the Bravais lattice

(15)$${\textbf{R}} = \sum\limits_{d = 1}^D {m_d}{{\textbf{a}}_d},\quad \quad {m_d} \in \left[{- \frac{{{N_d}}}{2},\frac{{{N_d}}}{2} - 1} \right],$$

that has volume ${V_{{\rm{crys}}}} = N{V_{{\rm{cell}}}}$, where $N = \prod\nolimits_{d = 1}^D {N_d}$ is the total number of unit cells, ${V_{{\rm{cell}}}}$ is the unit cell volume, and $D$ the lattice dimension. To avoid surface effects, periodic boundary conditions $\phi _n^{\textbf{k}}({\textbf{r}} + {N_d}{{\textbf{a}}_d}) = \phi _n^{\textbf{k}}({\textbf{r}})$ are enforced [53], which results in discretized crystal momenta in a BZ

(16)$$\!\!\!{\textbf{k}} = \sum\limits_{d = 1}^D \frac{{{n_d}}}{{{N_d}}}{{\textbf{b}}_d} = \sum\limits_{d = 1}^D {\kappa _d}\hat{{\textbf{b}}}_d,\;\; {n_d} \in \left[{- \frac{{{N_d}}}{2},\frac{{{N_d}}}{2} - 1} \right],\!$$

where we have defined the reduced reciprocal coordinates $2{\kappa _d} \in [- \left\| {{{\textbf{b}}_d}} \right\|,\left\| {{{\textbf{b}}_d}} \right\|]$ along the reciprocal vector directions with grid spacings $\Delta {\kappa _d} = \left\| {{{\textbf{b}}_d}} \right\|/{N_d}$, and $\hat{{\textbf{b}}}_d = {{\textbf{b}}_d}/\left\| {{{\textbf{b}}_d}} \right\|$ is an unit vector. Two useful orthogonality relations are

(17)$$\sum\limits_{{\textbf{k}} \in BZ} {e^{i{\textbf{k}} \cdot ({\textbf{R}} - {\textbf{R}}^\prime)}} = N{\delta _{{\textbf{R}},{\textbf{R}}^\prime}},$$

(18)$$\sum\limits_{\textbf{R}} {e^{- i({\textbf{q}} - {\textbf{k}}) \cdot {\textbf{R}}}} = N{\delta _{{\textbf{q}},{\textbf{k}}}},$$

and we define the crystal and unit cell inner products as

(19)$${\langle f|g\rangle _{{\rm{crys}}}} \equiv \int_{{\rm{crys}}} {f^*}({\textbf{r}})g({\textbf{r}}){\rm d}{\textbf{r}},$$

(20)$${\langle f|g\rangle _{{\rm{cell}}}} \equiv \int_{{\rm{cell}}} {f^*}({\textbf{r}})g({\textbf{r}}){\rm d}{\textbf{r}}.$$

As we will see in Section 4, the matrix elements of interest to solve the time-dependent problem are those of the momentum and position operators, $\hat{\textbf{p}}$ and $\hat{\textbf{r}}$, respectively. To provide an example of the steps involved in the evaluation of matrix elements, it can be seen that operators $\hat{O}(\hat{\textbf{p}})$ involving the momentum operator $\hat{\textbf{p}}$ are diagonal in the crystal momentum:

(21)$$\begin{split}&{{\langle \phi _m^{\textbf{q}}|\hat{O}(\hat{\textbf{p}})|\phi _n^{\textbf{k}}\rangle}_{{\rm{crys}}}}\\& = \int_{{\rm{crys}}} \phi _m^{{\textbf{q}}*}({\textbf{r}})\hat{O}(\hat{\textbf{p}})\phi _n^{\textbf{k}}({\textbf{r}}){\rm d}{\textbf{r}}\\ &= {\sum\limits_{\textbf{R}} \int_{{\rm{cell}}} u_m^{{\textbf{q}}*}({\textbf{r}}){e^{- i{\textbf{q}} \cdot ({\textbf{r}} + {\textbf{R}})}}\hat{O}(\hat{\textbf{p}})\left[{u_n^{\textbf{k}}({\textbf{r}}){e^{i{\textbf{k}} \cdot ({\textbf{r}} + {\textbf{R}})}}} \right]{\rm d}{\textbf{r}}}\\ &= {\sum\limits_{\textbf{R}} {e^{- i({\textbf{q}} - {\textbf{k}}) \cdot {\textbf{R}}}}\int_{{\rm{cell}}} u_m^{{\textbf{k}^\prime}*}({\textbf{r}}){e^{- i{\textbf{q}} \cdot {\textbf{r}}}}\hat{O}(\hat{\textbf{p}})\left[{u_n^{\textbf{k}}({\textbf{r}}){e^{i{\textbf{k}} \cdot {\textbf{r}}}}} \right]{\rm d}{\textbf{r}}}\\ &= {N{\delta _{{\textbf{q}},{\textbf{k}}}}\int_{{\rm{cell}}} \phi _m^{{\textbf{q}}*}({\textbf{r}})\hat{O}(\hat{\textbf{p}})\phi _n^{\textbf{k}}({\textbf{r}}){\rm d}{\textbf{r}}}\\& = N{\delta _{{\textbf{q}},{\textbf{k}}}}{{\langle \phi _m^{\textbf{q}}|\hat{O}(\hat{\textbf{p}})|\phi _n^{\textbf{k}}\rangle}_{{\rm{cell}}}},\end{split}$$

where we have used Eq. (18). For $\hat{O} \equiv \hat 1$ in Eq. (21),

(22)$${\langle \phi _m^{\textbf{q}}|\phi _n^{\textbf{k}}\rangle _{{\rm{crys}}}} = N{\delta _{{\textbf{q}},{\textbf{k}}}}{\langle u_m^{\textbf{q}}|u_n^{\textbf{k}}\rangle _{{\rm{cell}}}} \equiv N{\delta _{{\textbf{q}},{\textbf{k}}}}{\delta _{\textit{mn}}},$$

where the last equality represents the choice of orthonormality for the cell-periodic functions. Thus, in our convention, the resolution of identity is

(23)$$1 = {N^{- 1}}\sum\limits_{n{\textbf{k}}} |\phi _n^{\textbf{k}}\rangle \langle \phi _n^{\textbf{k}}{|_{{\rm{crys}}}} = \sum\limits_{n{\textbf{k}}} |u_n^{\textbf{k}}\rangle \langle u_n^{\textbf{k}}{|_{{\rm{cell}}}}.$$

For $\hat{O} \equiv \hat{\textbf{p}}$ in Eq. (21), we have

(24)$${\langle \phi _m^{\textbf{q}}|\hat{\textbf{p}}|\phi _n^{\textbf{k}}\rangle _{{\rm{crys}}}} = N{\delta _{{\textbf{q}},{\textbf{k}}}}{\textbf{p}}_{\textit{mn}}^{\textbf{k}},$$

where we have defined the momentum matrix elements ${\textbf{p}}_{\textit{mn}}^{\textbf{k}} = {\langle \phi _m^{\textbf{q}}|\hat{\textbf{p}}|\phi _n^{\textbf{k}}\rangle _{{\rm{cell}}}}$.

Similarly, to calculate the matrix elements of the position operator $\hat{\textbf{r}}$ in the Bloch basis, one would go through similar steps as in Eq. (21). However, the third equality would have an extra term involving the factor $\sum\nolimits_{\textbf{R}} {\textbf{R}}{e^{- i({\textbf{q}} - {\textbf{k}}) \cdot {\textbf{R}}}}$, which is difficult to evaluate and depends on the choice of the Bravais lattice. To circumvent this, we first consider the continuum limit where the supercell volume goes to infinity, such that ${\textbf{k}}$ sums turn into integrals, and the Kronecker delta into delta functions,

(25)$${N^{- 1}}\sum\limits_{{\textbf{k}} \in BZ} \to \frac{{{V_{{\rm{cell}}}}}}{{{{(2\pi)}^D}}}\int_{\textit{BZ}} {\rm d}{\textbf{k}},$$

(26)$$N{\delta _{{\textbf{q}},{\textbf{k}}}} \to \frac{{{{(2\pi)}^D}}}{{{V_{{\rm{cell}}}}}}\delta ({\textbf{q}} - {\textbf{k}}).$$

Hence, we mostly use the discrete notation. However, when ${\textbf{k}}$ derivatives are involved, the continuous case in Eq. (25) is implied. In contrast to the momentum operator in Eq. (21), the position operator couples Bloch states in ${\textbf{k}}$ space [67],

(27)$$\begin{split}{{{\langle \phi _m^{\textbf{q}}|\hat{\textbf{r}}|\phi _n^{\textbf{k}}\rangle}_{{\rm{crys}}}}}&= \int_{{\rm{crys}}} {\rm d}{\textbf{r}}\phi _m^{{\textbf{q}}*}({\textbf{r}}){\textbf{r}}\phi _n^{\textbf{k}}({\textbf{r}})\\ &= \int_{{\rm{crys}}} {\rm d}{\textbf{r}}\left\{{i{\nabla _{\textbf{q}}}\phi _m^{{\textbf{q}}*}({\textbf{r}}) - i{\nabla _{\textbf{q}}}u_m^{{\textbf{q}}*}({\textbf{r}}){e^{- i{\textbf{q}} \cdot {\textbf{r}}}}} \right\}\phi _n^{\textbf{k}}({\textbf{r}})\\ &= i{\nabla _{\textbf{q}}}\int_{{\rm{crys}}} {\rm d}{\textbf{r}}\phi _m^{{\textbf{q}}*}({\textbf{r}})\phi _n^{\textbf{k}}({\textbf{r}}) \\&\quad- i\int_{{\rm{crys}}} {\rm d}{\textbf{r}}{\nabla _{\textbf{q}}}u_m^{{\textbf{q}}*}({\textbf{r}}){e^{- i{\textbf{q}} \cdot {\textbf{r}}}}\phi _n^{\textbf{k}}({\textbf{r}})\\ &= iN{\delta _{m,n}}{\nabla _{\textbf{q}}}{\delta _{{\textbf{q}},{\textbf{k}}}} - i\sum\limits_{\textbf{R}} {e^{- ij({\textbf{q}} - {\textbf{k}}) \cdot {\textbf{R}}}}\\&\quad\times\int_{{\rm{cell}}} {\rm d}{\textbf{r}}{e^{- i({\textbf{q}} - {\textbf{k}}) \cdot {\textbf{r}}}}{\nabla _{\textbf{q}}}u_m^{{\textbf{q}}*}({\textbf{r}})u_n^{\textbf{k}}({\textbf{r}})\\ &= N(i{\delta _{m,n}}{\nabla _{\textbf{q}}} + {\textbf{d}}_{\textit{mn}}^{\textbf{k}}){\delta _{{\textbf{q}},{\textbf{k}}}},\end{split}$$

where in the last step we have used Eq. (18), performed partial integration of the second term, and defined the generalized dipole coupling as

(28)$${{\textbf{d}}_{\textit{mn}}^{\textbf{k}} \equiv i{{\langle u_m^{\textbf{k}}|{\nabla _{\textbf{k}}}|u_n^{\textbf{k}}\rangle}_{{\rm{cell}}}}.}$$

Note that our definition does not include the negative charge of the electron. By taking the derivative of the TISE, it can be shown that ${\textbf{d}}_{\textit{mn}}^{\textbf{k}} = - i{\textbf{p}}_{\textit{mn}}^{\textbf{k}}/(E_m^{\textbf{k}} - E_n^{\textbf{k}})$ for the nondiagonal elements. The diagonal elements are the Berry connections [see Eq. (8)],

(29)$${\cal A}_n^{\textbf{k}} \equiv {\textbf{d}}_{\textit{nn}}^{\textbf{k}}.$$

Clearly, to be able to calculate ${\textbf{k}}$ gradients in Eqs. (28) and (29), the cell-periodic functions $|u_n^{\textbf{k}}\rangle$ should be smooth and BZ periodic functions. Generally, ${\textbf{k}}$ gradients are often encountered in the description of dynamics in solids, and it is thus desirable to ensure that $|u_n^{\textbf{k}}\rangle$ fulfills the mentioned properties. This will be the topic of the subsequent subsections.

B. Random Structure Gauge

As discussed above, the Bloch states in the TISE of Eq. (14) are defined up to an arbitrary ${\textbf{k}}$ dependent phase factor ${e^{i\varphi _n^{\textbf{k}}}}$; i.e., there is a gauge degree of freedom. Since this gauge is related to the time-independent problem and the structure of the crystal, we will refer to it as the structure gauge. When Eq. (14) is solved by some diagonalization procedure, the value of $\varphi _n^{\textbf{k}}$ is usually random; we call this the random gauge, and denote the Bloch states in this gauge by $|\breve{\phi} _n^{\textbf{k}}\rangle$ and the corresponding cell-periodic functions by $|\breve{u} _n^{\textbf{k}}\rangle$. The random gauge, while being the simplest structure gauge, is clearly inapplicable in situations when ${\textbf{k}}$ derivatives of the states are needed.

We consider first the nondegenerate case, and will discuss the degenerate case in Section 3.F. The band energies have the BZ periodicity $E_n^{\textbf{k}} = E_n^{{\textbf{k}} + {\textbf{G}}}$, with ${\textbf{G}}$ as a reciprocal lattice vector, and the BZ can be regarded as a torus. With the gauge freedom, an ideal fixed gauge for our purposes is one in which the Bloch functions are BZ periodic because such a gauge would lead to BZ-periodic momentum and dipole matrix elements useful for numerical evaluations. The periodic gauge condition reads $|\phi _n^{{\textbf{k}} + {\textbf{G}}}\rangle = |\phi _n^{\textbf{k}}\rangle$, with ${\textbf{G}} = \sum\nolimits_{d = 1}^D {n_d}{{\textbf{b}}_d}$ an arbitrary reciprocal lattice vector, $D$ the dimension, ${n_d}$ integers, and ${{\textbf{b}}_d}$ the primitive reciprocal lattice vectors. This condition is equivalent to

(30)$$|\phi _n^{{\textbf{k}} + {{\textbf{b}}_d}}\rangle = |\phi _n^{\textbf{k}}\rangle .$$

Here, we outline a procedure to construct such a gauge, by first constructing a PT gauge, and afterward constructing the TPT gauge with the Berry phases along the reciprocal vectors distributed evenly across the BZ.

C. Parallel Transport Structure Gauge

In the PT gauge ($|\bar \phi _n^{\textbf{k}}\rangle$ and $|\bar u_n^{\textbf{k}}\rangle$), similar to what was discussed for Eq. (12) in Section 2, the scalar Berry connections along the reduced coordinates [see Eq. (16)] are forced to vanish,

(31)$$\bar {\cal A}_{n,{\kappa _d}}^{\textbf{k}} \equiv i{\langle \bar u_n^{\textbf{k}}|{\partial _{{\kappa _d}}}|\bar u_n^{\textbf{k}}\rangle _{{\rm{cell}}}} = 0,$$

where the integral is carried over a unit cell. Note that to fix the gauge in Eq. (31), the cell-periodic functions $|u_n^{\textbf{k}}\rangle$ are used over the Bloch states $|\phi _n^{\textbf{k}}\rangle$ to define the scalar Berry connections. Using the latter would have resulted in ambiguous integrals with fast oscillating integrands and dependence on the choice of spatial coordinate origin. When fixing the gauge, any condition can be used; it just turns out that it is more convenient to use the cell-periodic functions. The full Berry connection vector can be constructed from the scalar Berry connections as

(32)$${\cal A}_n^{\textbf{k}} \equiv i{\langle u_n^{\textbf{k}}|{\nabla _{\textbf{k}}}|u_n^{\textbf{k}}\rangle _{{\rm{cell}}}} = \sum\limits_{d,d^\prime}^D {\cal A}_{n,{\kappa _d}}^{\textbf{k}}{g^{d,d^\prime}}\hat{{\textbf{b}}}_{{d^\prime}},$$

with ${g^{d,d^\prime}}$ the inverse metric tensor. In the discrete picture pertinent to numerical evaluations, the Berry connections in the PT gauge can be evaluated as [57]

(33)$$\bar {\cal A}_{n,{\kappa _d}}^{\textbf{k}} = - {\left\| {\delta {\textbf{k}}} \right\|^{- 1}}{\rm{Im}}\ln {\langle \bar u_n^{\textbf{k}}|\bar u_n^{{\textbf{k}} + \delta {\textbf{k}}}\rangle _{{\rm{cell}}}},$$

with $\delta {\textbf{k}}$ a small displacement vector in the reciprocal space.

D. Periodic Structure Gauge

To construct the periodic gauge, note first that the periodic gauge condition for the Bloch states in Eq. (30) translates into the condition

(34)$$|u_n^{{\textbf{k}} + {{\textbf{b}}_d}}\rangle = {e^{- i{{\textbf{b}}_d} \cdot {\textbf{r}}}}|u_n^{\textbf{k}}\rangle$$

for the cell-periodic functions. The $|\bar u_n^{\textbf{k}}\rangle$ constructed in the PT gauge are smooth inside the first BZ, but they generally do not satisfy Eq. (34). For a closed path along ${{\textbf{b}}_d}$ that wraps around the BZ, a Berry phase is accumulated, the so-called Zak’s phase [77],

(35)$$\varphi _{n,{\kappa _d}}^B = \oint_{{\rm{BZ}}} {\bar {\cal A}_{n,{\kappa _d}}}{\rm d}{\kappa _d}.$$

It has been shown by Zak that for a 1D solid with inversion symmetry, this phase can only take on the values 0 or $\pi$.

Consider now the discrete case, and assume that we have constructed the PT gauge along the $\hat{{\textbf{b}}}_d$ direction in discrete steps starting from one end of the BZ at ${{\textbf{k}}_0}$, to the last point ${{\textbf{k}}_{{N_d} - 1}}$. From the PT gauge constraint in Eq. (31), the Berry connection along this path is identically zero. To close the loop of the integral in Eq. (35), we need to wrap around the BZ, and the Berry phase in the PT gauge is then calculated as

(36)$$\varphi _{n,{\kappa _d}}^B = - {\rm{Im}}\ln \langle \bar u_n^{{{\textbf{k}}_{{N_d} - 1}}}|{e^{- i{{\textbf{b}}_d} \cdot {\textbf{r}}}}|\bar u_n^{{{\textbf{k}}_0}}\rangle ,$$

where we have used Eq. (34), and $\Delta {\kappa _d}$ is the grid spacing of the reduced reciprocal coordinates defined in Eq. (16). Note that even though we defined the Berry phase in terms of the PT Berry connections, it is gauge-independent, in agreement with our discussions in Section 2. A periodic structure gauge is obtained by distributing this Berry phase evenly along the path onto the cell-periodic states in the PT gauge,

(37)$$|u_n^{\textbf{k}}\rangle \equiv {e^{- i\varphi _{n,{\kappa _d}}^B{\kappa _d}/\left\| {{{\textbf{b}}_d}} \right\|}}|\bar u_n^{\textbf{k}}\rangle .$$

For 2D or 3D systems ($D = 2,3$), the above procedure is then repeated along all the dimensions $d$. This gauge is denoted as the TPT gauge, and constitutes a periodic gauge with optimally smooth phase variation of the Bloch states [58].

It is important to mention that a globally smooth periodic gauge is only possible for topologically trivial systems. For example, the Chern theorem states that the Berry phase [see Eq. (9)] calculated over a closed surface $S$ is quantized as a integer multiple of $2\pi$, which in the case of the Bloch problem reads,

(38)$${C_n} = (2\pi {)^{- 1}}\mathop{{\iint }\mkern-21mu \bigcirc}\nolimits_S {\boldsymbol{\Omega} _n^{\textbf{k}}} \cdot {\rm d}{\textbf{S}},$$

with $\boldsymbol{\Omega} _n^{\textbf{k}}$ the Berry curvature [Eq. (10)], and ${C_n}$ is called the Chern number or the topological invariant. For a 2D system where the surface is the whole BZ, a nonzero Chern number (topologically nontrivial systems) presents a topological obstruction to the construction of a globally smooth periodic structure gauge [58]. For example, for topological systems, one can always construct the procedure for the TPT gauge along one dimension, but when one subsequently constructs the TPT gauge along the other dimension, the periodicity along the first dimension will break down somewhere in the BZ. As we will discuss later in Section 4 on the time-dependent problem, depending on the laser gauge, there are ways to circumvent the construction of such a global periodic gauge for HHG calculations.

E. Wannier Gauge

Up to this point, we have considered the Bloch states, which are entirely delocalized spatially. Often, however, it is useful to consider a spatially localized basis mimicking that of atomic or molecular orbitals. The construction of a periodic gauge opens up the possibility of constructing the Wannier basis [78,79], which consists of states localized on the individual lattice sites. The Wannier states are defined as the ${\textbf{k}}$ space Fourier transforms of the Bloch states,

(39)$$|w_n^{\textbf{R}}\rangle = {\cal F}\left\{{|\phi _n^{\textbf{k}}\rangle} \right\},\quad {\rm{with}}\quad {\cal F}\left\{\cdot \right\} \equiv \frac{{{V_{{\rm{cell}}}}}}{{{{(2\pi)}^D}}}\int_{{\rm{BZ}}} {e^{- i{\textbf{k}} \cdot {\textbf{R}}}}\{\cdot \} {\rm d}{\textbf{k}},$$

and ${\textbf{R}}$ is a Bravais lattice vector. As mentioned, the Wannier functions are localized, in the sense that $w_n^{\textbf{R}}({\textbf{r}}) = w_n^{\textbf{0}}({\textbf{r}} - {\textbf{R}})$, and $|w_n^{\textbf{R}}({\textbf{r}})| \to 0$ for $|{\textbf{r}} - {\textbf{R}}| \to \infty$.

Clearly, to perform the integral over the BZ in Eq. (39), the employed Bloch states in the integrand should be smooth and periodic; i.e., they should be precalculated in a periodic gauge. Since such a gauge choice is not unique, the resulting Wannier states depend on the choice of the periodic gauge. As the TPT structure gauge is the optimally smooth periodic gauge, the Wannier states constructed using the TPT gauge Bloch states are optimally localized (with a minimum spatial spread), and are called the maximally localized Wannier functions (MLWF) [59,80].

Independent of the gauge, the Wannier states are orthonormal ${\langle w_n^{\textbf{R}}|w_{{n^\prime}}^{{\textbf{R}}^\prime}\rangle _{{\rm{crys}}}} = {\delta _{{nn^\prime}}}{\delta _{{\textbf{R}},{\textbf{R}}^\prime}}$, and the Hilbert spaces spanned by the Bloch states and the Wannier states are identical. When working with Wannier states as a basis, the terminology “Wannier gauge” [59,81] is used, as the Wannier states are generally not energy eigenstates. Instead, the Hamiltonian matrix elements in the Wannier basis are the Fourier transforms of the band energies,

(40)$${\langle w_n^{\textbf{0}}|{\hat H_0}|w_n^{\textbf{R}}\rangle _{{\rm{crys}}}} = {\cal F}\left\{{E_n^{\textbf{k}}} \right\}.$$

The inverse transform of Eq. (40), $E_n^{\textbf{k}} = {\langle \phi _n^{\textbf{k}}|{\hat H_0}|\phi _n^{\textbf{k}}\rangle _{{\rm{crys}}}} = \sum\nolimits_{\textbf{R}} {e^{i{\textbf{k}} \cdot {\textbf{R}}}}{\langle w_n^{\textbf{0}}|{\hat H_0}|w_n^{\textbf{R}}\rangle _{{\rm{crys}}}}$, shows that the Wannier states considered as the localized basis in tight-binding models exactly reproduces the band energies. The position matrix elements in the Wannier gauge are well defined, and are just the Fourier transforms of the Berry connections,

(41)$${\langle w_n^{\textbf{0}}|\hat{\textbf{r}}|w_n^{\textbf{R}}\rangle _{{\rm{crys}}}} = {\cal F}\left\{{{\cal A}_n^{\textbf{k}}} \right\}.$$

The Wannier centers, ${\langle w_n^{\textbf{0}}|\hat{\textbf{r}}|w_n^{\textbf{0}}\rangle _{{\rm{crys}}}}$, are invariant with respect to a gauge change of the Bloch states, and in 1D are given by $a\varphi _n^B/(2\pi)$, with ${\varphi ^B}$ Zak’s phase and $a$ the lattice constant. However, as mentioned, the spread of the Wannier states are gauge-dependent.

F. Structure Gauges for Degenerate Bands

The gauge discussion in the previous subsections can be generalized to the degenerate case [58], which we sketch here for completeness. Let ${\{E_n^{\textbf{k}}\} _J}$ be a set of $J$ bands that are isolated from all other bands, in the sense that they have no degeneracies with the other bands anywhere in the BZ, but can have degeneracies within themselves. The Hilbert subspace spanned by the corresponding $J$ states $|u_n^{\textbf{k}}\rangle$ remains unchanged under an unitary transformation

(42)$$|\bar u_n^{\textbf{k}}\rangle = \sum\limits_{m = 1}^J U_{\textit{mn}}^{\textbf{k}}|u_m^{\textbf{k}}\rangle ,$$

and the objective is to pick ${U^{\textbf{k}}}$ such that the transformed states have the desired properties of our structure gauge.

We first discuss the construction of the PT gauge states ${\{\bar u_n^{\textbf{k}}\} _J}$ from a random gauge ${\{\breve{u} _n^{\textbf{k}}\} _J}$, and consider the discrete ${\textbf{k}}$ case, which is illustrative and useful for numerical applications. The procedure can be considered a generalization of Section 3.C. Starting at ${{\textbf{k}}_0}$ with $|\bar u_m^{{{\textbf{k}}_0}}\rangle \equiv |\breve{u} _m^{{{\textbf{k}}_0}}\rangle$, the overlap matrix with a neighboring point ${{\textbf{k}}_1} = {{\textbf{k}}_0} + \delta {\textbf{k}}$ can be decomposed by a singular value decomposition

(43)$$M_{\textit{mn}}^{{{\textbf{k}}_0}{{\textbf{k}}_1}} \equiv \langle \bar u_m^{{{\textbf{k}}_0}}|\breve{u} _n^{{{\textbf{k}}_1}}\rangle = (V\Sigma {W^\dagger}{)_{\textit{mn}}},$$

with $\Sigma$ a diagonal matrix with nonnegative values, and $V$, $W$ unitary matrices. The difference between $\Sigma$ and the identity matrix is a measure of the difference between the Hilbert space spanned by ${\{u_n^{{{\textbf{k}}_0}}\} _J}$ and ${\{u_n^{{{\textbf{k}}_1}}\} _J}$. We now define a new set of states by transforming the set of states from the random gauge,

(44)$$|\bar u_n^{{{\textbf{k}}_1}}\rangle = \sum\limits_{m = 1}^J {(W{V^\dagger})_{\textit{mn}}}|\breve{u} _m^{{{\textbf{k}}_1}}\rangle .$$

The overlap matrix can be shown using Eqs. (43) and (44) to be $\langle \bar u_m^{{{\textbf{k}}_0}}|\bar u_n^{{{\textbf{k}}_1}}\rangle = (V\Sigma {V^\dagger}{)_{\textit{mn}}}$ and is now Hermitian and positive definite, and is “as close as possible“ to the identity matrix, such that the constructed set of states ${\{\bar u_n^{{{\textbf{k}}_0}}\} _J}$ can be considered as “optimally aligned” to ${\{\bar u_n^{{{\textbf{k}}_1}}\} _J}$. The PT gauge of a local region in the BZ can now be constructed by repeating the procedure above to all other nearby points in the BZ.

We proceed to discuss a procedure to construct a periodic gauge in the multiband case, which is a generalization of the nondegenerate TPT gauge discussed in Section 3.D. First construct the PT gauge with the above procedure, starting from ${{\textbf{k}}_0}$ at one end of the BZ along the reciprocal vector direction $\hat{{\textbf{b}}}_d$ to the other end at ${{\textbf{k}}_{N - 1}}$. Even though the Hamiltonian at ${{\textbf{k}}_0}$ and ${{\textbf{k}}_N}$ are identical, the set of states ${\{\bar \phi _n^{{{\textbf{k}}_0}}\} _J}$ and ${\{\bar \phi _n^{{{\textbf{k}}_N}}\} _J}$ generally are not. In analogy with the nondegenerate case, when wrapping around the BZ, we can define the unitary overlap matrix

(45)$${{\cal U}_{\textit{mn}}} = (\Delta {\kappa _d}{)^{- 1}}\langle \bar u_m^{{{\textbf{k}}_{N - 1}}}|{e^{- i{{\textbf{b}}_d} \cdot {\textbf{r}}}}|\bar u_n^{{{\textbf{k}}_0}}\rangle ,$$

from which we can define a Berry phase

(46)$$\Phi _{{\kappa _d}}^B = - {\rm{Im}}\ln \det {\cal U} = \sum\limits_{n = 1}^J \varphi _{n,{\kappa _d}}^B.$$

The $\varphi _{n,{\kappa _d}}^B$ are the argument of the eigenvalues of ${\cal U}$, and interpreted as the Berry phases for the individual bands.

Similar to the single-band case discussed in Section 3.D, a periodic gauge can then be constructed from the PT gauge by dividing the Berry phases evenly along the path

(47)$$|u_n^{\textbf{k}}\rangle = {e^{- i\varphi _{n,{\kappa _d}}^B{\kappa _d}/\left\| {{{\textbf{b}}_d}} \right\|}}|\bar u_n^{\textbf{k}}\rangle ,$$

which is the TPT gauge for the multiband case. It should be noted that the states in the multiband PT and TPT gauges, ${\{\bar \phi _n^{\textbf{k}}\} _J}$ and ${\{\phi _n^{\textbf{k}}\} _J}$, are generally not energy eigenstates.

Once the multiband TPT gauge is constructed, the Wannier gauge in the multiband case can then be constructed from the TPT gauge as in the single-band case, i.e., using Eq. (39). The Wannier functions retain most of the important properties from the single-band case, such as localization in real space and lattice periodicity [59].

4. TIME-DEPENDENT PROBLEM AND LASER GAUGES

In this section, we discuss the different avenues to tackle the time-dependent problem of a solid interacting with a strong laser field, by treating the crystal quantum mechanically and the external field classically. Depending on the chosen laser gauge and basis, the resulting EOMs will have their own advantages and drawbacks. The section starts with the introduction to laser gauge freedom, and proceeds in Sections 4.A–4.C to obtain the relevant EOMs and equations for the microscopic current. Section 4.D compares the different time-dependent methods, in terms of both the numerical complexity and the interpretation of the physics, and Section 4.E describes the calculation of the HHG spectral and temporal profiles and provides an example.

The minimal coupling Hamiltonian for a nonrelativistic electron in a periodic potential interacting with an electromagnetic field reads [73,82]

(48)$$\hat H(t) = \frac{1}{2}{[\hat{\textbf{p}} + {\textbf{A}}({\textbf{r}},t)]^2} - \Phi ({\textbf{r}},t) + V({\textbf{r}}),$$

with ${\textbf{A}}({\textbf{r}},t)$ the vector potential, $\Phi ({\textbf{r}},t)$ the electric scalar potential, and the physical fields given by

(49)$${\textbf{F}}({\textbf{r}},t) = - \nabla \Phi ({\textbf{r}},t) - {\partial _t}{\textbf{A}}({\textbf{r}},t),$$

(50)$${\textbf{B}}({\textbf{r}},t) = \nabla \times {\textbf{A}}({\textbf{r}},t).$$

There is a gauge freedom in choosing ${\textbf{A}}$ and $\Phi$, as the physical fields remain invariant under the transformations

(51a)$${\textbf{A}}({\textbf{r}},t) \to {\textbf{A}}({\textbf{r}},t) + \nabla \Lambda ({\textbf{r}},t),$$(51b)$$\Phi ({\textbf{r}},t) \to \Phi ({\textbf{r}},t) - {\partial _t}\Lambda ({\textbf{r}},t),$$

with $\Lambda ({\textbf{r}},t)$ a differentiable real function. It is straightforward to show that under the gauge transform (51), the TDSE remains invariant if the wave function transforms as

(52)$$\Psi ({\textbf{r}},t) \to \Psi ^\prime ({\textbf{r}},t) = {e^{- i\Lambda ({\textbf{r}},t)}}\Psi ({\textbf{r}},t),$$

(53)$$\begin{split}\hat H(t) \to \hat H^\prime (t)& = \frac{1}{2}{[\hat{\textbf{p}} + {\textbf{A}}({\textbf{r}},t) + \nabla \Lambda ({\textbf{r}},t)]^2} \\&\quad- \Phi ({\textbf{r}},t) + {\partial _t}\Lambda ({\textbf{r}},t) + V({\textbf{r}}),\end{split}$$

where the Hamiltonian transform is explicitly written down in the second line, obtained by insertion of the gauge transform in Eq. (51) into the TDSE. Expectation values of physical observables such as the position $\langle \Psi (t)|\hat{\textbf{r}}|\Psi (t)\rangle$ and kinetic momenta $\langle \Psi (t)|[\hat{{\textbf{p}} + {\textbf{A}}(t)}]|\Psi (t)\rangle$ are gauge invariant under the transform.

The microscopic current operator is proportional to the kinetic momentum

(54)$$\hat{\textbf{j}}(t) = - \hat{\textbf{v}}(t) = i\left[\hat{{\textbf{r}}},\hat{H}(t) \right] = - [\hat{\textbf{p}} + {\textbf{A}}({\textbf{r}},t)],$$

where $[\cdot , \cdot]$ denotes a commutator. The wavelengths of driving fields used for HHG are in the order of micrometers, while the unit cell dimension is sub-nanometer (nm), so in the following we apply the dipole approximation ${\textbf{F}}(t) \equiv {\textbf{F}}({\textbf{r}},t)$ and ignore the magnetic field. In Section 6 we will discuss situations beyond the dipole approximation.

A. Velocity Gauge EOM in the Bloch Basis

We start by considering the dynamics in the VG, in which one fixes the gauge by choosing ${\Phi _{\rm{VG}}}({\textbf{r}},t) = 0$ and ${{\textbf{A}}_{\rm{VG}}}({\textbf{r}},t) = - \int {\textbf{F}}(t^\prime){\rm d}t^\prime \equiv {\textbf{A}}(t)$. As we will show, the VG can be advantangeous because it leads to EOMs that can be propagated separately for each ${\textbf{k}}$, which means that the construction of a periodic structure gauge prior to time propagation is unnecessary. The VG Hamiltonian reads

(55)$${\hat H_{\rm{VG}}}(t) = \frac{1}{2}{[\hat{\textbf{p}} + {\textbf{A}}(t)]^2} + V({\textbf{r}}).$$

The ${{\textbf{A}}^2}/2$ term can be transformed away by choosing $\Lambda ({\textbf{r}},t) = - \frac{1}{2}\int {\textbf{A}}{(t^\prime)^2}{\rm d}t^\prime $ in Eqs. (52) and (53). The resulting VG Hamiltonian form,

(56)$${\hat H_{\rm{VG}}}^\prime (t) = \hat T + V({\textbf{r}}) + \hat{\textbf{p}} \cdot {\textbf{A}}(t),$$

is often used instead of ${\hat H_{\rm{VG}}}$ in Eq. (56).

The VG TDSE $i|{\dot \Psi _{\rm{VG}}}(t)\rangle = \hat H{^\prime _{\rm{VG}}}(t)|{\Psi _{\rm{VG}}}(t)\rangle$ can be rewritten by first expanding the wave function ${\Psi _{\rm{VG}}}({\textbf{r}},t) = {N^{- 1}}\sum\nolimits_{m,{\textbf{k}} \in BZ} b_m^{\textbf{k}}(t)\phi _m^{\textbf{k}}({\textbf{r}})$ and then projecting onto the Bloch states, resulting in the EOM for the coefficients [41,83]

(57)$$i\dot b_m^{\textbf{k}}(t) = E_m^{\textbf{k}}b_m^{\textbf{k}}(t) + {\textbf{A}}(t) \cdot \sum\limits_n {\textbf{p}}_{\textit{mn}}^{\textbf{k}}b_n^{\textbf{k}}(t),$$

where we have used the property that the momentum matrix elements are diagonal in ${\textbf{k}}$ [Eq. (21)], and defined ${\textbf{p}}_{\textit{mn}}^{\textbf{k}} \equiv {\langle \phi _m^{\textbf{k}}|\hat{\textbf{p}}|\phi _n^{\textbf{k}}\rangle _{{\rm{cell}}}}$. We can define the density matrix operator $\hat g(t) = |{\Psi _{\rm{VG}}}(t)\rangle \langle {\Psi _{\rm{VG}}}(t)|$, such that the matrix elements $g_{\textit{mn}}^{\textbf{k}} \equiv b_m^{\textbf{k}}b_n^{{\textbf{k}}*}$ satisfy the EOM [56]

(58)$$\begin{split}i\dot g_{\textit{mn}}^{\textbf{k}}(t) &=\left({E_m^{\textbf{k}} - E_n^{\textbf{k}}} \right)g_{\textit{mn}}^{\textbf{k}}(t) + {\textbf{A}}(t) \\&\quad\cdot \sum\limits_l \left[{{\textbf{p}}_{\textit{ml}}^{\textbf{k}}g_{\textit{ln}}^{\textbf{k}}(t) - {\textbf{p}}_{\textit{ln}}^{\textbf{k}}g_{\textit{ml}}^{\textbf{k}}(t)} \right].\end{split}$$

These equations also can be derived directly from the Liouville–von Neumann equation $i{\dot {\hat g}}(t) = [\hat H(t),\hat g(t)]$. Finally, the relevant observable for HHG, the microscopic current, is evaluated in the VG as

(59)$$\begin{split}{{\textbf{j}}(t) = {\rm{Tr}}\left[{{\hat{{\textbf{j}}}_{\rm{VG}}}(t)\hat g(t)} \right] = - {\rm{Tr}}\left\{{\left[\hat{{\textbf{p}}} + {\textbf{A}}(t) \right]\hat g(t)} \right\}\\= - {N^{- 1}}\sum\limits_{{\textbf{k}} \in BZ} \sum\limits_{\textit{mn}} \left[{{\textbf{p}}_{\textit{mn}}^{\textbf{k}} + {\delta _{\textit{mn}}}{\textbf{A}}(t)} \right]g_{\textit{nm}}^{\textbf{k}}(t),}\end{split}$$

where we have used Eq. (54). We will discuss the advantages and drawbacks of time propagation in the VG in Section 4.D.

B. Length Gauge EOM in Bloch Basis

A different way to fix the laser gauge in the dipole approximation is to set ${\Phi _{\rm{LG}}}({\textbf{r}},t) = - {\textbf{r}} \cdot {\textbf{F}}(t)$ and ${{\textbf{A}}_{\rm{LG}}} = {\textbf{0}}$, which according to Eq. (49) results in the desired physical field ${\textbf{F}}(t)$. This gauge choice is denoted as the LG. The LG Hamiltonian is then, according to Eq. (48),

(60)$${\hat H_{\rm{LG}}}(t) = \hat T + V({\textbf{r}}) + {\textbf{r}} \cdot {\textbf{F}}(t).$$

This Hamiltonian also can be obtained from the VG Hamiltonian using the gauge transform ${\Lambda _{VG \to LG}} = - {\textbf{A}}(t) \cdot {\textbf{r}}$ in Eqs. (52) and (53).

We first rewrite the LG TDSE $i|{\dot \Psi _{\rm{LG}}}(t)\rangle = {\hat H_{\rm{LG}}}(t)|{\Psi _{\rm{LG}}}(t)\rangle$ by expanding ${\Psi _{\rm{LG}}}({\textbf{r}},t) = {N^{- 1}}\sum\nolimits_{m,{\textbf{k}} \in BZ} a_m^{\textbf{k}}(t)\phi _m^{\textbf{k}}({\textbf{r}})$ and projecting onto the Bloch states, resulting in the EOM [38,84,85]

(61)$$i\dot a_m^{\textbf{k}}(t) = E_m^{\textbf{k}}a_m^{\textbf{k}}(t) + {\textbf{F}}(t) \cdot \sum\limits_n {\textbf{d}}_{\textit{mn}}^{\textbf{k}}a_n^{\textbf{k}}(t) + i{\textbf{F}}(t) \cdot {\nabla _{\textbf{k}}}a_m^{\textbf{k}}(t).$$

In contrast to the VG in Eq. (57), the LG equations couple different crystal momenta to each other due to the last term involving the ${\textbf{k}}$ gradient. Defining the LG density matrix $\hat \rho (t) = |{\Psi _{\rm{LG}}}(t)\rangle \langle {\Psi _{\rm{LG}}}(t)|$ with matrix elements $\rho _{\textit{mn}}^{\textbf{k}} = a_m^{\textbf{k}}a_n^{{\textbf{k}}*}$, the density matrix EOM reads

(62)$$\begin{split}i\dot \rho _{\textit{mn}}^{\textbf{k}}(t) &=({E_m^{\textbf{k}} - E_n^{\textbf{k}}} )\rho _{\textit{mn}}^{\textbf{k}}(t) + {\textbf{F}}(t) \\&\quad\cdot \sum\limits_l \left[{{\textbf{d}}_{\textit{ml}}^{\textbf{k}}\rho _{\textit{ln}}^{\textbf{k}}(t) - {\textbf{d}}_{\textit{ln}}^{\textbf{k}}\rho _{\textit{ml}}^{\textbf{k}}(t)} \right] + i{\textbf{F}}(t) \cdot {\nabla _{\textbf{k}}}\rho _{\textit{mn}}^{\textbf{k}}(t).\end{split}$$

As in the VG, these LG equations also can be derived from the Liouville–von Neumann equation $i{\dot{ \hat \rho}} (t) = [{\hat H_{\rm{LG}}}(t),\hat \rho (t)]$ in the Bloch basis. In a many-body framework using a second quantization, a similar equation for the one-body reduced density matrix $\rho _{\textit{mn}}^{\textbf{k}}(t) \equiv \langle \Psi (t)|a_{n{\textbf{k}}}^\dagger {a_{m{\textbf{k}}}}|\Psi (t)\rangle$ can be derived (with $a_{n{\textbf{k}}}^\dagger$ creation and ${a_{m{\textbf{k}}}}$ annihilation operators), often denoted as the semiconductor Bloch equations (SBEs) [13,52–54,65,86,87]. Due to many-body couplings such as electron–electron and electron–phonon scattering, the SBEs will dephase, which at our level of theory is treated by adding a phenomenological dephasing term $(1 - {\delta _{\textit{mn}}})\rho _{\textit{mn}}^{\textbf{k}}/{T_2}$ on the right-hand side of Eq. (62). The SBEs using a phenomenological dephasing and neglecting excitons [53] are equivalent to the time-dependent density matrix equations presented here. A shorter dephasing time ${T_2}$ will result in less noisy HHG spectra, and ${T_2}$ is often chosen such that there is reasonable agreement between experiment and theory [38,88]. The phenomenogical dephasing, while being computationally convenient, can only give qualitative results, and more accurate treatments of dephasing is an active area of research [55,89–91].

For degenerate subspaces given in the periodic gauge (see Section 3.F), the rotated states of Eq. (42) are not energy eigenstates, and the EOM given in Eq. (62) cannot be applied. In this case, the Liouville–von Neumann expression should be used to reduce the correct propagation equations [65,92]. The correct propagation equations in the Wannier structure gauge (see Section 3.E) can also be derived using the same strategy [87,93].

In the LG, the current operator is $\hat{{\textbf{j}}}_{\rm{LG}}(t) = - \hat{\textbf{p}}$ [Eq. (54)], and the microscopic current can be evaluated as

(63)$${{\textbf{j}}(t) = {\rm{Tr}}[{{\hat{{\textbf{j}}}_{\rm{LG}}}(t)\hat \rho (t)}] = - {N^{- 1}}\sum\limits_{{\textbf{k}} \in BZ} \sum\limits_{\textit{mn}} {\textbf{p}}_{\textit{mn}}^{\textbf{k}}\rho _{\textit{nm}}^{\textbf{k}}(t).}$$

While Eq. (63) is useful to obtain the total current, more physical insight can be gained by a decomposition of the current into several terms with different physical meanings. While the decomposition into intraband and interband currents has been discussed extensively in the literature, we present a decomposition here that contains four terms each with its own physical interpretation [63,94]. We first split the position operator into an intraband and an interband component, $\hat{\textbf{r}} = \hat{{\textbf{r}}}^{{\rm{tra}}} + \hat{{\textbf{r}}}^{{\rm{ter}}}$, with the matrix elements in the Bloch basis [see Eqs. (27)–(29)]:

(64a)$${\langle \phi _m^{\textbf{k}}|\hat{{\textbf{r}}}^{{\rm{tra}}}|\phi _n^{\textbf{q}}\rangle _{{\rm{crys}}}} = N{\delta _{\textit{mn}}}({\cal A}_m^{\textbf{k}} + i{\nabla _{\textbf{k}}}){\delta _{{\textbf{kq}}}},$$(64b)$${\langle \phi _m^{\textbf{k}}|\hat{{\textbf{r}}}^{{\rm{ter}}}|\phi _n^{\textbf{q}}\rangle _{{\rm{crys}}}} = N(1 - {\delta _{\textit{mn}}}){\delta _{{\textbf{kq}}}}{\textbf{d}}_{\textit{mn}}^{\textbf{k}}.$$

The current can now be written

(65)$$\begin{split}{\textbf{j}}(t) &= {\rm{Tr}} \{{{\hat{{\textbf{j}}}_{\rm{LG}}}(t)\hat \rho (t)} \} = i{\rm{Tr}} \{{[\hat{\textbf{r}},{{\hat H}_{\rm{LG}}}(t)]\hat \rho (t)} \}\\ &= i{\rm{Tr}} ( \{[{\hat{{\textbf{r}}}^{{\rm{tra}}}},{{\hat H}_0}] + [{\hat{{\textbf{r}}}^{{\rm{tra}}}},{\textbf{F}}(t) \cdot {\hat{{\textbf{r}}}^{{\rm{tra}}}}] \\&\quad+ [{\hat{{\textbf{r}}}^{{\rm{tra}}}},{\textbf{F}}(t) \cdot {\hat{{\textbf{r}}}^{{\rm{ter}}}}] + [{\hat{{\textbf{r}}}^{{\rm{ter}}}},{{\hat H}_{\rm{LG}}}(t)] \}\hat \rho )\\ &\equiv {{{\textbf{j}}^{{\rm{tra}}}}(t) + {{\textbf{j}}^{{\rm{anom}}}}(t) + {{\textbf{j}}^{{\rm{mix}}}}(t) + {{\textbf{j}}^{{\rm{pol}}}}(t),}\end{split}$$

where we have inserted the LG Hamiltonian (60). After some tedious, but straightforward derivations using Eqs. (14) and (64), the four current components in Eq. (65) can be written as [63]

(66a)$${{\textbf{j}}^{{\rm{tra}}}}(t) = - {N^{- 1}}\sum\limits_{m{\textbf{k}}} {\nabla _{\textbf{k}}}E_m^{\textbf{k}}\rho _{\textit{mm}}^{\textbf{k}}(t),$$(66b)$${{\textbf{j}}^{{\rm{pol}}}}(t) = - {N^{- 1}}{\partial _t}\sum\limits_{m \ne n,{\textbf{k}}} {\textbf{d}}_{\textit{mn}}^{\textbf{k}}\rho _{\textit{nm}}^{\textbf{k}}(t),$$(66c)$${{\textbf{j}}^{{\rm{anom}}}}(t) = - {N^{- 1}}\sum\limits_{m{\textbf{k}}} \left[{{\textbf{F}}(t) \times \boldsymbol{\Omega} _m^{\textbf{k}}} \right]\rho _{\textit{mm}}^{\textbf{k}}(t),$$(66d)$$\begin{split}{\textbf{j}}^{{\rm{mix}}}(t) &= - {N^{- 1}}\sum\limits_\mu {F_\mu}(t)\sum\limits_{m \ne n,{\textbf{k}}}\\&\quad\times \left[{\left({{\nabla _{\textbf{k}}}d_{\mu ,mn}^{\textbf{k}}} \right) - i({\cal A}_m^{\textbf{k}} - {\cal A}_n^{\textbf{k}})d_{\mu ,mn}^{\textbf{k}}} \right]\rho _{\textit{nm}}^{\textbf{k}}(t),\end{split}$$

where ${F_\mu}(t)$ and $d_{\mu ,mn}^{\textbf{k}}$ are, respectively, the $\mu$th component of the field and dipole matrix elements. In Eq. (66a), the intraband current ${{\textbf{j}}^{{\rm{tra}}}}(t)$ is due to carrier transport within individual bands as is reflected in its dependence on the carrier group velocities ${\nabla _{\textbf{k}}}E_m^{\textbf{k}}$. The interband current ${{\textbf{j}}^{{\rm{pol}}}}(t)$ in Eq. (66b) originates from the coupling between the bands and is seen to be just the time-derivative of the polarization. The anomalous current ${{\textbf{j}}^{{\rm{anom}}}}(t)$ depends on the Berry curvature $\boldsymbol{\Omega} _m^{\textbf{k}} = {\nabla _{\textbf{k}}} \times {\cal A}_m^{\textbf{k}}$ [see Eq. (10)], and is perpendicular to the electric field ${\textbf{F}}(t)$. The mixture current ${{\textbf{j}}^{{\rm{mix}}}}(t)$ in Eq. (66b), as seen from the definition in Eq. (65), depends on the coupling between the interband and intraband position operators. The structure-gauge-invariant expression in the square parenthesis is also referred to as the generalized derivative [63], defined as ${(O_{\textit{mn}}^{\textbf{k}})_{;{\textbf{k}}}} \equiv {\nabla _{\textbf{k}}}O_{\textit{mn}}^{\textbf{k}} - i({\cal A}_m^{\textbf{k}} - {\cal A}_n^{\textbf{k}})O_{\textit{mn}}^{\textbf{k}}$. We note here that even though the ${{\textbf{j}}^{{\rm{anom}}}}(t)$ contribution to the current originated in the intraband component of the position operator in Eq. (65), it and the mixture term ${{\textbf{j}}^{{\rm{mix}}}}(t)$ depend explicitly on coherences between different bands and, as such, are interband in nature.

More generally, it is worth mentioning that in the literature, the decomposition in terms of the intraband [Eq. (66a)] and interband currents [Eq. (66b)] has been discussed extensively for HHG in solids [13,14,17,23,38,95–97]. The effect from the anomalous velocity and the Berry curvature has been investigated mostly separately [21,33,98]. The implication of the mixture terms in Eq. (66d) has largely been unexplored [94]. Often, the non-intraband current (containing the interband, anomalous, and mixture currents) is used interchangeably with the interband current, either due to semantics, or because the anomalous and mixture currents often are negligible. Sometimes, the anomalous current is also considered as a part of the intraband current [34,99], since it depends on the carrier density in a given band [see Eq. (66c)] and can be considered a band-specific current. Remember, however, that the Berry curvature is a geometric property stemming from the residual coupling between the bands [see Eq. (11)]. The interplay between the four current contributions investigated under a combined framework could be a potential avenue of future research.

C. EOM in the Adiabatic Houston Basis

Next, we discuss a popular propagation method that employs the Houston states, which, like the LG EOMs discussed in the previous subsection, can naturally incorporate a phenomenological description of dephasing. We first reformulate the time-dependent problem into one that draws parallel to the general adiabatic theory discussed in Section 2. First realize that we can define a transformed time-independent Hamiltonian, $\hat H_0^{\textbf{k}} \equiv {e^{- i{\textbf{k}} \cdot {\textbf{r}}}}{\hat H_0}{e^{i{\textbf{k}} \cdot {\textbf{r}}}}$, with the cell-periodic functions as eigenstates and the same band energies as the Bloch states:

(67)$$\hat H_0^{\textbf{k}}u_n^{\textbf{k}}({\textbf{r}}) = E_n^{\textbf{k}}u_n^{\textbf{k}}({\textbf{r}}).$$

For the time-dependent VG Hamiltonian in Eq. (55), the transformed Hamiltonian can easily be shown to satisfy

(68)$${\hat H^{\textbf{K}}}(t) \equiv {e^{- i{\textbf{K}} \cdot {\textbf{r}}}}{\hat H_{\rm{VG}}}(t){e^{i{\textbf{K}} \cdot {\textbf{r}}}} = \hat H_0^{{\textbf{K}} + {\textbf{A}}(t)}.$$

Using Eqs. (67) and (68), we see that a set of adiabatic states [Eq. (1)] of ${\hat H^{\textbf{K}}}(t)$ simply consists of the cell-periodic functions with shifted crystal momenta, $u_n^{{\textbf{K}} + {\textbf{A}}(t)}({\textbf{r}})$, and with corresponding adiabatic eigenenergies $E_n^{\textbf{K}}(t) = E_n^{{\textbf{K}} + {\textbf{A}}(t)}$. The solution to the adiabatic problem of the original VG Hamiltonian can now be written as

(69)$$\begin{split}{{{\hat H}^{\textbf{K}}}(t)u_n^{{\textbf{K}} + {\textbf{A}}(t)}({\textbf{r}}) ={E^{{\textbf{K}} + {\textbf{A}}(t)}}u_n^{{\textbf{K}} + {\textbf{A}}(t)}({\textbf{r}})}\\{\Leftrightarrow {{\hat H}_{\rm{VG}}}(t)h_n^{\textbf{K}}({\textbf{r}},t) = {E^{{\textbf{K}} + {\textbf{A}}(t)}}h_n^{\textbf{K}}({\textbf{r}},t)}\end{split}$$

with

(70)$$h_n^{\textbf{K}}({\textbf{r}},t) \equiv {e^{i{\textbf{k}} \cdot {\textbf{r}}}}u_n^{{\textbf{K}} + {\textbf{A}}(t)}({\textbf{r}}) = {e^{- i{\textbf{A}} \cdot {\textbf{r}}}}\phi _n^{{\textbf{K}} + {\textbf{A}}(t)}({\textbf{r}}).$$

This set of adiabatic states $h_n^{\textbf{K}}({\textbf{r}},t)$ are often referred to as Houston states in the literature [36,100], and we have used the capital letter ${\textbf{K}}$ for the crystal momenta to highlight that the Houston states are the accelerated Bloch states. All the properties of the adiabatic states discussed in Section 2 now automatically follow. If the structure gauge for the Bloch states is fixed, then the Houston states are completely well-determined and given by Eq. (70). However, we emphasize that the Houston states only represent one possible set of adiabatic states. At each instant of time $t$, a Houston state multiplied by an arbitrary ${\textbf{k}}$ dependent phase factor is another adiabatic state. To rephrase, an adiabatic state of ${\hat H_{\rm{VG}}}(t)$, $|n(t)\rangle$, is always related to a Houston state by an arbitrary phase factor

(71)$$\langle {\textbf{r}}|n(t)\rangle = {e^{i\varphi _n^{\textbf{k}}}}h_n^{\textbf{K}}({\textbf{r}},t).$$

Since the adiabatic states are defined as the instantaneous eigenstates of the Hamiltonian involving the laser–matter interaction, they can also be considered as the laser-field-dressed states.

Expanding the wave function in the Houston basis, ${\Psi _{\rm{VG}}}({\textbf{r}},t) = {N^{- 1}}\sum\nolimits_{m,{\textbf{k}}} c_m^{\textbf{K}}(t)h_m^{\textbf{K}}({\textbf{r}},t)$, the VG TDSE reads

(72)$$i\dot c_m^{\textbf{K}}(t) = E_m^{{\textbf{K}} + {\textbf{A}}(t)}c_m^{\textbf{K}}(t) + {\textbf{F}}(t) \cdot \sum\limits_n {\textbf{d}}_{\textit{mn}}^{{\textbf{K}} + {\textbf{A}}(t)}c_n^{\textbf{K}}(t),$$

where we have calculated the nonadiabatic couplings [see Eq. (4)]

(73)$$\begin{split}i{\langle h_m^{\textbf{Q}}(t)|\dot h_n^{\textbf{K}}(t)\rangle _{{\rm{crys}}}} &= Ni{\delta _{{\textbf{Q}},{\textbf{K}}}}{\langle u_m^{{\textbf{Q}} + {\textbf{A}}(t)}|\dot u_n^{{\textbf{K}} + {\textbf{A}}(t)}\rangle _{{\rm{cell}}}} \\[-3pt]&= - N{\delta _{{\textbf{Q}},{\textbf{K}}}}{\textbf{F}}(t) \cdot {\textbf{d}}_{\textit{mn}}^{{\textbf{K}} + {\textbf{A}}(t)}.\end{split}$$

The corresponding density matrix elements in the Houston basis, $\bar \rho _{\textit{mn}}^{\textbf{K}} = c_m^{\textbf{K}}c_n^{{\textbf{K}}*}$, evolve as

(74)$$\begin{split}i {\dot{\bar \rho}} _{\textit{mn}}^{\textbf{K}}(t)& = ({E_m^{{\textbf{K}} + {\textbf{A}}(t)} - E_n^{{\textbf{K}} + {\textbf{A}}(t)}})\bar \rho _{\textit{mn}}^{\textbf{K}}(t) + {\textbf{F}}(t) \\&\quad\cdot \sum\limits_l [{{\textbf{d}}_{\textit{ml}}^{{\textbf{K}} + {\textbf{A}}(t)}\bar \rho _{\textit{ln}}^{\textbf{K}}(t) - {\textbf{d}}_{\textit{ln}}^{{\textbf{K}} + {\textbf{A}}(t)}\bar \rho _{\textit{ml}}^{\textbf{K}}(t)}].\end{split}$$

Similarly, as in the previous subsection, for a degenerate subspace expressed in a periodic gauge, the Liouville–von Neumann equation must be used to reduce the relevant EOMs. Note that in the literature, due to the appearance of the dipole operator in the EOMs of Eq. (74), it has also been termed as the LG SBEs in the moving frame [38,55,56,101]. Indeed, the EOMs in the VG adiabatic basis [Eqs. (61) and (62)] can be obtained from the LG EOMs [Eqs. (72) and (74)] by the frame change ${\textbf{k}} = {\textbf{K}} + {\textbf{A}}(t)$. As in the LG case of Section 4.B, dephasing can be introduced by a phenomenological term $(1 - {\delta _{\textit{mn}}})\bar \rho _{\textit{mn}}^{\textbf{K}}/{T_2}$ on the right-hand side of Eq. (74). We note that one cannot add this dephasing term directly to EOM in the Bloch-basis VG [Eq. (58)] due to the severe mixing of the field-free states in the presence of the laser field.

The microscopic current is evaluated as

(75)$$\begin{split}{\textbf{j}}(t) &= {\rm{Tr}}[{\hat{{\textbf{j}}}_{\rm{VG}}(t)\hat g(t)} ] = - {\rm{Tr}}\{{[\hat{{\textbf{p}}} + {\textbf{A}}(t) ]\hat g(t)} \} \\&= - {N^{- 1}}\sum\limits_{mn{\textbf{K}}} {\textbf{p}}_{\textit{mn}}^{{\textbf{K}} + {\textbf{A}}(t)}\bar \rho _{\textit{nm}}^{\textbf{K}}(t),\end{split}$$

where we have used the identity ${\langle h_m^{\textbf{K}}(t)|[\hat{{\textbf{p}}} + {\textbf{A}}(t)]|h_n^{\textbf{Q}}(t)\rangle _{{\rm{crys}}}} = {\delta _{{\textbf{KQ}}}}{\textbf{p}}_{\textit{mn}}^{{\textbf{K}} + {\textbf{A}}(t)}$. The current can be split into an intraband and a non-intraband contribution, by splitting the diagonal and off-diagonal terms of the momentum matrix elements in Eq. (75), i.e., ${\textbf{j}}(t) = {{\textbf{j}}^{{\rm{tra}}}}(t) + {{\textbf{j}}^{{\rm{ter}}}}(t)$:

(76a)$${{\textbf{j}}^{{\rm{tra}}}}(t) = - {N^{- 1}}\sum\limits_{n{\textbf{K}}} {\nabla _{\textbf{K}}}E_n^{{\textbf{K}} + {\textbf{A}}(t)}\bar \rho _{\textit{nn}}^{\textbf{K}},$$(76b)$${{\textbf{j}}^{{\rm{ter}}}}(t) = - {N^{- 1}}\sum\limits_{m \ne n,{\textbf{K}}} {\textbf{p}}_{\textit{mn}}^{{\textbf{K}} + {\textbf{A}}(t)}\bar \rho _{\textit{nm}}^{\textbf{K}},$$

where the non-intraband current contains the interband, anomalous, and mixture contributions to the current [compare to Eqs. (65) and (66)].

D. Gauge Comparisons

The three different time-dependent propagation methods presented in Sections 4.A–4.C each have their own advantages and drawbacks. Here we provide a discussion on this topic. The calculation procedures consist, in broad stokes, of five steps that are summarized in the flowchart in Fig. 2. These steps are: (1) Obtain the band structure and its associated momentum coupling matrix elements. This could be done either by using a structure code or diagonalizing an analytical Hamiltonian by, e.g., the tight-binding approximation. As discussed in Section 3.B, the output will in general be in a random structure gauge. If required, we can then construct a periodic structure gauge using the procedures outlined in Sections 3.C–3.F. (2) Next, we propagate the time-dependent EOMs, using one of the laser gauges discussed in Sections 4.A–4.C. During the time propagation, every ${n_{\delta t}}$th steps, we can apply the dephasing and calculate the time-dependent currents. (4) After the time-propagation, with full knowledge of the time-dependent current, we can calculate the HHG spectrum and obtain the time-frequency information of the harmonics (see Section 4.E). (5) Finally, in bulk crystals, the microscopic dynamics can be coupled with Maxwell’s equations to account for the macroscopic propagation of the laser and harmonics through the bulk, a topic we will revisit in Section 6.

Fig. 2. Flowchart sketching the computational steps for HHG in solids, consisting of a structure calculation, the propagation of the microscopic dynamics, the inclusion of dephasing and current, the spectrum calculation, and the macroscopic propagation.

Download Full Size | PDF

Going through the flowchart in more detail, we start by considering a VG description of the dynamics (see right-hand side of the flowchart). As discussed in Section 4.A, the VG EOMs represented in the Bloch basis [see Eq. (58)] are diagonal in the crystal momenta ${\textbf{k}}$, meaning that each ${\textbf{k}}$ can be propagated independently. Since no fixed phase relationship between neighboring $|\phi _n^{\textbf{k}}\rangle$ are required, there is no need to construct the periodic structure gauge (Section 3.D) prior to the time propagation, and using momentum matrix elements ${\textbf{p}}_{\textit{mn}}^{\textbf{k}}$ easily obtained within a random structure gauge (see Section 3.B) is sufficient to calculate the total current in Eq. (59) (see Fig. 2). If desired, it is possible to include the phenomenological dephasing effect also into the VG Bloch basis calculation. This can be achieved during the time propagation by first transforming into a VG adiabatic basis at desired times (separated by intervals $\Delta t = {n_{\delta t}}{\delta _t}$), applying the dephasing as

(77)$$\bar \rho _{\textit{mn}}^{\textbf{k}} \to \bar \rho _{\textit{mn}}^{\textbf{k}}{e^{- \Delta t/{T_2}}},\quad m \ne n,$$

and then transforming back to the Bloch basis. Here, we present two paths for such a basis transform. One method is to directly transform to the Houston basis [36,100]

(78)$${g_{\textit{mn}}^{\textbf{k}}(t) = \sum\limits_{\textit{lk}} Q_{\textit{ml}}^{\textbf{k}}(t)\bar \rho _{\textit{lk}}^{\textbf{k}}Q_{\textit{nk}}^{{\textbf{k}}*}(t),}$$

with $Q_{\textit{mn}}^{\textbf{k}}(t) \equiv {\langle u_m^{\textbf{k}}|u_n^{{\textbf{k}} + {\textbf{A}}(t)}\rangle _{{\rm{cell}}}}$, and where we have used the resolution of identity and ${\langle \phi _m^{\textbf{q}}|h_n^{\textbf{k}}(t)\rangle _{{\rm{crys}}}} = {\delta _{{\textbf{qk}}}}Q_{\textit{mn}}^{\textbf{k}}(t)$. This approach was used, for example, in [56,64,85]. Evaluation of Eq. (78) presents several potential difficulties: It requires the direct knowledge of the cell-periodic functions $|u_n^{\textbf{k}}\rangle$; the computational complexity is high [56]; and the periodic structure gauge is required since $|u_n^{\textbf{k}}\rangle$ must be splined at ${\textbf{k}} + {\textbf{A}}(t)$, unless an analytical form of $|u_n^{\textbf{k}}\rangle$ is known, which is unlikely for real materials beyond simple model systems (dashed arrow in Fig. 2). A different method to transform to an adiabatic basis is to use the definition in Eq. (1), i.e., by diagonalizing the instantaneous Hamiltonian in the Bloch basis. As discussed in Eq. (71), a resulting adiabatic state will differ from a Houston state by an arbitrary phase factor, which will not affect the inclusion of the phenomenological dephasing. While instantaneous diagonalization has similar computational complexity as Eq. (78), a random structure gauge is sufficient since each ${\textbf{k}}$ can be treated independently. Note that such a forward and backward transform is not required at each propagation time step $\delta t$, but rather at time step $\Delta t = {n_{\delta t}}{\delta _t}$, with ${n_{\delta t}} \gg 1$, reducing the computational complexity. In the VG Bloch basis, compared to the other two discussed methods, more bands are often required to achieve convergence [56,63–65,102], but at the same time degenerate bands can naturally be treated. In a future publication, we will go into more detail about treating HHG in the VG Bloch basis and decomposition of the current.

On the left-hand side of the flowchart, the LG approach involves the coupling of different crystal momenta by the ${\nabla _{\textbf{k}}}$ term, as discussed in Section 4.B, due to the non-diagonal nature of the position operator $\hat{\textbf{r}}$ in the Bloch basis. Hence, prior to the time propagation, the construction of a periodic structure gauge is required. The phenomenological dephasing term can be naturally included during propagation, and the total current can be decomposed into the four terms given in Eq. (65). For the VG EOMs in the Houston basis described in Section 4.B, a periodic structure gauge is also required since the dipole couplings must be splined at crystal momenta ${\textbf{K}} + {\textbf{A}}(t)$. The calculation of the numerical gradient or the splining generally requires a much finer ${\textbf{k}}$-space sampling than the VG Bloch scheme to obtain the same level of convergence, which severely increases the computational effort [56,65]. Similar to the LG case, the current can actually be decomposed into intraband and non-intraband terms, as presented in Eq. (76). Compared to the VG Bloch basis case, the LG EOM in the Bloch basis and VG EOMs in the Houston basis often require a smaller number of bands for convergence. Note that we have mostly discussed the construction of a globally smooth structure gauge, with the advantage that the gauge must only be constructed once, prior to the time propagation. In the literature, a method exists [65] that aims to construct a locally smooth structure gauge during the time propagation of the LG SBEs, which has recently been applied to HHG [92]. For electric fields with changing polarization directions, the gauge construction in all three Cartesian directions must be applied at each time step, potentially severely increasing the computational complexity. The key advantages of our three discussed propagation schemes are briefly summarized in Table 1.

Table 1. Relative Advantages of the VG Bloch Basis Propagation Scheme Compared to the LG Bloch Basis or Houston Basis Propagation Schemes

View Table

E. HHG Calculation and Time Profiles

The HHG spectral yield is proportional to the spectral intensity of the current, and is calculated as Larmor’s formula [103]

(79)$$S(\omega) \propto {\omega ^2}{\left\| {{\textbf{j}}(\omega)} \right\|^2},$$

with ${\textbf{j}}(\omega)$ the Fourier transform of the time-dependent current

(80)$${\textbf{j}}(\omega) = (2\pi {)^{- 1/2}}\int_{- \infty}^\infty {\textbf{j}}(t){e^{i\omega t}}{\rm d}t.$$

Often, a window function (mask) $w(t)$ is multiplied onto the current ${\textbf{j}}(t)$ to smoothly reduce it at large times (mimics current decay from scattering or other decay mechanisms) such that the integral can be done at finite times. For long pulses (${\sim}10$ optical cycles), the qualitative features of the HHG spectrum do not depend on the mask chosen, as long as the emission near the peak of the external field is retained. The emission intensity for harmonics polarized along a direction $\hat{\textbf{n}}$ is given by

(81)$${S_{\hat{{\textbf{n}}}}}(\omega) \propto {\omega ^2}{\left| {{\textbf{j}}(\omega) \cdot \hat{\textbf{n}}} \right|^2}.$$

Separate spectra for the decomposed currents can be obtained, for example, as ${S^{{\rm{tra}}}}(\omega) \propto {\omega ^2}{\left\| {{{\textbf{j}}^{{\rm{tra}}}}(\omega)} \right\|^2}$. Note, however, that generally $S(\omega) \ne {S^{{\rm{tra}}}}(\omega) + {S^{{\rm{anom}}}}(\omega) + {S^{{\rm{mix}}}}(\omega) + {S^{{\rm{pol}}}}(\omega)$ due to the existence of cross terms such as $j_\mu ^{{\rm{tra}}}(\omega)j_\mu ^{{\rm{pol}}}(\omega)$. The decomposition of the spectrum only makes sense at those $\omega$ where the cross terms are negligible. In practice, this is almost always the case.

More information on the time-frequency characteristics of the emitted harmonics can be obtained by a continuous wavelet transform (CWT) of the current

(82)$${S(t,a) \propto {a^{- 2}}{{\left\| {\int_{- \infty}^\infty {\textbf{j}}(t^\prime){W^*}\left[{\frac{{t^\prime - t}}{a}} \right]{\rm d}t^\prime} \right\|}^2}}$$

with $t$ a time-translation variable, $a$ a frequency-scaling variable, and $W(x)$ a mother wavelet. In our calculations, we use the Morlet–Grossman mother wavelet, with $W(x) = ({\sigma ^2}\pi {)^{- 1/4}}{e^{- i\Omega x}}{e^{- {x^2}/(2{\sigma ^2})}}$, where $\sigma$ is the standard deviation, and $\Omega$ the center frequency of the mother wavelet. Insertion of $x \equiv (t^\prime - t)/a$ into $W(x)$ shows that the frequency of the daughter wavelet is scaled as $\omega \equiv \Omega /a$, while the time-spread is scaled as ${\sigma _{{\rm{time}}}} = \sigma a$. Similarly, the daughter wavelet in the Fourier space has the spread ${\sigma _{{\rm{freq}}}} = (\sigma\! a{)^{- 1}}$. This scaling means that higher time resolution, along with lower frequency resolution, is employed for the higher frequencies in the HHG spectrum. This can be an advantage over a windowed Fourier transform where the time resolution and frequency resolution is fixed. The CWT can be efficiently evaluated [104] using the convolution theorem on the integral in Eq. (82), resulting in

(83)$${S(t,a) \propto {{\| {{{\rm FT}^{- 1}}[{{\textbf{j}}(f){W^*}(af)}]}\|}^2},}$$

where $W(f)$ denotes the Fourier transform of the mother wavelet. Eq. (83) is efficient, since ${\textbf{j}}(f)$ is independent of $a$, and $W(f)$ is known analytically. As discussed, the frequency is related to $a$ by $\omega = \Omega /a$, and $a$ is sampled dyadically.

We now provide a specific example of an HHG calculation. We consider a monolayer of hexagonal boron nitride (hBN), with the band structure calculated using the pseudopotential method detailed in [56,102]. A two-band model is considered, consisting of the highest valence band and the lowest conduction band, with the hexagonal BZ having minimum band gaps at the $K$ high-symmetry points with energy 7.8 eV. The periodic TPT structure gauge is constructed following the procedure outlined in Section 3.D, which provides us with smooth and periodic Berry connections and dipole couplings. Prior to the time propagation, the valence band is assumed fully occupied and the conduction band empty. We irradiate the monolayer with a laser linearly polarized along the armchair direction $\hat{\textbf{x}}$ [indicated by the red arrow in the inset of Fig. 2(a)], with the vector potential of the form

(84)$${\textbf{A}}(t) = {A_0}\mathop {\cos}\nolimits^2 \left[{\frac{{\pi t}}{{2\tau}}} \right]\sin ({\omega _0}t)\hat{\textbf{x}},\quad t \in [- \tau ,\tau],$$

where ${A_0} = 0.35$ ($I = 3.5 \;{\rm{TW}}/{\rm{cm}}^2$), $\omega = 0.0285$ ($\lambda = 1600 \;{\rm nm}$), and $\tau = 58.7\;{\rm fs}$. For the time propagation, we employ Eq. (74), i.e., the SBEs in the moving frame. A Monkhorst–Pack mesh is used with total number of ${\textbf{K}}$ discretization points $300 \times 300 = 9 \times {10^4}$ and the dephasing time is chosen ${T_2} = 5\;{\rm fs}$. Per the discussion in Section 4.D, if the simulations are done in the VG Bloch scheme, a much sparser ${\textbf{k}}$ grid of $100 \times 100$ is sufficient to obtain a spectrum with approximately the same level of accuracy, but using more bands. The HHG spectrum for harmonics polarized along the arm-chair direction $\hat{\textbf{x}}$ [blue arrow, inset of Fig. 3(a)] is calculated using Eq. (81) with $\hat{\textbf{n}} = \hat{\textbf{x}}$, and shown in Fig. 3(a). Harmonics up to the 28th order are observed, consisting of both even- and odd-order harmonics. The spectrum in Fig. 3(a) is further separated into the intraband and non-intraband contributions from the microscopic current, following Eq. (76). Below the band gap energy (dashed line), the odd harmonics are dominated by the intraband contribution, originating from the carrier transport in the individual bands. Above the band gap, however, the non-intraband contribution dominates, originating from the coupling between the bands. The even-order harmonics below the band gap are due to the anomalous current, which is part of the non-intraband contribution.

Fig. 3. (a) HHG spectrum for hBN irradiated by a 1600 nm, ${3.5}\;{\rm{TW/c}}{{\rm{m}}^2}$, 58.7 fs pulse. The inset shows a sketch of the crystal structure, the driving laser polarization direction (red arrow along the armchair direction) and the detected HHG polarization direction (blue arrow). The horizontal dashed line shows the position of the minimum band gap. (b) Time-frequency profiles for the harmonics. The labels A and B mark two different features in the spectrum. The results from the semiclassical recollision model (see Section 5) are superimposed, with the gray points originating from recollisions of electron-hole pairs created near an ${M_1}$ symmetry point (an $M$ point located on the line going through $\Gamma$ along the armchair direction), while the purple points are from electron-hole pairs created near the other $M$ points and $K$ points.

Download Full Size | PDF

The density plot in Fig. 3(b) shows the time-frequency profiles of the harmonic emissions, obtained using Eq. (83). The intraband harmonics below the band gap have a broad time profile, with the highest-order intraband harmonics [around harmonic 5 (H5)] emitted near $t = n/2$ optical cycles ($n$ an integer), corresponding to the zeros of the vector potential ${\textbf{A}}(t)$ in Eq. (84). This can be understood from the dependence of the intraband current on the time-dependent band structure $E_n^{{\textbf{K}} + {\textbf{A}}(t)}$ [Eq. (76a)], where the harmonic emissions occur at times corresponding to the largest band curvature (largest rate of change in the group velocity) [41]. In the case of hBN, this is near the $K$ and $M$ high-symmetry points.

The non-intraband harmonics above the band gap energy in Fig. 3(b) have a distinct emission profile compared to the intraband case. The strongest non-intraband harmonic at H13 [see also Fig. 3(a)] are emitted at around $t = n/2$ optical cycles, and we have labeled the structure in the time profiles as “A.” The time profiles for the harmonics emitted above 17th order have a bow-like structure, with the highest-order harmonics emitted at around ${t_1} = (0.25 + n/2)$ optical cycles, which we have labeled “B.” For $t \lt {t_1}$, the slope of the time-profile is positive, corresponding to a positive chirp, while for $t \gt {t_1}$ the emissions are negatively chirped. This bow-like structure is similar to the time-frequency profile for HHG in gases, where every energy below the cut-off harmonic is emitted twice, corresponding to the short and long trajectories [105]. Using shorter dephasing times ${T_2}$ will suppress the long trajectories and result in a more well-resolved HHG spectra. In Section 5, we will introduce the semiclassical recollision model for non-intraband harmonics, and discuss that the structure A and B originate from recollisions of electron-hole pairs created near the different symmetry point in the BZ.

5. SADDLE-POINT EQUATIONS AND THE RECOLLISION MODEL FOR HHG

In the previous section, we have presented methods to obtain relevant observables such as the microscopic current and the HHG spectrum, by time-propagating the relevant EOM. These black-box calculations can be considered numerical experiments that contain all the relevant information, but they are often too complex to gain physical insights, with everything intermingled. In this section, we discuss methods that can help us gain physical intuition and understanding, especially on the emission dynamics of HHG.

A. Saddle-Point Method for HHG

We consider a two-band model with nondegenerate bands that include an initially filled valence band denoted by the band index $v$ and an empty conduction band denoted by $c$. For concreteness, we work in the VG and in the Houston basis from Section 4.C. The EOM in Eq. (74) (also referred as the SBEs in the moving frame) reduce to

(85a)$${\dot {\bar \rho}} _{\textit{vv}}^{\textbf{K}}(t) = i{\textbf{F}}(t) \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(t)}}\bar \rho _{\textit{vc}}^{\textbf{K}}(t) + {\rm{c.c.}},$$(85b)$${\dot {\bar \rho}} _{\textit{cc}}^{\textbf{K}}(t) = - i{\textbf{F}} \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(t)}}\bar \rho _{\textit{vc}}^{\textbf{K}}(t) + {\rm{c.c.}},$$(85c)$$\begin{split}{\dot{\bar \rho}} _{\textit{cv}}^{\textbf{K}}(t)& =\left[{- i\omega _g^{{\textbf{K}} + {\textbf{A}}(t)} - i{\textbf{F}}(t) \cdot \Delta {{\cal A}^{{\textbf{K}} + {\textbf{A}}(t)}}} \right]\bar \rho _{\textit{cv}}^{\textbf{K}}(t)\\ &\quad- i\left[{\bar \rho _{\textit{vv}}^{\textbf{K}}(t) - \bar \rho _{\textit{cc}}^{\textbf{K}}(t)} \right]{\textbf{F}}(t) \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(t)}},\end{split}$$

with $\omega _g^{\textbf{K}} \equiv E_c^{\textbf{K}} - E_v^{\textbf{K}}$ the band gap and $\Delta {{\cal A}^{\textbf{K}}} \equiv {\cal A}_c^{\textbf{K}} - {\cal A}_v^{\textbf{K}}$ the Berry connection difference. For HHG in semiconductors and insulators, the population transfer to the conduction band is small, and we make the approximation $\bar \rho _{\textit{vv}}^{\textbf{K}} - \bar \rho _{\textit{cc}}^{\textbf{K}} \approx 1$. The formal solutions to Eq. (85) now read

(86a)$$\bar \rho _{\textit{vv}}^{\textbf{K}}(t) = i\int^t {\rm d}s{\textbf{F}}(s) \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(s)}}\bar \rho _{\textit{vc}}^{\textbf{K}}(s) + {\rm{c.c.}},$$(86b)$$\bar \rho _{\textit{cc}}^{\textbf{K}}(t) = - i\int^t {\rm d}s{\textbf{F}}(s) \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(s)}}\bar \rho _{\textit{vc}}^{\textbf{K}}(s) + {\rm{c.c.}},$$(86c)$$\begin{split}\bar \rho _{\textit{cv}}^{\textbf{K}}(t)& = - i\int^t {\rm d}s{\textbf{F}}(s) \cdot {{\textbf{d}}^{{\textbf{K}} + {\textbf{A}}(s)}}{e^{- T_2^{- 1}(t - s)}}\\&\quad \times {e^{- i\int_s^t \left[{\omega _g^{{\textbf{K}} + {\textbf{A}}(t^\prime)} + {\textbf{F}}(t^\prime) \cdot \Delta {{\cal A}^{{\textbf{K}} + {\textbf{A}}(t^\prime)}}} \right]{\rm d}t^\prime}},\end{split}$$

which can easily be checked by insertion. We are interested in the above-band gap harmonics, which are dominated by the non-intraband contribution, as illustrated by the example shown in Fig. 3. Insertion of Eq. (86) into Eq. (76b), and transforming into the fixed frame ${\textbf{k}} \equiv {\textbf{K}} + {\textbf{A}}(t)$ results in

(87)$${j_\mu ^{{\rm{ter}}}(t) = {N^{- 1}}\sum\limits_{\textbf{k}} R_\mu ^{\textbf{k}}\int^t {T^{\kappa (t,s)}}{e^{- i{S^\mu}({\textbf{k}},t,s)}}{\rm d}s + {\rm{c.c.}}},$$

with $\mu= \{x,y,z\}$ the Cartesian indices, ${T^{\kappa (t,s)}} = |{\textbf{F}}(s) \cdot {{\textbf{d}}^{\kappa (t,s)}}|$ the transition matrix element, $R_\mu ^{\textbf{k}} = \omega _g^{\textbf{k}}|d_\mu ^{\textbf{k}}|$ the recombination dipole, $\kappa (t,t^\prime) = {\textbf{k}} - {\textbf{A}}(t) + {\textbf{A}}(t^\prime)$ the time-dependent crystal momentum, and c.c. stands for complex conjugate. The times $s$ and $t$ can be interpreted as the excitation and emission times, respectively. The accumulated phase in Eq. (87) is

(88)$$\begin{split}{S^\mu}({\textbf{k}},t,s) = \int_s^t[{\omega _g^{\kappa (t,t^\prime)} + {\textbf{F}}(t^\prime) \cdot \Delta {{\cal A}^{\kappa (t,t^\prime)}}}]{\rm d}t^\prime+ {{\alpha ^{{\textbf{k}},\mu}} - {\beta ^{\kappa (t,s)}}},\end{split}$$

with ${\alpha ^{{\textbf{k}},\mu}} \equiv \arg (d_\mu ^{\textbf{k}})$ the transition-dipole phases [101,106], and ${\beta ^{\kappa (t,s)}} \equiv \arg [{\textbf{F}}(s) \cdot {{\textbf{d}}^{\kappa (t,s)}}]$. Note that inclusion of ${\alpha ^{{\textbf{k}},\mu}}$ and ${\beta ^{\kappa (t,s)}}$ in ${S^\mu}({\textbf{k}},t,s)$ results in it being structure-gauge invariant [101,107].

We are interested in the frequency-resolved non-intraband current, $j_\mu ^{{\rm{ter}}}(\omega) = \int_{- \infty}^\infty {\rm d}t{e^{i\omega t}}j_\mu ^{{\rm{ter}}}(t)$, which contributes to the non-intraband HHG spectrum. The saddle-point approximation in the absence of simple poles consists of only including the stationary (saddle) points of the phase factor involving ${S^\mu}({\textbf{k}},t,s) - \omega t$, since other contributions lead to highly oscillatory terms in the integrand of $j_\mu ^{{\rm{ter}}}(\omega)$. Taking the partial derivatives with respect to the three integration variables ${\textbf{k}}$, $s$, and $t$, the saddle point conditions are [101,107]

(89a)$$\omega _g^{\kappa (t,s)} + {\textbf{F}}(s) \cdot {{\cal Q}^{\kappa (t,s)}} = 0,$$(89b)$$\Delta {{\textbf{R}}^\mu} \equiv \Delta {\textbf{r}} - {{\cal D}^{{\textbf{k}},\mu}} + {{\cal Q}^{\kappa (t,s)}} = {\textbf{0}},$$(89c)$$\omega _g^{\textbf{k}} + {\textbf{F}}(t) \cdot \left[{{{\cal Q}^{\kappa (t,s)}} + \Delta {\textbf{r}}} \right] = \omega ,$$

with the electron-hole separation vector and group velocities

(90a)$$\Delta {\textbf{r}} \equiv \int_s^t \left[{{\textbf{v}}_c^{\kappa (t,t^\prime)} - {\textbf{v}}_v^{\kappa (t,t^\prime)}} \right]{\rm d}t^\prime ,$$(90b)$${\textbf{v}}_n^{\kappa (t,t^\prime)} \equiv {\nabla _{\textbf{k}}}E_n^{\kappa (t,t^\prime)} + {\textbf{F}}(t^\prime) \times \boldsymbol{\Omega} _n^{\kappa (t,t^\prime)},$$

and the structure-gauge invariant quantities

(91a)$${{\cal D}^{{\textbf{k}},\mu}} \equiv \Delta {{\cal A}^{\textbf{k}}} - {\nabla _{\textbf{k}}}{\alpha ^{{\textbf{k}},\mu}},$$(91b)$${{\cal Q}^{\textbf{k}}} \equiv \Delta {{\cal A}^{\textbf{k}}} - {\nabla _{\textbf{k}}}{\beta ^{\textbf{k}}}.$$

We denote a saddle point, i.e., a solution to Eqs. (89a)–(89c), as $\{\bar{\textbf{k}},\bar t,\bar s\}$. The essense of the interband HHG process is contained in the interpretation of Eqs. (89a)–(89c) in terms of the following three steps: an electron-hole pair is created by tunnel excitation at time $\bar s$ and with the crystal momentum ${{\textbf{k}}_0} \equiv \bar{\textbf{k}} - {\textbf{A}}(\bar t) - {\textbf{A}}(\bar s)$; the hole and electron are accelerated by the laser with the instantaneous group velocities $v_v^{\bar{\textbf{k}} - {\textbf{A}}(\bar t) + {\textbf{A}}(t^\prime)}$ and $v_c^{\bar{\textbf{k}} - {\textbf{A}}(\bar t) + {\textbf{A}}(t^\prime)}$, respectively; and the electron-hole pair recombine at time $\bar t$ with final crystal momentum $\bar{\textbf{k}}$ and relative distance $\Delta {\textbf{r}}$, with the simultaneous emission of high-harmonics with energy $\omega$.

The saddle-point conditions cannot generally be satisfied by purely real values of ${\textbf{k}}$, $t$ and $s$, and must generally be solved by analytic continuation of the parameter space into the complex plane. For example, in band-gap materials where ${{\cal D}^{{\textbf{k}},\mu}}$ and ${{\cal Q}^{\textbf{k}}}$ are negligible, Eq. (89a) requires the band gap to be zero, which clearly cannot be satisfied for real values of ${\textbf{k}}$; similarly, Eq. (89b) requires $\Delta {\textbf{r}} = {\textbf{0}}$, which is often too restrictive for real-valued saddle points. Thus, a fully rigorous quantum solution would involve complex-valued saddle points and constitute a monumental numerical task, and we mark this as another interesting direction of research for HHG in solids. Recent progress has been made toward approximately solving the saddle-point equations for reduced-dimensionality model systems in [108,109]. It remains to be seen whether such formalisms can treat electron-hole pairs starting from different regions in the BZ. In the next subsection, we will describe a semiclassical solution to the saddle-point conditions that can provide temporal and spectral insights into the mechanism of interband HHG.

The saddle-point approximation to the frequency-resolved non-intraband current in Eq. (87) involves deforming the integration contours from the real axes into the complex planes, going through the complex saddle points along the path of steepest descent [110] defined by constant ${\rm{Re}}[{{S^\mu}({\textbf{k}},t,s) - \omega t}]$. The approximation reads in the present case [32,109]

(92)$${j_\mu ^{{\rm{ter}}}(\omega) \propto \sum\limits_{\bar{\textbf{k}},\bar t,\bar s} \frac{{R_\mu ^{\bar{\textbf{k}}}{T^{\,\bar{\textbf{k}} - {\textbf{A}}(\bar t) + {\textbf{A}}(\bar s)}}{e^{- i\left[{{S^\mu}(\bar{\textbf{k}},\bar t,\bar s) - \omega \bar t} \right]}}}}{{\sqrt {\det [{\partial ^2}{S^\mu}(\bar{\textbf{k}},\bar t,\bar s)]}}} + {\rm{c}}{\rm{.c}}.(\omega \to - \omega)},$$

where we have used the notation ${\partial ^2}{S^\mu}$ for the Hessian and the second term denotes complex conjugate of first term and with $\omega \to - \omega$. In systems where ${{\cal D}^{{\textbf{k}},\mu}}$ and ${{\cal Q}^{\textbf{k}}}$ can be chosen to be zero, the Hessian is often proportional to $\parallel {\nabla _{\textbf{k}}}\omega _g^{\bar{\textbf{k}}}\parallel$ [32]. This makes it clear that for small values of $\parallel {\nabla _{\textbf{k}}}\omega _g^{\bar{\textbf{k}}}\parallel$, i.e., when the valence and conduction bands have similar slopes at the saddle point $\bar{\textbf{k}}$, the harmonic yield at the corresponding harmonic energy $\omega = \omega _g^{\bar{\textbf{k}}}$ is expected to be greatly enhanced, a phenomenon termed as “spectral singularities” in [32].

B. Semiclassical Recollision Model

As mentioned, the full quantum solution to the saddle-point equations in Eq. (89) presents a monumental task, and a solution has only been attempted for very simple model systems. A lot of physical insight and intuition, however, can be gained from a semiclassical solution [14,38,101,107,111,112] to the saddle-point equations, which we present below and sketched in Fig. 1.

In the semiclassical procedure, an initial tunneling time $s$ and an initial crystal momentum ${{\textbf{k}}_0}$ is picked, where equation (89a) determines ${{\textbf{k}}_0} \equiv \kappa (t,s)$ at tunneling time $s$. Since this equation generally cannot be satisfied for real-valued ${{\textbf{k}}_0}$ [the second term of Eq. (89a) is generally small], we choose ${{\textbf{k}}_0}$ to be at or close to a high-symmetry point with a small band gap, which corresponds to a high tunneling probability. The tunneling time is picked in an optical cycle of the pulse $s \in [- T,0]$.

The next step in the semiclassical solution involves the propagation of the integral in Eq. (90a), i.e., calculating the electron and hole classical real-space motions for $t^\prime \in [s,t]$ with the time-dependent crystal momentum $\kappa (t,t^\prime) = {{\textbf{k}}_0} \,+ {\textbf{A}}(t^\prime) - {\textbf{A}}(s)$ and the group velocities given in Eq. (90b). Equation (89b) denotes the electron-hole recollision condition, i.e., when the electron-hole distance $\left\| {\Delta {\textbf{r}}} \right\|$ is equal to $\left\| {{{\cal D}^{{\textbf{k}},\mu}} - {{\cal Q}^{{{\textbf{k}}_0}}}} \right\|$, where the latter quantity is generally small in the systems considered by us. For real-valued saddle points, this condition is often too restrictive due to the complicated band dispersions in a crystal. This is in contrast to the free-electron dispersion relevant for HHG in gases. We thus relax the recollision condition. At each $t^\prime $ during the time-propagation, we calculate $\Delta {{\textbf{R}}^\mu}$, and record a semiclassical recollision event if: (i) $\left\| {\Delta {{\textbf{R}}^\mu}} \right\|$ as a function of $t^\prime $ is a local minimum and (ii) $\left\| {\Delta {{\textbf{R}}^\mu}} \right\| \lt {R_0}$ is fulfilled, with ${R_0}$ a preset recollision threshold value. The threshold ${R_0}$ is chosen as the minimum ${R_0}$ such that the recollision-model results agree with the time profiles obtained from the quantum simulation results. Often, ${R_0}$ can be multiple times greater than the lattice constants of the crystal. The last saddle point in Eq. (89c) determines the electron-hole recollision energy $\omega$, and thus the energy of the emitted harmonics.

As discussed above, the relaxation of the electron-hole recollision condition means that we allow for imperfect recollisions where $\left\| {\Delta {{\textbf{R}}^\mu}} \right\| \ne {\textbf{0}}$ [107,112–114]. One consequence of such an recollision event is that the harmonic energy $\omega$ attains an extra electron-hole-pair polarization energy ${\textbf{F}}(\bar t) \cdot \Delta {\textbf{r}}$ [112]. Physically, this energy constitutes the potential energy of the electric dipole comprised of the electron-hole pair at the time of recollision. Clearly, imperfect recollisions will occur whenever the different Cartesian components of $\Delta {{\textbf{R}}^\mu}$ are nonzero or zero at different times, i.e., whenever the direction of motion of the time-dependent crystal momentum ${\partial _{{t^\prime}}}\kappa (t,t^\prime) = - {\textbf{F}}(t^\prime)$ is not along the instantaneous group velocities ${\textbf{v}}_n^{\kappa (t,t^\prime)}$. This condition can, e.g., be satisfied for systems driven by elliptically polarized pulses or in systems with large Berry curvatures. The nonzero recollision distance is a consequence of the spatially extended electron and hole wavepackets, which in periodic systems can span over several unit cells [112,113,115]. At the time of recollision, even if the electron-hole centers are displaced, the wave packets can still overlap and ensure a recollision event. Indeed, the minimum value of the chosen recollision threshold ${R_0}$ at which the semiclassical recollision model and quantum emission profiles agree is a qualitative measure of the width of the wave packets. Finally, even though the presented semiclassical model assumes the birth of the electron-hole pair with zero spatial displacement, recent work [108,116] has indicated that the electron and hole can emerge displaced in real-space after tunneling.

In the monolayer–hBN–HHG example presented in Section 4.E, the semiclassical result for the recollision energy versus recollision time is plotted in Fig. 3(b) as points superimposed on the time-frequency profiles. The recollision events represented by the dark gray points originate from electron-hole pairs initially created in a disk of radius $\Delta k = 0.1$ around a ${M_1}$ symmetry points, while the purple points originate from electron-hole pairs initially created near the other $M$ points and $K$ points. Note that an ${M_1}$ point is defined as an $M$ point located on the line going through the $\Gamma$ point along the armchair direction. The agreement of the semiclassical result with the time-frequency profile is evident, and more physical understanding is gained as the semiclassical model is able to attribute the different structures in the density plot (denoted as structures A and B in Section 4.E) to different tunneling sites in reciprocal space. Since many regions in the BZ are relevant for HHG in solids, the time-frequency profiles and subcycle emission dynamics can be quite complicated for different materials and laser polarizations [32,107,111,114,116]. At the same time, however, the fact that electron-hole pairs created in different BZ regions leads to distinct harmonic time-frequency characteristics can potentially facilitate the all-optical reconstruction of the band structure in the whole BZ, and not only near the minimum band gap as demonstrated in [31]. We note that although the expanded recollision model discussed here is able to qualitatively explain the HHG emission mechanism, there are differences between the semiclassical and the quantum results, as evidenced in Fig. 3(b) by the deviations in the cutoff energies and emission times for structure B. Recent research has discussed such shortcomings in the semiclassical models and potential improvements [108,112,117,118].

Note that a natural procedure to obtain a full solution to the saddle-point equations (89) would be to first pick a harmonic frequency $\omega$, and then find the desirable saddle points by some numerical procedure. In contrast, in the approximate semiclassical solution, we pick an initial tunneling time $s$ and crystal momentum ${{\textbf{k}}_0}$, and check the above-described recollision conditions on-the-fly during the classical time propagation of a trajectory. If the recollision conditions are satisfied, we record it as a recollision event and say that we have found a saddle-point solution $\{\bar{\textbf{k}},\bar t,\bar s\}$. Due to the deterministic nature of the classical trajectories, it is numerically manageable to pick all combinations of the tunneling time in an optical cycle $s \in [- T,0]$ and the crystal momentum ${{\textbf{k}}_0}$ in a sphere around the high-symmetry points, as well as propagate the trajectories up to two optical cycles after tunneling $t \in [s,s + 2T]$.

6. MACROSCOPIC PROPAGATION

Up until this point, we have mostly discussed ways to treat the microscopic dynamics responsible for HHG in solids. However, when electromagnetic radiation propagates through a medium, the medium responds by not only generating new radiation, as discussed above, but sometimes also by modifying the propagating radiation. In the field-intensity regime of pulsed lasers, the propagation of the electromagnetic fields can be described classically, i.e., governed by Maxwell’s equations,

(93a)$$\nabla \cdot {\textbf{F}} = \epsilon _0^{- 1}\rho ,$$(93b)$$\nabla \cdot {\textbf{B}} = 0,$$(93c)$$\nabla \times {\textbf{F}} = - {\partial _t}{\textbf{B}},$$(93d)$$\nabla \times {\textbf{B}} = {\mu _0}\left({{\textbf{J}} + {\epsilon _0}{\partial _t}{\textbf{F}}} \right),$$

with ${\epsilon _0}$ the vacuum permittivity, ${\mu _0}$ the vacuum permeability, $\rho$ the (induced) charge density, and ${\textbf{J}}$ the current density proportional to the microscopic current ${\textbf{j}}$ from Section 4. We thus consider here the case where there are no free charges or currents in the medium. Generally, all sources and fields in Eq. (93) are dependent on space ${\textbf{r}}$ and time $t$, which we have omitted for notational clarity. For small intensities, the medium response is linear with respect to the fields, and governed by the linear constitutive equations involving linear susceptibilities. At higher intensities, the nonlinear response can be modeled by a series expansion in terms of the nonlinear susceptibilities. When such a series expansion fails, as is often the case for the intensity regimes responsible for the highly nonlinear recollision mechanisms and the non-intraband currents, one must resort to numerical solutions of the microscopic response, which was the topic in the previous sections.

In terms of the vector potential ${\textbf{A}}$ and scalar potential $\Phi$ defined in Eq. (49), the inhomogeneous wave equations can easily be obtained from Eqs. (93a) and (93c) by taking the curl of ${\textbf{B}}$,

(94)$$\begin{split}{{\nabla ^2}\Phi + {\partial _t}\big({\nabla \cdot {\textbf{A}}} \big) = - \epsilon _0^{- 1}\rho}\\{{c^{- 2}}\partial _t^2{\textbf{A}} - {\nabla ^2}{\textbf{A}} + \nabla ({c^{- 2}}{\partial _t}\Phi + \nabla \cdot {\textbf{A}}) = {\mu _0}{\textbf{J}}.}\end{split}$$

In principle, the source terms can be obtained from the microscopic calculations at each time and position, and the final coupled Maxwell–Schrödinger equations should then be solved. However, due to the insurmountable computational complexities involved, appropriate approximations are usually required. One often-used approximation is to define two spatial scales with two numerical grids: one macroscopic scale with coordinate ${\textbf{R}}$ on the order of the electromagnetic wavelength $\lambda$ to treat the pulse propagation, and a smaller scale ${\textbf{r}}$ on the order of unit cells ($\ll \lambda$), where a local dipole approximation for the microscopic physics can be made. One advantage of using the scalar and vector potentials instead of the physical fields is the flexibility to choose different laser gauges in the two different scales. One of the gauges in Section 4, for example, can be used for the microscopic physics, while a gauge with $\Phi \equiv 0$ can be used for the macroscopic scale [72]. To further reduce the problem, an approximate 1D propagation scheme can be used, where the laser pulse propagates along $Z$ (perpendicular to a crystal surface), such that the vector potential in the local dipole approximation reads ${{\textbf{A}}_{\textbf{R}}}(t) = {A_Z}(t)\hat{\textbf{e}}$, with $\hat{\textbf{e}}$ the laser polarization direction. The symmetry of the setup is assumed to be such that the generated current has the same polarization as the laser, and ${{\textbf{J}}_{\textbf{R}}}(t) = {J_Z}(t)\hat{\textbf{e}}$. The wave equation reduces in this case to

(95)$${c^{- 2}}\partial _t^2{A_Z} - \partial _Z^2{A_Z} = {\mu _0}{J_Z},$$

with ${J_Z}$ the current obtained from the microscopic calculation. Such an approach has, for example, been used in [55] to show that the HHG spectrum after propagation through a bulk crystal exhibits a much “cleaner” spectrum (more well-resolved harmonic peaks) compared to the purely microscopic result, which is attributed to the destructive interference between electron-hole recollision events at different recombination times along the propagation path. Very recently, the optimal thickness for HHG in silicon thin films has been investigated [119]. Note, however, that this scheme does not include the radial variation of the laser and thus does not include many nonlinear optical effects such as self-focusing. In addition, the second-order derivative of the propagation coordinate $Z$ in Eq. (95) requires a dense discretization along $Z$, with the maximum realizable propagation distance being a couple of micrometers [55,72,119].

A different approach is to consider the wave equation for the physical fields

(96a)$${c^{- 2}}\partial _t^2{\textbf{F}} - {\nabla ^2}{\textbf{F}} = - \varepsilon _0^{- 1}\nabla \rho - {\mu _0}{\partial _t}{\textbf{J}},$$(96b)$${c^{- 2}}\partial _t^2{\textbf{B}} - {\nabla ^2}{\textbf{B}} = {\mu _0}(\nabla \times {\textbf{J}}).$$

For electric fields slowly varying in the transverse dimensions, the scalar approximation for the wave in Eq. (96) can be made [68]. For nonrelativistic carriers, the magnetic field ${\textbf{B}}$ can be neglected, and the wave equation reduces to

(97)$${c^{- 2}}\partial _t^2F - {\nabla ^2}F = - {\mu _0}{\partial _t}J,$$

where the right-hand side of Eq. (97) again represents the source terms generated by the medium in response to the propagating field. This equation has a second-order derivative along the propagation direction $z$, which requires significant numerical effort. The computational effort can be significantly reduced by first transforming to a frame that moves at the speed of light, and then ignoring the second derivative with respect to the new $z$ compared to $\frac{2}{c}{\partial _t}{\partial _z}$ and the transverse derivative

(98)$${\left[{\nabla _ \bot ^2 - \frac{2}{c}{\partial _t}{\partial _z}} \right]F = {\mu _0}{\partial _t}{J_0}.}$$

This approximation is usually termed the slowly evolving wave approximation [68–71,120] and also implies ignoring backward-propagating waves. Eq. (98) can be conveniently solved in the spectral domain.

Many similar variations of the envelope propagation in Eq. (98) exist in the literature, where the current or polarization often is separated into linear and nonlinear parts $J = {J^{(1)}} + {J^{\rm{NL}}}$ and the linear part is expressed in terms of the linear susceptibilities and permittivities. One example is the unidirectional pulse propagation equation (UPPE), which assumes that the nonlinear response is purely due to the forward-propagating wave [69,121]. Recently, the UPPE was applied to HHG in solids to show that for propagation lengths longer than the laser wavelength, the propagation significantly reduces the HHG yield and can potentially be responsible for the short dephasing times ${T_2}$ used in microscopic simulations to match with experiments [90]. The UPPE has also in its full vectorial form been used to show that the term in Eq. (96) involving $\epsilon _0^{- 1}\rho = \nabla \cdot {\textbf{F}}$ in some cases can have a pronounced effect on the self-focusing of ultrashort pulses [69].

A. Spatio-Spectral Properties of Solid-State HHG

In this section, we give an example of a spatio-spectral analysis of solid-state HHG. The conceptual sketch of our numerical setup is shown in Fig. 4: An incoming laser beam propagates normal to a very thin crystal sample along the $z$ axis, hits the sample, and generates high-order harmonics; the beam and the generated harmonics propagate from the sample (near field) toward an opaque screen located a great distance away from the sample (far field). The opaque screen has a circular aparture that filters away the parts of the beam that has a large spatial divergence. Optionally, beyond the circular aperture, a focusing lens can be installed to focus the filtered spectrum back to conditions that mimic the near field. The spatially resolved HHG spectrum and time profiles often provide a separation of different contributions to the HHG process as we will discuss below.

Fig. 4. Sketch showing a numerical or (potentially) experimental setup to probe the spatio-spectral properties of solid-state HHG in which an incoming laser beam hits a thin sample and produces HHG in the near field. The fields are then propagated from the near field to the far field, where a spatial filter, for example, can select the on-axis radiation, essentially filtering out contributions from undesired electron-hole recollisions in the near-field generation process (see text).

Download Full Size | PDF

For the incident beam, we assume a Gaussian beam ${\textbf{F}}({\textbf{r}}) = {F_0}(r,z){e^{\textit{ikz}}}\hat{\textbf{e}}$:

(99)$${F_0}(r,z) = C\frac{{{w_0}}}{{w(z)}}{e^{- \frac{{{r^2}}}{{w{{(z)}^2}}}}}{e^{i\left[{\frac{{k{r^2}}}{{2R(z)}} - \mathop {\tan}\nolimits^{- 1} \left({\frac{z}{{{z_0}}}} \right)} \right]}},$$

where $z$ is the propagation direction, $r$ is the radial coordinate measured from the beam axis, $k = 2\pi n/\lambda$ is the wave number, $w(z) = {w_0}\sqrt {1 + {z^2}/z_0^2}$ is the beam waist, ${z_0} = \pi w_0^2n/\lambda$ is the Rayleigh range, $R(z) = z + z_0^2/z$ is the radius of curvature, and $C$ is a constant that determines the pulse amplitude.

We take the crystal sample to be a monolayer of hBN (see Section 4.E), placed at the Gaussian beam focus $z = 0$. The pulse parameters are chosen to be the same as in Section 4.E, i.e., 1600 nm, ${3.5}\;{\rm{TW/c}}{{\rm{m}}^2}$ peak intensity, and 58.7 fs pulse duration. At the sample, the radial intensity profile of the beam is ${I_0}(r,0) = {| C |^2}{e^{- 2{r^2}/w_0^2}}$, with ${w_0} = 20$ µm, and we use it to calculate the microscopic response at each $r$, resulting in the frequency- and radial-dependent near field ${F_{\rm{near}}}(\omega ,r) \propto \omega J(\omega ,r)$. The microscopic dynamics is solved using the VG EOM covered in Section 4.A (no dephasing), including one valence and nine conduction bands, and the pulse parameters are chosen to be the same as in Fig. 3. The HHG spectrum, ${| {{F_{\rm{near}}}(\omega ,r)} |^2}$, is plotted in Fig. 5(b). For harmonics in the interval H10-H20, the spectrum shows complicated behavior compared to the interval H20–H30, reflecting the fact that the former interval has contributions from electron-hole recombinations originating from different BZ symmetry points, as well as the interference of the short and long trajectories, as shown in Fig. 3 and the accompanying main text. The total near-field spectrum is given by ${S_{\rm{near}}}(\omega) \propto \int_0^\infty {| {{F_{\rm{near}}}(\omega ,r)} |^2}r{\rm d}r$, which is shown in Fig. 5(a) by the black curve. Again, all harmonics in the interval H20–H30 are well-resolved, while only some of the odd harmonics in the interval H10-20 are visible. Note that since we are dealing with a monolayer material, the macroscopic propagation in the sample is not required.

Fig. 5. Near- and far-field HHG spectra, for a monolayer of hBN irradiated by a Gaussian beam with waist ${w_0} = 20$ µm and pulse parameters 1600 nm, ${3.5}\;{\rm{TW/c}}{{\rm{m}}^2}$, 58.7 fs. (a) Near-field, far-field, and filtered far-field radially integrated HHG spectra. (b) Radially resolved near-field HHG spectrum. (c) Radially resolved far-field HHG spectrum. The vertical dashed line marks the minimum band gap. The gray dotted curve plots the divergence radius in the far field expected for Gaussian beams with different carrier frequencies.

Download Full Size | PDF

In an actual experiment, the photodetector is placed far from the sample, and thus the spectra are only detected in the far field. As we will see, the far-field spectrum carries information on the underlying HHG process at the near field. For a Gaussian beam profile known at $z = 0$, using the paraxial wave equation, the beam profile at $z \gt 0$ is given as a Hankel transform [122–124]:

(100)$$\begin{split}\bar F(\omega ,r,z) &= - \frac{{ik{e^{i\frac{{k{r^2}}}{{2z}}}}{e^{\textit{ikz}}}}}{z}\int_0^\infty {F_{\rm{near}}}(\omega ,r)\\&\quad\times{e^{i\frac{{k{{r^\prime}^2}}}{{2z}}}}{J_0}\left({\frac{{krr^\prime}}{z}} \right)r^\prime {\rm d}r^\prime ,\end{split}$$

with ${J_0}$ the zeroth order Bessel function and the integration is performed over the near-field radial coordinate $r^\prime $. Equation (100) is the same formula as for the diffraction from a circular aparture in the Fresnel approximation. We calculate the far-field spectrum as ${F_{\rm{far}}}(\omega ,r) \equiv \bar F(\omega ,r,z = L)$, with $L = 1 \;{\rm m}$. The total integrated spectrum is accordingly ${S_{\rm{far}}}(\omega) \propto \int_0^\infty {| {{F_{\rm{far}}}(\omega ,r)} |^2}r{\rm d}r$.

Figure 5(c) shows the $r$ dependent far-field spectrum $|{F_{\rm{far}}}(\omega ,r{)|^2}$. As observed by comparing the $r$ axis between Figs. 5(a) and 5(b), the beam has diverged appreciably. For the low-order harmonics H1–H9, the divergence decreases with increasing harmonic order, following the divergence angle for a Gaussian beam $\theta = 2c/(n\omega {w_0})$ [122], which is plotted as the gray dotted line. Starting at around the band gap energy (${\sim}{\rm H}10$), the divergence angle of the harmonic radiation has a local maximum, and then decreases again for larger frequencies. This behavior indicates that the harmonics below and above the band gap are dominated by two distinct microscopic HHG mechanisms. This is in agreement with our discussion in Section 4.E, which showed that the harmonics below the band gap are dominated by the intraband current, while the harmonics above the band gap originate primarily from the non-intraband current. The spatiotemporal profile of the HHG clearly encodes this information. Due to energy conservation, the total far-field harmonic spectrum ${S_{\rm{far}}}(\omega)$ is identical to the near-field spectrum ${S_{\rm{near}}}(\omega)$, as shown by comparing the green dashed line and the black line in Fig. 5(a).

For the interband harmonics, the radial divergence is related to the radial dependence of the accumulated phase of an electron-hole pair between tunneling time $\bar s$ and recombination time $\bar t$ in the near field, i.e., ${S^\mu}(\bar{\textbf{k}},\bar t,\bar s)$ from Eqs. (87) and (88). The steeper the radial phase function (defined as the accumulated phase as a function of $r$) for a harmonic, the larger its radial profile becomes in the far field [124]. Consequently, the harmonics detected near the beam axis $r = 0$ in the far field mostly originate from recollided electron-hole pairs in the near field with small travel times between tunneling and recollision. We can isolate this near-axis spectrum and related trajectories by placing a filter of radius 1 cm in the far field, as shown in Fig. 4, and the resulting spectrum ${S_{\rm{filter}}}(\omega) \propto \int_0^{{r_{\rm{filter}}}} {| {{F_{\rm{far}}}(r,\omega)} |^2}r{\rm d}r$ is shown in Fig. 5(a) as the red dotted line. In the filtered spectrum, more harmonics are discernible; e.g., H14 and H16 have become visible.

Clearly, the addition of the radial degree of freedom provides additional understanding of the underlying recollision dynamics of solid-state HHG, and can also potentially provide realistic experimental pathways to probe the HHG process. For example, a far-field spectrum was experimentally measured in [17] to confirm the spatial coherence of the generated extreme ultraviolet high-harmonic radiation in ${\rm{Si}}{{\rm{O}}_2}$. In [88], interferometry of the dipole phase in HHG was performed by experimentally measuring the far-field spectrum of two overlapping beams in the near field.

7. SUMMARY AND OUTLOOK

In this tutorial, we have given a hands-on introduction to the theory of HHG in solids. In Section 2, we discussed the adiabatic states, the Berry phase, and related concepts. In Section 3, we described the time-independent problem of a crystalline solid and methods for structure gauge constructions for both nondegerate and degenerate cases. In Section 4, we covered approaches to describe the microscopic HHG mechanism in different laser gauges and structure gauges, and pointed out the advantages and drawbacks for the different methods, as well as provided an example HHG analysis for a monolayer material. In Section 5, we discussed the saddle-point approximation to HHG. The semiclassical solutions to the saddle-point equations, also termed the recollision model, could reveal spectro-temporal information about the HHG process. In Section 6, we formulated ways to describe the macroscopic HHG process, which involved the coupling of the microscopic dynamics to Maxwell’s equations. We provided an example of a monolayer irradiated by a Gaussian beam and discussed the spatio-spectral properties of the HHG detected in the far field.

Solid-state HHG is a rapidly expanding field, and many emerging theory trends are emerging that are beyond the scope of this tutorial. These include HHG in topological insulators [34,48,125], doped and amorphous systems [126,127], and strongly correlated systems [47,128–130]. We hope that this work gives newcomers to get an overview of the topic, as well as provides the necessary tools to perform simulations and stimulate the development of new theories.

Funding

National Science Foundation (PHY-1713671, PHY-2110317); Air Force Office of Scientific Research (FA9550-16-1-0013; supported parts of the development described in Sections 3 and 4).

Acknowledgment

Portions of this research were conducted with high-performance computational resources provided by the Louisiana Optical Network Infrastructure (http://www.loni.org). Author Lun Yue thanks Shicheng Jiang and Francois Mauger for useful discussions.

Disclosures

The authors declare no conflicts of interest.

Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. P. B. Corkum and F. Krausz, “Attosecond science,” Nat. Phys. 3, 381–387 (2007). [CrossRef]

2. F. Krausz and M. Ivanov, “Attosecond physics,” Rev. Mod. Phys. 81, 163–234 (2009). [CrossRef]

3. A. Flettner, T. Pfeifer, D. Walter, C. Winterfeldt, C. Spielmann, and G. Gerber, “High-harmonic generation and plasma radiation from water microdroplets,” Appl. Phys. B 77, 747–751 (2003). [CrossRef]

4. H. G. Kurz, D. S. Steingrube, D. Ristau, M. Lein, U. Morgner, and M. Kovačev, “High-order-harmonic generation from dense water microdroplets,” Phys. Rev. A 87, 063811 (2013). [CrossRef]

5. T. T. Luu, Z. Yin, A. Jain, T. Gaumnitz, Y. Pertot, J. Ma, and H. J. Wörner, “Extreme–ultraviolet high–harmonic generation in liquids,” Nat. Commun. 9, 3723 (2018). [CrossRef]

6. J. Inga, H. Martin, R. Dominik, P. Michael, J. Denis, P. Conaill, A. von Conta, S. Axel, and W. H. Jakob, “Attosecond spectroscopy of liquid water,” Science 369, 974–979 (2020). [CrossRef]

7. N. H. Burnett, H. A. Baldis, M. C. Richardson, and G. D. Enright, “Harmonic generation in CO₂ laser target interaction,” Appl. Phys. Lett. 31, 172–174 (1977). [CrossRef]

8. R. L. Carman, D. W. Forslund, and J. M. Kindel, “Visible harmonic emission as a way of measuring profile steepening,” Phys. Rev. Lett. 46, 29–32 (1981). [CrossRef]

9. D. von der Linde, T. Engers, G. Jenke, P. Agostini, G. Grillon, E. Nibbering, A. Mysyrowicz, and A. Antonetti, “Generation of high-order harmonics from solid surfaces by intense femtosecond laser pulses,” Phys. Rev. A 52, R25–R27 (1995). [CrossRef]

10. P. A. Norreys, M. Zepf, S. Moustaizis, A. P. Fews, J. Zhang, P. Lee, M. Bakarezos, C. N. Danson, A. Dyson, P. Gibbon, P. Loukakos, D. Neely, F. N. Walsh, J. S. Wark, and A. E. Dangor, “Efficient extreme UV harmonics generated from picosecond laser pulse interactions with solid targets,” Phys. Rev. Lett. 76, 1832–1835 (1996). [CrossRef]

11. A. H. Chin, O. G. Calderón, and J. Kono, “Extreme midinfrared nonlinear optics in semiconductors,” Phys. Rev. Lett. 86, 3292–3295 (2001). [CrossRef]

12. S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, “Observation of high-order harmonic generation in a bulk crystal,” Nat. Phys. 7, 138–141 (2011). [CrossRef]

13. O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, “Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations,” Nat. Photonics 8, 119–123 (2014). [CrossRef]

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

15. Y. S. You, D. A. Reis, and S. Ghimire, “Anisotropic high-harmonic generation in bulk crystals,” Nat. Phys. 13, 345–349 (2017). [CrossRef]

16. S. Jiang, S. Gholam-Mirzaei, E. Crites, J. E. Beetar, M. Singh, R. Lu, M. Chini, and C. D. Lin, “Crystal symmetry and polarization of high-order harmonics in ZnO,” J. Phys. B 52, 225601 (2019). [CrossRef]

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

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

19. G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne, K. J. Schafer, M. B. Gaarde, and D. A. Reis, “Solid-state harmonics beyond the atomic limit,” Nature 534, 520–523 (2016). [CrossRef]

20. N. Yoshikawa, T. Tamaya, and K. Tanaka, “High-harmonic generation in graphene enhanced by elliptically polarized light excitation,” Science 356, 736–738 (2017). [CrossRef]

21. H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, “High-harmonic generation from an atomically thin semiconductor,” Nat. Phys. 13, 262–265 (2017). [CrossRef]

22. H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, “Extremely efficient terahertz high-harmonic generation in graphene by hot Dirac fermions,” Nature 561, 507–511 (2018). [CrossRef]

23. N. Yoshikawa, K. Nagai, K. Uchida, Y. Takaguchi, S. Sasaki, Y. Miyata, and K. Tanaka, “Interband resonant high-harmonic generation by valley polarized electron-hole pairs,” Nat. Commun. 10, 3709 (2019). [CrossRef]

24. S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, “High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure,” Nat. Commun. 7, 13105 (2016). [CrossRef]

25. G. Vampa, B. G. Ghamsari, S. S. Mousavi, T. J. Hammond, A. Olivieri, E. Lisicka-Skrek, A. Y. Naumov, D. M. Villeneuve, A. Staudte, P. Berini, and P. B. Corkum, “Plasmon-enhanced high-harmonic generation from silicon,” Nat. Phys. 13, 659–662 (2017). [CrossRef]

26. M. Sivis, M. Taucer, G. Vampa, K. Johnston, A. Staudte, A. Y. Naumov, D. M. Villeneuve, C. Ropers, and P. B. Corkum, “Tailored semiconductors for high-harmonic optoelectronics,” Science 357, 303–306 (2017). [CrossRef]

27. Y. S. You, Y. Yin, Y. Wu, A. Chew, X. Ren, F. Zhuang, S. Gholam-Mirzaei, M. Chini, Z. Chang, and S. Ghimire, “High-harmonic generation in amorphous solids,” Nat. Commun. 8, 724 (2017). [CrossRef]

28. V. E. Nefedova, S. Fröhlich, F. Navarrete, N. Tancogne-Dejean, D. Franz, A. Hamdou, S. Kaassamani, D. Gauthier, R. Nicolas, G. Jargot, M. Hanna, P. Georges, M. F. Ciappina, U. Thumm, W. Boutu, and H. Merdji, “Enhanced extreme ultraviolet high-harmonic generation from chromium-doped magnesium oxide,” Appl. Phys. Lett. 118, 201103 (2021). [CrossRef]

29. Y. Bai, F. Fei, S. Wang, N. Li, X. Li, F. Song, R. Li, Z. Xu, and P. Liu, “High-harmonic generation from topological surface states,” Nat. Phys. 17, 311–315 (2021). [CrossRef]

30. C. P. Schmid, L. Weigl, P. Grössing, V. Junk, C. Gorini, S. Schlauderer, S. Ito, M. Meierhofer, N. Hofmann, D. Afanasiev, J. Crewse, K. A. Kokh, O. E. Tereshchenko, J. Güdde, F. Evers, J. Wilhelm, K. Richter, U. Höfer, and R. Huber, “Tunable non-integer high-harmonic generation in a topological insulator,” Nature 593, 385–390 (2021). [CrossRef]

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

32. A. J. Uzan, G. Orenstein, Á. Jiménez-Galán, C. McDonald, R. E. F. Silva, B. D. Bruner, N. D. Klimkin, V. Blanchet, T. Arusi-Parpar, M. Krüger, A. N. Rubtsov, O. Smirnova, M. Ivanov, B. Yan, T. Brabec, and N. Dudovich, “Attosecond spectral singularities in solid-state high-harmonic generation,” Nat. Photonics 14, 183–187 (2020). [CrossRef]

33. T. T. Luu and H. J. Wörner, “Measurement of the berry curvature of solids using high-harmonic spectroscopy,” Nat. Commun. 9, 916 (2018). [CrossRef]

34. A. Chacón, D. Kim, W. Zhu, S. P. Kelly, A. Dauphin, E. Pisanty, A. S. Maxwell, A. Picón, M. F. Ciappina, D. E. Kim, C. Ticknor, A. Saxena, and M. Lewenstein, “Circular dichroism in higher-order harmonic generation: Heralding topological phases and transitions in Chern insulators,” Phys. Rev. B 102, 134115 (2020). [CrossRef]

35. F. Bloch, “Über die quantenmechanik der elektronen in kristallgittern,” Z. Phys. 52, 555–600 (1929). [CrossRef]

36. W. V. Houston, “Acceleration of electrons in a crystal lattice,” Phys. Rev. 57, 184–186 (1940). [CrossRef]

37. S. Ghimire, A. D. DiChiara, E. Sistrunk, G. Ndabashimiye, U. B. Szafruga, A. Mohammad, P. Agostini, L. F. DiMauro, and D. A. Reis, “Generation and propagation of high-order harmonics in crystals,” Phys. Rev. A 85, 043836 (2012). [CrossRef]

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

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

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

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

42. Z. Guan, X.-X. Zhou, and X.-B. Bian, “High-order-harmonic generation from periodic potentials driven by few-cycle laser pulses,” Phys. Rev. A 93, 033852 (2016). [CrossRef]

43. T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto, M. Kuwata-Gonokami, and K. L. Ishikawa, “Trajectory analysis of high-order-harmonic generation from periodic crystals,” Phys. Rev. A 95, 043416 (2017). [CrossRef]

44. L. Li, P. Lan, X. Zhu, T. Huang, Q. Zhang, M. Lein, and P. Lu, “Reciprocal-space-trajectory perspective on high-harmonic generation in solids,” Phys. Rev. Lett. 122, 193901 (2019). [CrossRef]

45. E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,” Phys. Rev. Lett. 52, 997–1000 (1984). [CrossRef]

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

47. N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, “Ultrafast modification of Hubbard U in a strongly correlated material: Ab initio high-harmonic generation in NiO,” Phys. Rev. Lett. 121, 097402 (2018). [CrossRef]

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

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

50. S. V. B. Jensen and L. B. Madsen, “Edge-state and bulklike laser-induced correlation effects in high-harmonic generation from a linear chain,” Phys. Rev. B 104, 054309 (2021). [CrossRef]

51. O. Neufeld, N. Tancogne-Dejean, U. De Giovannini, H. Hübener, and A. Rubio, “Light-driven extremely nonlinear bulk photogalvanic currents,” Phys. Rev. Lett. 127, 126601 (2021). [CrossRef]

52. D. Golde, T. Meier, and S. W. Koch, “High harmonics generated in semiconductor nanostructures by the coupled dynamics of optical inter- and intraband excitations,” Phys. Rev. B 77, 075330 (2008). [CrossRef]

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

54. M. Kira and S. W. Koch, Semiconductor Quantum Optics (Cambridge University, 2012).

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

56. L. Yue and M. B. Gaarde, “Structure gauges and laser gauges for the semiconductor Bloch equations in high-order harmonic generation in solids,” Phys. Rev. A 101, 053411 (2020). [CrossRef]

57. R. Resta, “Macroscopic polarization in crystalline dielectrics: the geometric phase approach,” Rev. Mod. Phys. 66, 899–915 (1994). [CrossRef]

58. D. Vanderbilt, Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators (Cambridge University, 2018).

59. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Maximally localized Wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012). [CrossRef]

60. M. V. Berry, “Quantal phase factors accompanying adiabatic changes,” Proc. R. Soc. Lond. A 392, 45–57 (1984). [CrossRef]

61. R. D. King-Smith and D. Vanderbilt, “Theory of polarization of crystalline solids,” Phys. Rev. B 47, 1651–1654 (1993). [CrossRef]

62. D. J. Thouless, “Quantization of particle transport,” Phys. Rev. B 27, 6083–6087 (1983). [CrossRef]

63. C. Aversa and J. E. Sipe, “Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis,” Phys. Rev. B 52, 14636–14645 (1995). [CrossRef]

64. P. Földi, “Gauge invariance and interpretation of interband and intraband processes in high-order harmonic generation from bulk solids,” Phys. Rev. B 96, 035112 (2017). [CrossRef]

65. K. S. Virk and J. E. Sipe, “Semiconductor optics in length gauge: A general numerical approach,” Phys. Rev. B 76, 035213 (2007). [CrossRef]

66. A. Taghizadeh and T. G. Pedersen, “Gauge invariance of excitonic linear and nonlinear optical response,” Phys. Rev. B 97, 205432 (2018). [CrossRef]

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

68. T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,” Rev. Mod. Phys. 72, 545–591 (2000). [CrossRef]

69. M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional optical pulse propagation equation,” Phys. Rev. Lett. 89, 283902 (2002). [CrossRef]

70. R. W. Boyd, Nonlinear Optics, 3rd ed. (Academic, 2008).

71. M. B. Gaarde, C. Buth, J. L. Tate, and K. J. Schafer, “Transient absorption and reshaping of ultrafast XUV light by laser-dressed helium,” Phys. Rev. A 83, 013419 (2011). [CrossRef]

72. K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe, and G. F. Bertsch, “Time-dependent density functional theory for strong electromagnetic fields in crystalline solids,” Phys. Rev. B 85, 045134 (2012). [CrossRef]

73. J. J. Sakurai and E. D. Commins, Modern Quantum Mechanics, Revised Edition (Addison-Wesley, 1994).

74. D. Xiao, M.-C. Chang, and Q. Niu, “Berry phase effects on electronic properties,” Rev. Mod. Phys. 82, 1959 (2010). [CrossRef]

75. L. B. Madsen, “Different forms of laser–matter interaction operators and expansion in adiabatic states,” Eur. Phys. J. Spec. Top. (2021). [CrossRef]

76. M. Gradhand, D. V. Fedorov, F. Pientka, P. Zahn, I. Mertig, and B. L. Györffy, “First-principle calculations of the Berry curvature of Bloch states for charge and spin transport of electrons,” J. Phys. Condens. Matter 24, 213202 (2012). [CrossRef]

77. J. Zak, “Berry’s phase for energy bands in solids,” Phys. Rev. Lett. 62, 2747–2750 (1989). [CrossRef]

78. G. H. Wannier, “The structure of electronic excitation levels in insulating crystals,” Phys. Rev. 52, 191–197 (1937). [CrossRef]

79. G. H. Wannier, “Dynamics of band electrons in electric and magnetic fields,” Rev. Mod. Phys. 34, 645–655 (1962). [CrossRef]

80. W. Kohn, “Analytic properties of Bloch waves and Wannier functions,” Phys. Rev. 115, 809–821 (1959). [CrossRef]

81. X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, “Ab initio calculation of the anomalous Hall conductivity by Wannier interpolation,” Phys. Rev. B 74, 195118 (2006). [CrossRef]

82. B. H. Bransden and C. J. Joachain, Physics of Atoms and Molecules (Pearson Education, 2003).

83. M. Korbman, S. Y. Kruchinin, and V. S. Yakovlev, “Quantum beats in the polarization response of a dielectric to intense few-cycle laser pulses,” New J. Phys. 15, 013006 (2013). [CrossRef]

84. I. Souza, J. Íñiguez, and D. Vanderbilt, “Dynamics of Berry-phase polarization in time-dependent electric fields,” Phys. Rev. B 69, 085106 (2004). [CrossRef]

85. G. Ernotte, T. J. Hammond, and M. Taucer, “A gauge-invariant formulation of interband and intraband currents in solids,” Phys. Rev. B 98, 235202 (2018). [CrossRef]

86. J. E. Sipe and A. I. Shkrebtii, “Second-order optical response in semiconductors,” Phys. Rev. B 61, 5337–5352 (2000). [CrossRef]

87. G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes, and N. M. R. Peres, “Gauge covariances and nonlinear optical responses,” Phys. Rev. B 96, 035431 (2017). [CrossRef]

88. J. Lu, E. F. Cunningham, Y. S. You, D. A. Reis, and S. Ghimire, “Interferometry of dipole phase in high harmonics from solids,” Nat. Photonics 13, 96–100 (2019). [CrossRef]

89. I. Floss, C. Lemell, K. Yabana, and J. Burgdörfer, “Incorporating decoherence into solid-state time-dependent density functional theory,” Phys. Rev. B 99, 224301 (2019). [CrossRef]

90. I. Kilen, M. Kolesik, J. Hader, J. V. Moloney, U. Huttner, M. K. Hagen, and S. W. Koch, “Propagation induced dephasing in semiconductor high-harmonic generation,” Phys. Rev. Lett. 125, 083901 (2020). [CrossRef]

91. T.-Y. Du, D. Tang, and X.-B. Bian, “Subcycle interference in high-order harmonic generation from solids,” Phys. Rev. A 98, 063416 (2018). [CrossRef]

92. L. H. Thong, C. Ngo, H. T. Duc, X. Song, and T. Meier, “Microscopic analysis of high harmonic generation in semiconductors with degenerate bands,” Phys. Rev. B 103, 085201 (2021). [CrossRef]

93. R. E. F. Silva, F. Martín, and M. Ivanov, “High harmonic generation in crystals using maximally localized Wannier functions,” Phys. Rev. B 100, 195201 (2019). [CrossRef]

94. J. Wilhelm, P. Grössing, A. Seith, J. Crewse, M. Nitsch, L. Weigl, C. Schmid, and F. Evers, “Semiconductor Bloch-equations formalism: Derivation and application to high-harmonic generation from Dirac fermions,” Phys. Rev. B 103, 125419 (2021). [CrossRef]

95. Z. Wang, H. Park, Y. H. Lai, J. Xu, C. I. Blaga, F. Yang, P. Agostini, and L. F. DiMauro, “The roles of photo-carrier doping and driving wavelength in high harmonic generation from a semiconductor,” Nat. Commun. 8, 1686 (2017). [CrossRef]

96. K. Kaneshima, Y. Shinohara, K. Takeuchi, N. Ishii, K. Imasaka, T. Kaji, S. Ashihara, K. L. Ishikawa, and J. Itatani, “Polarization-resolved study of high harmonics from bulk semiconductors,” Phys. Rev. Lett. 120, 243903 (2018). [CrossRef]

97. N. Klemke, N. Tancogne-Dejean, G. M. Rossi, Y. Yang, F. Scheiba, R. E. Mainz, G. Di Sciacca, A. Rubio, F. X. Kärtner, and O. D. Mücke, “Polarization-state-resolved high-harmonic spectroscopy of solids,” Nat. Commun. 10, 1319 (2019). [CrossRef]

98. H. B. Banks, Q. Wu, D. C. Valovcin, S. Mack, A. C. Gossard, L. Pfeiffer, R.-B. Liu, and M. S. Sherwin, “Dynamical birefringence: Electron-hole recollisions as probes of Berry curvature,” Phys. Rev. X 7, 041042 (2017). [CrossRef]

99. R. E. F. Silva, Á. Jiménez-Galán, B. Amorim, O. Smirnova, and M. Ivanov, “Topological strong-field physics on sub-laser-cycle timescale,” Nat. Photonics 13, 849–854 (2019). [CrossRef]

100. J. B. Krieger and G. J. Iafrate, “Time evolution of Bloch electrons in a homogeneous electric field,” Phys. Rev. B 33, 5494–5500 (1986). [CrossRef]

101. J. Li, X. Zhang, S. Fu, Y. Feng, B. Hu, and H. Du, “Phase invariance of the semiconductor Bloch equations,” Phys. Rev. A 100, 043404 (2019). [CrossRef]

102. A. Taghizadeh, F. Hipolito, and T. G. Pedersen, “Linear and nonlinear optical response of crystals using length and velocity gauges: Effect of basis truncation,” Phys. Rev. B 96, 195413 (2017). [CrossRef]

103. J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, 1999).

104. C. Chandre, S. Wiggins, and T. Uzer, “Time-frequency analysis of chaotic systems,” Physica D 181, 171–196 (2003). [CrossRef]

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

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

107. L. Yue and M. B. Gaarde, “Expanded view of electron-hole recollisions in solid-state high-order harmonic generation: Full-Brillouin-zone tunneling and imperfect recollisions,” Phys. Rev. A 103, 063105 (2021). [CrossRef]

108. A. M. Parks, G. Ernotte, A. Thorpe, C. R. McDonald, P. B. Corkum, M. Taucer, and T. Brabec, “Wannier quasi-classical approach to high harmonic generation in semiconductors,” Optica 7, 1764–1772 (2020). [CrossRef]

109. F. Navarrete, M. F. Ciappina, and U. Thumm, “Crystal-momentum-resolved contributions to high-order harmonic generation in solids,” Phys. Rev. A 100, 033405 (2019). [CrossRef]

110. G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists (International Edition), 6th ed. (Academic, 2005).

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

112. L. Yue and M. B. Gaarde, “Imperfect recollisions in high-harmonic generation in solids,” Phys. Rev. Lett. 124, 153204 (2020). [CrossRef]

113. J. A. Crosse, X. Xu, M. S. Sherwin, and R. B. Liu, “Theory of low-power ultra-broadband terahertz sideband generation in bi-layer graphene,” Nat. Commun. 5, 4854 (2014). [CrossRef]

114. X. Zhang, J. Li, Z. Zhou, S. Yue, H. Du, L. Fu, and H.-G. Luo, “Ellipticity dependence transition induced by dynamical Bloch oscillations,” Phys. Rev. B 99, 014304 (2019). [CrossRef]

115. S. Y. Kruchinin, F. Krausz, and V. S. Yakovlev, “Colloquium: Strong-field phenomena in periodic systems,” Rev. Mod. Phys. 90, 021002 (2018). [CrossRef]

116. C. Yu, U. Saalmann, and J. M. Rost, “High harmonics from backscattering of delocalized electrons,” arXiv:2102.11208 (2021).

117. J. Chen, Q. Xia, and L. Fu, “Spectral caustics of high-order harmonics in one-dimensional periodic crystals,” Opt. Lett. 46, 2248–2251 (2021). [CrossRef]

118. L. Li, P. Lan, X. Zhu, and P. Lu, “Huygens-Fresnel picture for high harmonic generation in solids,” Phys. Rev. Lett. 127, 223201 (2021). [CrossRef]

119. S. Yamada and K. Yabana, “Determining the optimum thickness for high harmonic generation from nanoscale thin films: An ab initio computational study,” Phys. Rev. B 103, 155426 (2021). [CrossRef]

120. T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. 78, 3282–3285 (1997). [CrossRef]

121. M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70, 036604 (2004). [CrossRef]

122. P. W. Milonni and J. H. Eberly, Laser Physics (Wiley, 2010).

123. J. Peatross and M. Ware, “Physics of light and optics: A free online textbook,” in Frontiers in Optics 2010/Laser Science XXVI (Optical Society of America, 2010), paper JWA64.

124. C. Q. Abadie, M. Wu, and M. B. Gaarde, “Spatiotemporal filtering of high harmonics in solids,” Opt. Lett. 43, 5339–5342 (2018). [CrossRef]

125. C. Jürß and D. Bauer, “Helicity flip of high-order harmonic photons in Haldane nanoribbons,” Phys. Rev. A 102, 043105 (2020). [CrossRef]

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

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

128. R. E. F. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M. Ivanov, “High-harmonic spectroscopy of ultrafast many-body dynamics in strongly correlated systems,” Nat. Photonics 12, 266–270 (2018). [CrossRef]

129. Y. Murakami, S. Takayoshi, A. Koga, and P. Werner, “High-harmonic generation in one-dimensional Mott insulators,” Phys. Rev. B 103, 035110 (2021). [CrossRef]

130. C. Orthodoxou, A. Zaïr, and G. H. Booth, “High harmonic generation in two-dimensional Mott insulators,” npj Quantum Mater. 6, 76 (2021). [CrossRef]

LG Bloch or Houston Scheme	VG Bloch Scheme
Small number of bands for convergence	No periodic structure gauge requirement
Easier inclusion of dephasing	Less number of $k$ points for convergence
Easier current decomposition	Easier treatment of degenerate bands

Introduction to theory of high-harmonic generation in solids: tutorial

Abstract

1. INTRODUCTION

2. ADIABATIC STATES, BERRY CONNECTIONS, AND CURVATURES

A. Adiabatic States

B. Berry Phase, Connection, and Curvature

3. TIME-INDEPENDENT PROBLEM AND STRUCTURE GAUGES

A. Bloch States and Coupling Matrix Elements

B. Random Structure Gauge

C. Parallel Transport Structure Gauge

D. Periodic Structure Gauge

E. Wannier Gauge

F. Structure Gauges for Degenerate Bands

4. TIME-DEPENDENT PROBLEM AND LASER GAUGES

A. Velocity Gauge EOM in the Bloch Basis

B. Length Gauge EOM in Bloch Basis

C. EOM in the Adiabatic Houston Basis

D. Gauge Comparisons

E. HHG Calculation and Time Profiles

5. SADDLE-POINT EQUATIONS AND THE RECOLLISION MODEL FOR HHG

A. Saddle-Point Method for HHG

B. Semiclassical Recollision Model

6. MACROSCOPIC PROPAGATION

A. Spatio-Spectral Properties of Solid-State HHG

7. SUMMARY AND OUTLOOK

Funding

Acknowledgment

Disclosures

Data Availability

REFERENCES

Data Availability

Cited By

Figures (5)

Tables (1)

Equations (118)

Journal of the Optical Society of America B