High-harmonic generation in solids

TL;DR

This paper combines analytical and numerical methods to study high-harmonic generation in solids, revealing the importance of entire Brillouin zone contributions and providing predictions that align with experimental data.

Contribution

It introduces a model using Houston states and a Kronig-Penney potential to analyze HHG, including multiband effects and BZ contributions, advancing understanding of solid-state HHG mechanisms.

Findings

High-harmonic spectra depend on entire first Brillouin zone contributions.

Derived cutoff harmonic orders match analytical predictions and experimental data.

Multiband calculations provide detailed insights into intra- and interband HHG contributions.

Abstract

We analytically and numerically investigate the emission of high-harmonic radiation from model solids by intense few-cycle mid-infrared laser pulses. In single-active-electron approximation, we expand the active electron's wavefunction in a basis of adiabatic Houston states and describe the solid's electronic band structure in terms of an adjustable Kronig-Penney model potential. For high-harmonic generation (HHG) from MgO crystals, we examine spectra from two-band and converged multiband numerical calculations. We discuss the characteristics of intra- and interband contributions to the HHG spectrum for computations including initial crystal momenta either from the $\Gamma-$point at the center of the first Brioullin zone (BZ) only or from the entire first BZ, demonstrating relevant contributions from the entire first BZ. From the numerically calculated spectra we derive cutoff harmonic orders as a function of the laser peak intensity that compare favorably with our analytical stationary-phase-approximation predictions and published data.

Equations

$$\left[\frac{1}{2}{\left(\hat{p}+A(t)\right)}^2+V(x)\right]|\psi (t)⟩=i\frac{\partial }{\partial t}|\psi (t)⟩,$$

$$A(t)$$

$$\displaystyle=\frac{A_{0}}{2T}\sin(\omega_{0}t)\left\{\begin{tabular}[]{cl}$t$&\;, $0\leq t\leq 2T\;\,$\\ $2T$&\;, $2T\leq t\leq 8T\;\,$\\ $(10T-t)$&\;, $8T\leq t\leq 10T\;\,$\end{tabular}\right..$$

$$\cos (ak)=\cos (a\sqrt{2{\epsilon }_{nk}})+\frac{V_0}{\sqrt{2{\epsilon }_{nk}}}\sin (a\sqrt{2{\epsilon }_{nk}})$$

$$|{\psi }_k(t)⟩=e^{-iA(t)x}\sum _n{B}_{nk}(t)e^{\left[-i{\int }_0^t{\epsilon }_{n\kappa (t^{\prime })}d{t}^{\prime }\right]}|{\varphi }_{n\kappa (t)}⟩,$$

$$i{\dot{B}}_{nk}(t)$$

$$\times \exp \left[i{\int }_0^t\Delta {\epsilon }_{\kappa (t)}^{n{n}^{\prime }}d{t}^{\prime }\right].$$

$$D_{\kappa (t)}^{n{n}^{\prime }}$$

$$D_{\kappa (t)}^{nn}$$

$$P_{\kappa (t)}^{n{n}^{\prime }}=\frac{1}{a}{\int }_0^a{\varphi }_{n\kappa (t)}^{\ast }(x)\frac{1}{i}\frac{\partial }{\partial x}{\varphi }_{n^{\prime }\kappa (t)}(x)dx.$$

$$P_{\kappa (t)}^{nn}=\frac{\partial {\epsilon }_{n\kappa (t)}}{\partial \kappa (t)},$$

$$\kappa (t)=k+A(t),$$

$$\left[\frac{{\hat{p}}^2}{2}+V(x)\right]|{\varphi }_{n\kappa (t)}⟩={\epsilon }_{n\kappa (t)}|{\varphi }_{n\kappa (t)}⟩.$$

$$\hat{\rho }(t)=|{\psi }_k(t)⟩⟨{\psi }_k(t)|,$$

$${\rho }_k^{n{n}^{\prime }}(t)$$

$$=⟨{\varphi }_{n\kappa (t)}|{\hat{\rho }}_k|{\varphi }_{n^{\prime }\kappa (t)}⟩e^{i{\int }_0^t\Delta {\epsilon }_{\kappa (t^{\prime })}^{n{n}^{\prime }}}d{t}^{\prime }$$

$$i{\dot{\rho }}_k^{n{n}^{\prime }}(t)$$

$$-D_{\kappa (t)}^{n^{\prime }{n}^{\prime \prime }}{e}^{i{\int }_0^t\Delta {\epsilon }_{\kappa (t^{\prime })}^{n^{\prime \prime }{n}^{\prime }}d{t}^{\prime }}{\rho }_k^{n{n}^{\prime \prime }}(t)].$$

$$J_k(t)=J_k^{ra}(t)+J_k^{er}(t),$$

$$J_k^{ra}(t)$$

$$J_k^{er}(t)$$

$$Y_k(\omega )$$

$$=Y_k^{ra}(\omega )+Y_k^{er}(\omega )+2{\hat{J}}_k^{er}(\omega ){\hat{J}}_k^{ra}(\omega ),$$

$$Y\left(\omega \right)={\left|{\int }_{-\infty }^{\infty }dt{e}^{-i\omega t}{\int }_{BZ}dk{J}_k(t)\right|}^2,$$

$$J_k^{ra}(t)$$

$$J_k^{er}(t)$$

$$\Delta J_k^{er}(t)$$

$${\bar{J}}_k^{er}(t)$$

$$S(k,t)={\int }_0^t\Delta {\epsilon }_{\kappa (t^{\prime })}^{cv}d{t}^{\prime }.$$

$$i{\dot{\rho }}_k^{cv}(t)=e^{iS(k,t)}\Delta {\rho }_k^{vc}(t)E(t)D_{\kappa (t)}^{vc},$$

$$\Delta {\rho }_k^{vc}(t)=\left[{\rho }_k^{vv}(t)-{\rho }_k^{cc}(t)\right],$$

$$i{\dot{\rho }}_k^{cc}(t)=E(t)D_{\kappa (t)}^{vc}{e}^{-iS(k,t)}{\rho }^{cv}(t)-\text{c.c.}.$$

$$J_k^{ra}(t)=-P_{\kappa (t)}^{vv}+{\rho }_k^{cc}(t)\Delta P_{\kappa (t)}^{vc},$$

Full text

††preprint: APS/123-QED

I Introduction

Exposed to intense laser fields gases and solids emit a spectrum of radiation that is strongly enhanced at and near frequencies corresponding to multiples of the driving laser frequency. Over the past two decades, this high-order harmonic generation (HHG) process has been carefully investigated in atomic gases, and the underlying generation mechanism - the emission of radiation by laser-electric-field-driven rescattered electrons - is well understood Le et al. (2009); Krausz and Ivanov (2009). While solids are being discussed theoretically for decades in view of their large electronic density possibly enabling the design of high-intensity sources of harmonic radiation Plaja and Roso-Franco (1992); Protopapas et al. (1997); Faisal and Kamiński (1996); Faisal and Kaminski (1997), HHG from solids has remained a matter of debate Ghimire et al. (2011); Wu et al. (2015); Vampa et al. (2014). Experimentally, it was first carefully scrutinized less than a decade ago by Ghimire et al.  Ghimire et al. (2011). Understanding the mechanisms of HHG in solids is an area of emerging research interest and part of the ongoing diversification of attosecond science from the study of atoms and molecules to more complex nanoparticles Li et al. (2016, 2017, 2018a) and solids Leone et al. (2014); Thumm et al. (2015); Ciappina et al. (2017); Kruchinin et al. (2018). This extension of attosecond science holds promise for promoting the development of novel table-top intense high-frequency radiation sources and our understanding of the light-induced electron dynamics in solids, a prerequisite for improved ultrafast electron-optical switches Luu and Wörner (2016); Ghimire and Reis (2019).

Compared to HHG in atomic gases, theoretical investigations of solid HHG have indicated striking new effects, such as multiple plateaus Ndabashimiye et al. (2016); Wu et al. (2015) in HHG spectra and a linear dependence of the HHG cutoff frequency on the peak electric-field strengths of the driving laser Wu et al. (2015); Vampa et al. (2015); Osika et al. (2017); Li et al. (2018b). These characteristics have been revealed by numerically solving either the time-dependent Schrödinger (TDSE) in single-active-electron (SAE) approximation Mücke (2011); Hawkins et al. (2015); Wu et al. (2015); Li et al. (2018b) or semiconductor Bloch equations (SBEs) Golde et al. (2008); Schubert et al. (2014); Meier et al. (1994); Vampa et al. (2014); Jiang et al. (2017, 2018); Li et al. (2018b). SAE-TDSE-based numerical models have employed basis-set expansions of the active electron’s wavefunction using either static Korbman et al. (2013) or adiabatic Krieger and Iafrate (1985); Wu et al. (2015) Bloch states.

SAE solutions of the TDSE can be expressed in terms of density matrices for convenient comparison with the SBE approach that introduces a phenomenological dephasing time to account for relaxation processes Haug and Koch (1998). An advantage of working within the SAE-TDSE framework is that the computational time for solving a system with $n$ electronic bands scales linearly with $n$, while it scales as ${n(n-1)/2}$ in SBE calculations Luu and Wörner (2016). For ${n>20}$ this leads to approximately one order of magnitude difference in computation time (cf., Ref. Wu et al. (2015), where 51 bands are included for solving the TDSE within a static Bloch basis). This is of relevance at high intensities of the driving laser, where calculations with a large number of bands are required to reveal the multiplateau structure of converged HHG spectra  Wu et al. (2015). SAE models have successfully explained the main features of HHG in atomic gases Le et al. (2009); Krausz and Ivanov (2009), which has motivated their transfer to describing HHG in solids Wu et al. (2015); Ndabashimiye et al. (2016).

In this work we develop a numerical model for solving the TDSE in SAE approximation employing a basis-set expansion of the electronic wavefunction in so-called “Houston states” that vary adiabatically with the instantaneous driving-laser electromagnetic field Wu et al. (2015); Lindefelt et al. (2004); Krieger and Iafrate (1985). The use of an adiabatic basis is advantageous for gaining physical insight into the underlying basic mechanisms for HHG in solids, since it allows the identification of two distinct processes, intra- and an interband emission, that operate in different spectral regions. In calculating HHG spectra from solids, we pay attention to the crystal momentum of the initial state and scrutinize contributions from different crystal momenta in the entire first Brillouin zone (BZ). In particular, we find that HHG spectra calculated by only including the initial crystal momentum at the center of the BZ (the $\Gamma-$point) Wu et al. (2015) noticeably differ from calculations that include initial crystal momenta in the entire BZ Vampa et al. (2014). It has been recently proposed that the range of crystal momenta used in SAE calculations should be used as an adjustable parameter for getting converged numerical HHG spectra Li et al. (2018b). In the present work we revisit this suggestion and analytically determine, based on a stationary-phase approximation, the range of initial crystal momenta needed for the computation of HHG spectra at a given accuracy.

We organized this paper as follows. In Sec. II we describe our theoretical framework. In particular, in Sec. II.1 we solve the TDSE by expanding the active electron’s wavefunction in an adiabatic basis (Sec. II.1.1), compare our approach with a density-matrix formulation of solid HHG (Sec. II.1.2), and show how the observables of interest in this work, intra- and interband yields, are retrieved from our numerical results (Sec. II.1.3). In Sec. II.2 we discuss solid HHG for a simplified two-band system distinuishing intra- (Sec. II.2.1) and interband (Sec. II.2.2) yields. In Sec. II.3 we analyze interband HHG within a stationary-phase approximation (Sec. II.3.1), which allows us to estimate the relevant range of crystal momenta $k$ (Sec. II.3.2) that need to be included when adding HHG contributions from different $k$ in the first BZ (Sec. II.3.3). In Sec. III we present and discuss our numerical results for HHG in a model MgO crystal. First, based on simplified two-band calculations, we analyze in Sec. III.1 $k$-resolved spectra for a specific field strength of the driving-laser pulse (Sec. III.1.1), field-strength-dependent spectra (Sec. III.1.2), and cut-off harmonic orders as a function of the field strengths (Sec. III.1.3). Next, in Sec. III.2, we analyze field-strength-dependent HHG spectra (Sec. III.2.1) and cutoff harmonic orders (Sec. III.2.2) for calculations that are converged in the number of included electronic bands. In several appendices, we add details of our theoretical analysis. We use atomic units (${q_{e}=m_{e}=\hbar=1}$) throughout this work, unless specified otherwise.

II Theory

II.1 Single-active-electron solution of the TDSE

We solve the TDSE,

[TABLE]

subject to the interaction of the active electron with both, a one-dimensional infinitely extended solid and an infrared (IR) external laser field ${E(t)}$. We represent the solid by a periodic potential ${V(x)=V(x+a)}$, with lattice constant $a$, and the laser field by a 10-cycle “flat-top” vector potential:

[TABLE]

${\hat{p}=-i\frac{\partial}{\partial x}}$ denotes the momentum operator and $A_{0}$, $\omega_{0}$, and ${T=2\pi/\omega_{0}}$ the external vector-potential amplitude, frequency, and period, respectively. Since in our numerical simulation (Sec. III) the driving-laser wavelength is three orders of magnitude larger than $a$ and the classical excursion range of the active electron in the laser field, in solving Eq.(1) we can safely invoke the dipole approximation, ${A(x,\,t)\approx A(t)}$.

We model ${V(x)}$ as a Kronig-Penney potential Kronig et al. (1931); Atkins and Friedman (2005), which yields the dispersion relation

[TABLE]

for the valence (${n=v}$) and conduction bands (${n=c}$). The potential strength $V_{0}$ is adjusted to match the electronic band structure of the solid. While Eq. (3) needs to be solved numerically, the Kronig-Penney model potential provides a convenient basis set of 2-fold-degenerate orthonormal eigenstates and allows us to calculate transition matrix elements in a closed analytical form (see Sec. II.1.1 below and Appendix A).

II.1.1 Expansion in Houston states

Expanding solutions of Eq. (1),

[TABLE]

in terms of Houston states ${\left|\phi_{n\kappa(t)}\right>}$, results in the set of coupled equations

[TABLE]

In Eq. (5) we define the energy difference between Houston states Krieger and Iafrate (1985); Wu et al. (2015) (also referred to as “band-gap energy”) ${\Delta\varepsilon_{\kappa(t)}^{nn^{\prime}}=\varepsilon_{n\kappa(t)}-\varepsilon_{n^{\prime}\kappa(t)}}$ and the transition dipole moments (TDMs)

[TABLE]

The TDMs are given in terms of the momentum-operator matrix elements Lindefelt et al. (2004)

[TABLE]

The diagonal elements,

[TABLE]

are related to the band energy $\varepsilon_{n{\kappa(t)}}$ and correspond to the group velocity of an electron wave packet in band $n$. Since the elements ${P^{nn^{\prime}}_{\kappa(t)}}$ are real, the TDMs satisfy ${D^{nn^{\prime}}_{\kappa(t)}=-D^{n^{\prime}n}_{\kappa(t)}}$ and ${\left(D^{nn^{\prime}}_{\kappa(t)}\right)^{*}=D^{n^{\prime}n}_{\kappa(t)}}$ (Appendix A).

Houston states are adiabatic in the field-dressed time-dependent crystal momentum

[TABLE]

and solve the Schrödinger equation

[TABLE]

They can be viewed as adiabatic Bloch states, with $\kappa(t)$ replacing the Bloch-state crystal momentum $k$. ${\kappa(t)}$ thus parameterizes the field-driven electronic evolution under the influence of the IR pulse out of an initial state with crystal momentum $k$ (Fig. 1). As for ordinary Bloch functions Ashcroft and Mermin (1976), Houston states with different initial (field-free) crystal momentum $k$ or different band indices $n$ are not coupled by the Hamiltonian in Eq. (10), yet evolve differently. Since Houston states for different initial momenta $k$ explore the electronic band $n$ distinctively, they will be referred to as “$k-$channels” in this work.

II.1.2 Density-matrix formulation

The electronic evolution described in Eq. (5) can also be expressed in terms of the density operator,

[TABLE]

in the pure state ${\left|\psi_{k}(t)\right>}$ given by Eq. (4). The matrix elements

[TABLE]

represent band populations for ${n=n^{\prime}}$. Their time evolution is obtained by solving

[TABLE]

II.1.3 Intra- and interband yield

The electronic current in each ${k-}$channel,

[TABLE]

consists of intra- and interband contributions,

[TABLE]

and

[TABLE]

respectively. It defines the spectral HHG yield from a given $k-$channel,

[TABLE]

which, in addition to the intra- [$Y_{k}^{ra}(\omega)$] and interband [$Y_{k}^{er}(\omega)$] yields, includes the interference term $2\hat{J}_{k}^{er}(\omega)\hat{J}_{k}^{ra}(\omega)$.

For the dielectric solid analyzed in this work, the Fermi energy lies in the band gap between the highest occupied band (the valence-band) and the lowest unoccupied band (the first conduction band). Since the valence-band is fully occupied, we obtain the total HHG yield as

[TABLE]

and corresponding expressions for the total intra- [${Y^{ra}(\omega)}$] and interband HHG yields [${Y^{er}(\omega)}$], after including current contributions from all ${k-}$channels, i.e. , from all crystal momenta $k$ in the first BZ.

II.2 HHG mechanism for a two-band system

In this subsection, we restrict the theory developed in Sec. II.1 to two electronic bands: the valence and first conduction band, even though more than two bands are required and included in our converged numerical calculations in Sec. III below. Designating the valence and conduction bands with superscripts $v$ and $c$, respectively, the intra- and interband currents in Eqs. (12) and (13) simplify to

[TABLE]

As detailed in Appendix C, Eq. (17) can be written as

[TABLE]

where we define

[TABLE]

and the action

[TABLE]

Applying Eq. (11) to the two-band system, we obtain

[TABLE]

with the population difference

[TABLE]

and

[TABLE]

Both, the SBE model for HHG in Refs. Vampa et al. (2014); Jiang et al. (2017); Li et al. (2018b) and the expansion of TSDE solutions in a Houston basis, use an adiabatic basis. Note that the SBE model in Refs. Vampa et al. (2014); Jiang et al. (2017); Li et al. (2018b) adds the damping term ${\rho^{cv}_{k}(t)/T_{2}}$ to Eq. (21), with an adjustable damping time $T_{2}$ for the interband coupling Haug and Koch (1998). In addition, it does not include the term ${\Delta J_{k}^{er}(t)}$ in Eq. (18). Therefore, the SBE model Vampa et al. (2014) and our SAE-TSDE cannot be expected to yield identical HHG spectra, not even for ${T_{2}\rightarrow\infty}$.

II.2.1 Intraband yield

For a two-band system, norm preservation demands ${\rho^{vv}_{k}(t)=1-\rho^{cc}_{k}(t)}$, such that the intraband current in Eq. (12) simplifies to

[TABLE]

resulting in the interband yield

[TABLE]

where

[TABLE]

At low laser intensities, the conduction-band population is very small (${\rho^{cc}_{k}(t)\ll 1}$), and, for every ${k-}$channel, the by far dominant contribution to the intraband HHG spectrum is generated by the term ${P^{vv}_{\kappa(t)}}$ in the intraband current. With increasing driving-field intensities, the second term (${\sim\rho^{cc}_{k}}$) gains significance.

Upon integration over the first BZ, the first term in Eq. (24) vanishes, due to symmetry, since ${P^{vv}_{k}}$ is an odd function of $k$ (Fig. 1). The integrated current is therefore given by

[TABLE]

leading to the intraband HHG yield for the two-band system

[TABLE]

At low intensities, even though each ${k-}$channel defines a significant intraband current, due to the symmetry-related cancellation of contributions from crystal momenta ${+k}$ and ${-k}$, the integrated intraband HHG yield is comparatively small.

II.2.2 Interband yield

Calculation of the HHG yield for a specific $k$-channel, according to Eq. (14), requires the Fourier-transformed current

[TABLE]

Replacing ${\rho^{cv}_{k}(t)}$, obtained by integrating Eq. (21), in Eq. (18), the interband current and HHG yield for the two-band system become

[TABLE]

and

[TABLE]

respectively.

II.3 Approximate evaluation of the interband yield

While the numerical HHG spectra discussed in Sec. III below are calculated based on the theory outlined in Secs. II.1 and II.2, we apply in this subsection additional approximations to the interband current in order to derive analytical expressions that reveal additional physical properties of the interband HHG process in solids. We restrict this analysis to vector-potential amplitudes ${A_{0}<\pi/a}$, for which the range of ${\kappa(t)}$ is limited by the width, ${2\pi/a}$, of one BZ (Fig. 1).

II.3.1 Stationary-phase approximation

Setting the TDM in Eq. (30) equal to its value at the band center (${k=0}$) Vampa et al. (2014, 2015), ${D^{vc}_{\kappa(t)}\approx D^{vc}_{0}}$, remembering that according to Eq. (6) ${\left(D^{vc}_{0}\right)^{2}=-\left|D^{vc}_{0}\right|^{2}}$, and applying the frozen valence-band approximation Keldysh (1965) ${|\Delta\rho_{k}^{vc}(t)|\approx 1}$ Keldysh (1965), Eq. (30) simplifies to

[TABLE]

where

[TABLE]

Since ${e^{i\sigma_{\omega}(k,\,t,\,t^{\prime})}}$ rapidly oscillates as a function of $t$ and ${t^{\prime}}$, the integrals in Eq. (32) are dominated by contributions at times ${t=t_{e}}$ and ${t^{\prime}=t_{s}}$ when the phase ${\sigma_{\omega}(k,\,t,\,t^{\prime})}$ is stationary, i.e.,

[TABLE]

At these times we find

[TABLE]

Evaluating Eq. (32) in stationary-phase approximation Lewenstein et al. (1994); Vampa et al. (2014); Wong (2001) now results in

[TABLE]

with the Hessian matrix

[TABLE]

Equation (34) implies that for each ${k-}$channel the cutoff frequency for interband HHG becomes

[TABLE]

Equation (35) cannot be fulfilled in our Kronig-Penney model for real-valued times and energies because the bandgap is nonzero across the entire BZ. We designate the maximal cutoff energy as ${\Delta\varepsilon_{max}^{cv}}$. Complying with Eq. (35) requests allowing for complex-valued times, energies, and crystal momenta $k$. Designating the complex crystal momentum as $K$, we analytically continue Eq. (35) and the transcendental Eq. (3) into the complex $K-$plane. Even though complex roots ${K_{s}=K(t_{s})}$ of Eq. (35) can only be obtained numerically, we get further insight into the interband HHG process by Taylor-expanding about ${K=0}$,

[TABLE]

While the expansion coefficients are in general complex, we show in Appendix B that ${\frac{\text{d}\Delta\varepsilon_{0}^{cv}}{\text{d}K}=0}$ and ${\Im\left[\frac{\text{d}^{2}\Delta\varepsilon_{0}^{cv}}{\text{d}K^{2}}\right]=0}$, such that

[TABLE]

where $\Im$ stands for “imaginary part of”. We evaluate the real-valued reduced effective mass,

[TABLE]

at the band center,

[TABLE]

in terms of the valence- and conduction-band effective masses

[TABLE]

Note that the derivatives in Eq. (39) are taken along the real axis and ${\varepsilon_{nk}=\Re[\varepsilon_{nK}]}$, where $\Re$ stands for “real part of”. Since the first term and the coefficient of the quadratic term in Eq. (38) are real, the roots $K_{s}$ are purely imaginary and correspond to interband transitions at the $\Gamma-$point (${k=0}$).

In Appendix D we show that, for a continuous-wave driving field of the form ${A_{0}\sin(\omega_{0}t)}$ and for ${\omega_{0}\sqrt{2\Delta\varepsilon_{0}^{cv}m_{0}^{*}}<E_{0}}$, the interband HHG yield is given by

[TABLE]

Applying the effective mass theorem Ashcroft and Mermin (1976),

[TABLE]

we obtain

[TABLE]

This allows us to rewrite Eq. (40) as

[TABLE]

As expected, the interband HHG yield increases with decreasing band gap ${\Delta\varepsilon_{0}^{cv}}$ and increasing TDM ${\left|D^{vc}_{0}\right|}$.

II.3.2 Relevant k-range for the interband HHG

The repeated-zone scheme in Fig. 1 shows that, for ${k>A_{0}}$ or ${k<-A_{0}}$, $0<|\kappa(t)|<2\pi/a$. Hence, $|\kappa(t)|$ does not reach the limits [math] and $2\pi/a$, Eq. (35) cannot be satisfied, and contributions to HHG are restricted to the interval ${k\in(-A_{0},\,A_{0})}$ within the first BZ. For the symmetrical dispersion relation of the Kronig-Penney model (Fig. 1), the interband HHG yield in Eq. (42) is largest at the $\Gamma-$point.

With the upper limit for the integrated interband yield

[TABLE]

and since ${\left[1-(k/A_{0})^{2}\right]^{1/2}\leq 1}$ in Eq. (42), we find the contribution of each $k$-channel to the interband yield to be limited by

[TABLE]

We can now estimate the range ${k_{max}(A_{0})}$ of initial-state crystal momenta ${k\in\left(-\frac{k_{max}(A_{0})}{2},\frac{k_{max}(A_{0})}{2}\right)}$ that yield relative contributions larger than ${10^{-N}}$ times the maximal yield $Y^{er}_{0}(\omega)$,

[TABLE]

as

[TABLE]

II.3.3 Integrated interband yield

In this subsection we examine the net contribution to the interband current from all $k$-channels in the first BZ to obtain the observable integrated interband yield ${Y^{er}(\omega)}$. For numerical applications, the numerical effort can possibly be reduced by limiting the integration range to ${k_{max}(A_{0})}$, as determined in the previous subsection. Whether this is possible depends on the specific laser parameters and solid electronic structure.

Equation (29) leads to the interband current, integrated over the whole BZ,

[TABLE]

For constant TDMs and a frozen valence-band population, we arrive at the saddle-point conditions

[TABLE]

with ${\sigma_{\omega}(k,\,t^{\prime},\,t)}$ defined in Eq. (33). These conditions imply

[TABLE]

where ${\kappa_{r}(t)=k_{r}+A(t)}$, and allow us to determine the roots $k_{r}$, $t_{s}$, and $t_{e}$ numerically.

Equations (47) and (49) are equivalent to Eqs. (35) and (34). Condition (48) arises due to the integration over the BZ. It expresses the requirement of the excited photoelectron wave packet, moving with group velocity ${P^{cc}_{k_{r}+A(t)}}$, and the residual hole wave packet, propagating with group velocity $P^{vv}_{k_{r}+A(t)}$, to recombine at time $t_{e}$ after their birth at time $t_{s}$, while emitting a photon with energy ${\Delta\varepsilon_{\kappa_{r}(t_{e})}^{cv}}$. We show in Appendix E that, since ${\Im[t_{s}]\propto\omega_{0}\sqrt{2\Delta\varepsilon^{cv}_{0}m_{0}^{*}}}$ and the matrix elements $P^{nn}_{k}$ are odd functions of $k$, Eq. (48) can be approximated as

[TABLE]

Even though within the Kronig-Penney model the band-gap energy ${\Delta\varepsilon_{k}^{cv}}$ grows continuously from the center to the edge of the first BZ, the $k-$channel at ${k=k_{r}^{c}}$ with the largest frequency,

[TABLE]

lies in the range ${k_{r}\in(-A_{0},\,A_{0})}$, as seen in Sec. II.3.2, and needs to be determined numerically (see Sec. III.2.1 for a specific numerical example).

In analogy to Eq. (36), we obtain the Fourier-transformed net interband current as

[TABLE]

This expression is equivalent to

[TABLE]

where

[TABLE]

resulting in the approximated integrated interband yield

[TABLE]

III Numerical Results

For our numerical applications of the theoretical model described in Sec. II, we adopt the laser wavelength (3250 nm) and pulse duration (10 optical cycles) of the experiment by Ghimire et al.  Ghimire et al. (2011) and the temporal pulse profile given by Eq. (2). We model the electronic structure of MgO based on the Kronig-Penney model potential [Eq. (56)] with an interlayer spacing of ${a=8}$ a.u. Xu and Ching (1991) and adjust the potential strength to ${V_{0}=22.345}$ eV in order to reproduce the bandgap energy between the valence and conduction band at the $\Gamma-$point, 4.19 eV, obtained by Xu and Ching Xu and Ching (1991) from a full-dimensionality orthogonalized-linear-combination-of-atomic-orbitals (OLCAO) calculation within the framework of density-functional theory (DFT) in local density approximation (LDA). These values of the Kronig-Penney potential parameters result in a local bandgap at the BZ edge (at ${k=\pm\pi/a}$) of 13.6 eV, in very good agreement with the local bandgap at the $X$ point of 13 eV computed by Xu and Ching (Fig. 2).

In this section, we perform a systematic numerical study of the contributions to the HHG spectrum from different $k$ channels within the entire first BZ. We first discuss results obtained by restricting the external-field-driven electron dynamics to the valence and conduction band in Sec. III.1, before presenting converged HHG spectra obtained by including up to 13 electronic bands of MgO in Sec. III.2. For all calculations we employed a fourth-order Runge-Kutta algorithm to numerically solve Eq. (5) at 400 equally spaced ${k-}$points in the first BZ.

III.1 HHG spectra in two-band approximation

In order to understand the basic mechanisms of intra- and interband HHG in solids, we complement our theoretical analysis of HHG by a two-band system in Sec. II.2 above, with the numerical solution of Eq. (5), restricted to the lowest (highest) conduction (valence) band of MgO. While this approach can only provide acceptable results at moderate intensities of the driving field (${A_{0}<\pi/a}$), it fails at higher field strengths, for which the inclusion of more bands is mandatory (Sec. III.2).

III.1.1 Lattice-momentum-resolved contributions to HHG

To reveal the characteristics of HHG spectra, including all $k$-channels in the first BZ, we performed a calculation at a field strength of ${0.13\,\text{V/\AA}}$. Figures 3 (a) - 3 (c) display two-band HHG spectra as a function of the lattice momentum $k$. While, according to Eq. (14), the total yield is not equal to the sum of the intra- and interband yield, it is instructive to examine intra- and interband spectra separately. The comparison of Figs. 3(a) and 3(b) shows that intraband harmonics dominate HHG at harmonic photon energies below the band gap, ${\Delta\varepsilon_{0}^{cv}\approx 11\omega_{0}}$, while above this threshold mainly interband harmonics contribute to the total HHG yield in Fig. 3(c). The observation that the band gap establishes a threshold between intra- and interband HHG is consistent with energy conservation, requiring electronic probability density to cross the local band gap, ${\Delta\varepsilon_{\kappa(t)}^{cv}}$, before contributing to the interband current and thus to the interband HHG.

In view of Eq. (24) and the semiclassical interpretation of $P^{vv}_{\kappa(t)}$ as the photoelectron group velocity in the valence band, at low electric-field strengths we expect the intraband yield [Eq. (25)] to be dominated by the Fourier transform of $P^{vv}_{\kappa(t)}$ (Fig. 4). Indeed, below the interband-gap threshold (11th harmonic), the yields in Figs. 3(a) and 4 very closely resemble each other and are identical with respect to their zeros along the $k$-axis for every given harmonic order: both spectra exhibit even and odd harmonics that vary over the first BZ. At the $\Gamma-$point and edge of the first BZ both spectra only include odd harmonics.

Below the interband-gap threshold, the interband spectrum in Fig. 3(b) bears some similarity with the intraband yield in Fig. 3(a), but has significantly lower yields and a different distribution of yield nodes along the $k$-axis, at all harmonic orders. Above the interband-gap threshold, the interband yield has a rich $k$-dependent structure of even and odd harmonics. According to Eq. (37), the spectral range of interband high-order harmonics above the interband-gap threshold at the 11th harmonic is limited by ${\omega^{c}(k,\,A_{0})}$ . This $k$-dependent upper limit is indicated in Fig. 3(b) by the V-shaped solid black line. According to Eq. (45), contributions to the interband yield from $k$-channels that are $N$ orders of magnitude smaller than the maximal yield at the $Gamma$-point lie in the range ${\pm k_{max}(A_{0})/2}$. The V-shaped solid black line shows their onset for ${N=5}$.

We note that our numerical yields in Fig. 3, including contributions to HHG for $k$-channels within the entire first BZ, are incompatible with the assumption in previous studies Ghimire et al. (2012); Wu et al. (2015) that only a small part of the first BZ near the ${\Gamma-}$point contributes to HHG in solids. Even though for the one-dimensional model solid investigated here computing time is not an issue, the limit imposed by Eq. (45), and its numerical validation in Fig. 3, in addition to providing physical insight into the HHG process, is relevant for reducing the computational effort in multi-band HHG calculations based on a three-dimensional representation of the solid.

While the $k$-channel-resolved HHG spectrum in Fig. 3(c) includes even and odd high-order harmonics and depends in a rather complex way on the harmonic photon energy and lattice momentum $k$, including contributions to HHG from the entire first BZ according to Eq. (15), results in the comparatively simple total HHG spectrum shown in Fig. 3(d). As expected due to the inversion symmetry about the $\Gamma-$point of the Kronig-Penney band structure (Fig. 1), the total HHG spectrum in Fig. 3(d) is strongly dominated by odd harmonics.

III.1.2 Field-strength dependence of HHG spectra

Figure 5 shows the total HHG yield (including intra- and interband HHG) for the ${k=0}$ channel, i.e., only including the ${\Gamma-}$point at the center of the first BZ, for peak electric-field strengths ${0.05\,\text{V/\AA}<E_{0}<0.3\,\text{V/\AA}}$. This field-strength range corresponds to the laser-peak-intensity range ${3.3\times 10^{10}\,\text{W/cm}^{2}<\,I_{0}\,<1.19\times 10^{12}\,\text{W/cm}^{2}}$. We chose the upper limit of this interval to slightly exceed the field strength ${(\pi/a)\omega_{0}^{-1}\approx 0.283\,\text{V/\AA}}$. At this field strength, the vector-potential amplitude is ${A_{0}=\pi/a}$, such that the field-dressed time-dependent crystal momentum ${\kappa(t)=k+A(t)=A(t)}$ defined in Eq. (9) explores the entire first BZ (Fig. 1) within one optical cycle in the plateau of the laser pulse and, thus, the entire range of local band gaps. We selected the lower limit of the field strength in view of the approximations made in Sec. II.2, which are based on a series expansion in the parameter ${\gamma=\omega_{0}\sqrt{2\Delta\varepsilon_{0}^{cv}m^{*}_{0}}/E_{0}}$. Requesting ${\gamma<1}$, implies for MgO ${E_{0}>0.1\,\text{V/\AA}}$, slightly above the lower limit of the assumed range of field strengths.

The comparison of the intra- and interband yields in Figs. 5(a) and  5(b) shows that the intraband emission dominates the total yield in Fig. 5(c) below and interband emission above the band-gap threshold near the 11th harmonic, for the entire considered range of electric-field strengths. The HHG cutoff is thus determined by interband emission and displayed as superimposed white lines in Figs. 5(b) and  5(c). We analyze the HHG cutoff behavior in more detail in Sec. III.1.3 below. Over the range of displayed electric-field strengths the intra- and interband spectra are dominated by odd harmonics, with slim traces of even harmonics. Small even and non-integer HHG yields were also noticed in an SBE-based calculation by Li et al.  Li et al. (2018b) and explained in terms of the combined effect of the external-field and time dependence of the TDM $D^{vc}_{\kappa(t)}$. We speculate that this effect may be enhanced due to our inclusion of the term defined by Eq. (19) in Eq. (18). This term is absent in the SBE model.

Relaxing the ${k=0}$-channel ($\Gamma$-point) emission restriction, Fig. 6 shows yields obtained after integrating $k$ over the first BZ. As for the $\Gamma$-point-only yield in Fig. 5, the comparison of the intra-, interband, and total yields in Figs. 6(a), 6(b), and 6(c), respectively, reveals for the entire shown laser-peak-intensity range that below the lowest band-gap threshold, near the 11th harmonic, the yields are dominated by intraband emission, while above this threshold practically only interband emission occurs. Even though the spectra include traces of even and non-integer harmonics (as they do for $\Gamma-$point-only emission), the contrast between odd and even harmonic yields is larger than for $\Gamma$-point-only emission.

As discussed in the preceding Sec. III.1.1, intraband emission at the $\Gamma-$point is directly related to the valence-band group velocity. It yields intense odd and even harmonics, as seen in Figs. 3 and 4. On the other hand, as given by Eq. (27), upon integration over the first BZ the valence-band term in Eq. (25) cancels (by symmetry), and high-yield emission requires the laser-electric field to be strong enough to effectively promote electrons to the conduction band. This explains the low yield at low intensities in Figs. 6(b) and 6(c), as compared to the corresponding $\Gamma-$point-only yields in Fig. 5. The experimental investigation of HHG below the lowest band-gap threshold and for low to moderate peak intensities might thus resolve the range ${k_{max}(A_{0})}$ of carriers involved in HHG for a specific substrate. However, it remains to be explored to what extent experimental focal-volume effects, i.e., averaging of the laser intensity profile, prevent the accurate determination of ${k_{max}(A_{0})}$.

III.1.3 Field-strength dependence of the HHG cutoff

In this subsection we analyze the field-strength dependence of the HHG cutoff frequency obtained including both, $\Gamma-$point-only emission and $k$-channels from the entire first BZ. As the comparison of interband and total yields in Figs. 5 and 6 reveals, inclusion of the entire first BZ increases the highest generated frequencies, as compared to $\Gamma-$point-only emission. For each given peak electric-field strength, we visually determine the cutoff as the HHG order at which the yield starts to decline exponentially (as indicated in Fig. 3(d)), for both $\Gamma-$point-only and BZ-integrated yields. This leads to the intensity-dependent cutoff shown by the markers in Fig. 7. For the shown range of harmonic orders, inclusion of all $k$-channels results in cutoffs (indicated as triangular markers) in good agreement with our saddle-point prediction (green dotted line) that lie 8 - 10 harmonic orders above the cutoff for $\Gamma-$point-only emission (square markers). The cutoff orders predicted for the $\Gamma$-point-only emission by our 2-band TDSE calculations and our saddle-point analysis (red dashed line) are in fair agreement with the 2-band TDSE results of Wu et al.  Wu et al. (2015) (solid blue blue line), which we adapted from Fig. 3(b) in Ref. Wu et al. (2015). In agreement with Wu et al. , we find that the cutoff increases approximately linearly with laser peak-electric-field strength over the displayed range of harmonic orders, albeit with a noticeably smaller slope.

To better understand the field-strength dependence of the HHG cutoff, we resort to the saddle-point analysis of the HHG process in the 2-band approximation discussed in Secs. II.3.1 and II.3.3. Referring to HHG including all $k$-channels in the first BZ, we numerically solve the saddle-point Eq. (50) for $t_{e}$ and subsequently Eq. (49) for each $k_{r}$. This yields the frequency curve $\Delta\varepsilon_{\kappa_{r}(t_{e})}^{cv}$, from which we get the maximal frequency $\omega^{c}_{BZ}$ and cutoff harmonic order $\omega^{c}_{BZ}/\omega_{0}$. The red line in Fig. 8, shows the harmonic energy ${\Delta\varepsilon_{\kappa_{r}(t_{e})}^{cv}}$ as a function of values for $k_{r}$ that contribute to the BZ-integrated interband yield for a peak field strength of the driving laser of ${0.13\text{V/\AA}}$. As discussed in Sec. II.3.3, ${|k_{r}|<A_{0}}$ (indicated by vertical thin blue lines). The maximum harmonic energy is obtained at crystal momenta ${\pm k_{r}^{c}}$ indicated by the dotted green lines. $k$-channels with ${k=\pm k_{r}^{c}}$ thus yield the maximum cutoff energy at the given laser field strength. The V-shaped black solid line in Fig. 8 shows the maximum energy ${\omega^{c}(k,A_{0})}$, i.e., the electron-hole-pair-recombination energy in each $k-$channel given by Eq. (37) [cf., Fig. 3(b)]. This energy depends on the maximum local band gap in each channel and exceeds the cutoff energy in the BZ-integrated yield $\Delta\varepsilon_{\kappa_{r}(t_{e})}^{cv}$ (red line). Since $\omega^{c}(k_{r},A_{0})\approx\Delta\varepsilon_{\kappa_{r}(t_{e})}^{cv}$ the ${k=\pm k_{r}^{c}}$ channels have the largest contribution the HHG yield.

III.2 Multiband spectra

III.2.1 Field-strength dependence of HHG spectra

Figures 9 and 10 show multiband HHG spectra for both the $\Gamma$-point-only and BZ-integrated calculations, respectively. Including the lowest 13 bands of the Kronig-Penney model MgO crystal, we find these spectra to have converged. The convergence of the results was determined from the evaluation, in the $k=0$ channel, of the population $\rho^{MM}_{0}(t)$ in each new band $M$ we added to the calculation. Based on our convergence criterium $\rho^{MM}_{0}(t)<10^{-12}\rho^{vv}_{0}(0)$, we find $M=13$ for a maximum peak field strength of ${E_{0}=0.3\,\text{V/\AA}}$. We included the same number (13) of bands for calculations at lower field strengths and for $k\neq 0$. Comparison with the two-band yields of Figs. 5 and 6 shows that the spectral region below the first band-gap threshold is dominated by the electron dynamics in the lowest two bands. For harmonic orders below $\approx 35$, the shape of the multiband spectra and their laser-electric-field dependence largely resemble the two-band spectra. Our numerical tests showed that the inclusion of more than two bands gradually improves the agreement with fully converged spectra in the shown spectral range.

Above the $\approx$ 35th harmonic and at the highest field strengths shown in Figs. 9 and 10, a new plateau emerges, which we attribute to contributions from the second and third conduction band. If only the ${k=0}$ channel is included, the second plateau emerges at a higher field strength, $E_{0}>\pi\omega_{0}/a$, than in calculations including the entire first BZ. This field strength corresponds to the vector potential for which the cutoff frequency acquires its maximum value ${\omega^{c}(0,\,E_{0}/\omega_{0})=\Delta\varepsilon^{vc}_{max}}$ given by Eq. (37).

III.2.2 Field-strength dependence of the HHG cutoff

Figure 11 shows the field-strength dependence of the HHG cutoff we obtained from 13-band TDSE calculations by either including $\Gamma-$point-only emission (square markers) or $k$-channels from the entire first BZ (triangular markers), in comparison with the multi-band calculations from the literature. The solid blue line shows 51-band TDSE results for emission from the $\Gamma-$point-only, adapted from Fig. 3(a) in Ref. Wu et al. (2015). The orange dashed line shows HHG cutoff orders obtained by the direct propagation of the one-dimensional TDSE for $\Gamma-$point-only emission by Li et al.  Li et al. (2018b). The black dash-dotted line is adapted from the 6-band SBE calculation, including $k$-channels from the entire first BZ, of Li et al. .

For $\Gamma-$point only emission yields our stationary-phase approximation results in Fig. 11 (red dashed line) compare well with our 13-band calculations (square markers). Both are in reasonable agreement with the 51-band TDSE calculation of Wu et al.  Wu et al. (2015) (blue line). The cutoff harmonic orders we obtain in stationary-phase approximation are red-shifted by about three harmonic orders relative to the results of Li et al.  Li et al. (2018b) (orange dashed line), but approximately match the slope of the field-strength-dependent cutoff increase found by Li et al. . The cutoff harmonic orders predicted by our BZ-integrated 13-band calculation (triangular markers) agree well with our analytical stationary-phase approximation (green dotted line) and the 6-band SBE calculation of Li et al. (black dashed-dotted line). We note that the calculation of Li et al. includes a heuristic dephasing time, in contrast to our approach, where, apart from the adjusted potential strenghth of the Kronig-Penney model potential, no ad hoc parameters are introduced.

IV Conclusions

We investigated intra- and interband HHG from a solid within the SAE approximation, based on a one-dimensional Kronig-Penney model of the solid’s electronic band structure. We expanded the active electron wavefunction in a basis of adiabatically field-dressed Bloch states (Houston states). We theoretically analyzed within a stationary-phase approximation and numerically evaluated in two-band calculations, contribution to the HHG yield from specific initial-state crystal momenta $k$ in the first BZ. This revealed essential contributions to the HHG yield from non-zero crystal momenta in the first BZ. We found complementary contributions to the HHG yield from intra- and interband emission. While intraband emission dominates below the threshold for interband excitation, interband emission almost exclusively determines the HHG yield above this threshold. Applying a stationary-phase approximation, we derived closed analytical expressions for the HHG yields and cutoff dependence on the laser-electric-field strength. This allow us, for example, to estimate the loss of accuracy in calculating HHG yields that are induced by restricting $k$ to a small interval near the $\Gamma-$point in the first BZ, as compared to BZ-integrated yields.

By adjusting the Kronig-Penney model potential to the published electronic-structure data for MgO crystals, we validated our theoretical analysis against numerically calculated $k$-resolved and $k$-integrated HHG yields and cutoff frequencies for a range of laser intensities. These analytical predictions and numerical calculations are in very good agreement with each other. They are also in qualitative and fair quantitative agreement with different numerical calculations from other authors. Investigating the characteristics of HHG for initial crystal momenta in the entire first BZ, we find significant HHG yields at even high-order harmonic orders for ${k\neq 0}$ (off the $\Gamma-$point). Even harmonics are particularly prominent for the lowest intraband harmonics. As observed experimentally and expected due to the symmetric dispersion of the model MgO solid with respect to the BZ-zone center, our BZ-integrated yields predominantly contain odd harmonics. In contrast to HHG in gaseous atoms, we confirm the cutoff harmonic order to increase linearly with the peak electric-field strength of the driving laser.

Our numerical calculations for MgO confirm our stationary-phase-approximation results that restraining the crystal momentum to an arbitrarily small range of crystal momenta close to the $\Gamma-$point does not adequately model HHG in solids driven by intense laser fields. Strongly depending on the laser intensity, in calculations with more than two bands, we find a second plateau above the cutoff HHG order we determined in two-band calculations.

Appendix A Kronig-Penney basis functions

The band structure of the Kronig-Penney model for periodic delta potentials,

[TABLE]

is given by Eq. (3) of the main text and the eigenfunctions

[TABLE]

with the normalization factor

[TABLE]

and corresponding eigenenergies

[TABLE]

For crystal momenta ${k=0}$ and even $n$, and for ${k=\pi/a}$ and odd $n$, we have ${\alpha_{nk}=n\pi/a}$.

We adjust the global phase of the eigenfunctions to yield real matrix elements in Eq. (7) of the main text,

[TABLE]

for ${n\neq n^{\prime}}$, and real diagonal matrix elements

[TABLE]

Appendix B Analytical continuation of the Kronig-Penney basis

In our analysis of the interband current, we extend the energies in Eq. (57) to the complex ${K-}$plane, defining ${K=k+ik_{i}}$, where ${k_{i}}$ is the imaginary part of the complex crystal momentum, and request

[TABLE]

Equation (3) and its $K$-derivative at ${K=0}$yield

[TABLE]

After eliminating ${V_{0}}$ in these two equations, we see that the term in square brackets of the second equation cannot be zero, since ${2\pi\leq a\alpha_{n0}<3\pi}$ for ${n=v\,,c}$. This implies

[TABLE]

and shows that

[TABLE]

From the second derivative of Eq. (3) with respect to $K$ and Eq. (58), we obtain

[TABLE]

As noted in the main text, the coefficients in front of the second derivative are real. They do not vanish, since $2\pi\leq a\alpha_{n0}<3\pi$ for ${n=v\,,c}$. From this we conclude that

[TABLE]

and, consequently,

[TABLE]

Appendix C Alternative expression for the interband current

We here derive the interband current in Eq. (18) of in the main text. Defining the function

[TABLE]

we employ Eq. (6) to obtain its derivative

[TABLE]

Using ${\left(D^{vc}_{\kappa(t)}\right)^{2}=-\left|D^{vc}_{\kappa(t)}\right|^{2}}$ and substitution of ${\dot{\rho}^{cv}_{k}(t)}$ from Eq. (21) results in

[TABLE]

Since the second term in this equation is purely imaginary, after replacing the ${\eta^{er}_{k}(t)}$ from Eq. (60), Eq. (17) can be written as

[TABLE]

Appendix D Saddle-point approximation for the $k$-channel interband yield

We evaluate the action in Eq. (36) by splitting its integral representation into two parts,

[TABLE]

with ${t_{0}=\Re[t_{s}]}$. Along the integration path illustrated in Fig. 12, the first integral is real.

To evaluate the second integral, we need to calculate ${t_{s}}$. For this purpose, we consider the continuum-wave vector potential ${A_{0}\sin(\omega_{0}t)}$, for which we find

[TABLE]

with

[TABLE]

and ${E_{0}=A_{0}\,\omega_{0}}$. For ${\gamma<1}$ and ${k<A_{0}<\pi/a}$, Taylor expansion of Eq. (63) results in

[TABLE]

which will be a fair approximation if ${\sqrt{1-(k/A_{0})^{2}}<\gamma}$.

We can now carry out the integrations in Eq. (62) and obtain

[TABLE]

where

[TABLE]

To find the pre-exponential factor in Eq. (36), we first point out that the cross derivatives ${\partial^{2}\sigma_{\omega}(k,\,t,\,t^{\prime})/\partial t\partial t^{\prime}}$ and ${\partial^{2}\sigma_{\omega}(k,\,t,\,t^{\prime})/\partial t^{\prime}\partial t}$ vanish, such that

[TABLE]

Following Milos̆evic et al. (2002); Lewenstein et al. (1994), only keeping the regular solution in Eq. (64), ${\Im[t_{s}]>0}$, we obtain

[TABLE]

The interband HHG yield is thus given by

[TABLE]

In order to get a more explicit equation in terms of the band and laser field parameters, we approximate ${E(t_{s})}$ to first order in $\gamma$,

[TABLE]

and apply the approximations used to derive ${\Im[S(k,\,t_{s})]}$ in Eq. (65). This leads to the interband yield

[TABLE]

Appendix E Simplification of Eq. (48)

We here approximate Eq. (48), proceeding in a similar way as in Appendix D, by splitting Eq. (48) in a real and a complex integral (Fig. 12)

[TABLE]

with

[TABLE]

Following the derivation of Eq. (38), we find

[TABLE]

Using ${\kappa(t)=k+A(t)}$, we perform the integral in $I_{1}$ by contour integration in the complex plane to obtain

[TABLE]

where we used that, according Eq. (38), ${K(t_{s})}$ is imaginary. Thus, to second order in $\gamma$, ${I_{1}\approx 0}$ and Eq. (48) can be approximated as the real integral

[TABLE]

This expression is more suitable for numerical calculations than Eq. (48).

Acknowledgements

This work was supported in part by the Air Force Office of Scientific Research under award number FA9550-17-1-0369, NSF Grants No. 1464417 and No. 1802085, the “High Field Initiative” (CZ.02.1.01/0.0/0.0/15_003/0000449) of European Regional Development Fund (HIFI), and the project “Advanced research using high intensity laser produced photons and particles” (CZ.02.1.01/0.0/0.0/16_019/0000789) through the European Regional Development Fund (ADONIS). Any opinions, finding, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the United States Air Force.