# Microscopic analysis of high harmonic generation in semiconductors with
degenerate bands

Le Huu Thong  Ho Chi Minh City Institute of Physics, Vietnam Academy of
Science and Technology, 1 Mac Dinh Chi Street, District 1, Ho Chi Minh City,
Vietnam  Department of Physics, University of Science, Vietnam National
University Ho Chi Minh City, 227 Nguyen Van Cu Street, District 5, Ho Chi Minh
City, Vietnam  Cong Ngo  Department of Physics and CeOPP, Universität
Paderborn, Warburger Strasse 100, D-33098 Paderborn, Germany  Huynh Thanh Duc
Ho Chi Minh City Institute of Physics, Vietnam Academy of Science and
Technology, 1 Mac Dinh Chi Street, District 1, Ho Chi Minh City, Vietnam
Xiaohong Song  Research Center for Advanced Optics and Photoelectronics,
Department of Physics, College of Science, Shantou University, Shantou,
Guangdong 515063, China  Torsten Meier  Department of Physics and CeOPP,
Universität Paderborn, Warburger Strasse 100, D-33098 Paderborn, Germany

(August 27, 2024)

###### Abstract

Based on the multiband semiconductor Bloch equations a microscopic approach to
high-harmonic generation in crystalline solids which is able to properly
describe degenerate bands and band crossings is presented and analyzed. It is
well-known that numerical band structure calculations typically provide
electronic wave functions with an undetermined k-dependent phase which results
in matrix elements which contain arbitrary k-dependent phases. In addition,
such approaches usually mix degenerate bands and bands with an energy
difference smaller than the numerical precision in an arbitrary way for each
point in k-space. These ambiguities are problematic if one considers the
dynamics induced by electric fields since the matrix elements of the position
operator involve a derivative of the wave functions with respect to k. When
the light-matter interaction is described in the length gauge, the problem of
arbitrary phases and degenerate subspace mixing of Bloch states is solved by
adopting a smooth gauge transformation along the field direction. The results
obtained within this method are validated by comparing with calculations in
the velocity gauge. Although we obtain in both gauges the same overall result,
the length gauge is advantageous since it converges with a smaller number of
bands and thus requires significantly less numerical effort than the velocity
gauge. Also an unique distinction between inter- and intraband contributions
and thus an instructive physical interpretation is possible in the length
gauge whereas in the velocity gauge this is unclear. The computed
polarization-direction-dependent high-harmonic spectra agree well with
experimental data reported for GaAs. Furthermore, it is demonstrated that,
under proper conditions, the Berry curvature is largely responsible for the
even-order harmonics which are polarized perpendicular to the driving field.

##  I Introduction

Since its discovery in 2010 [1], high-harmonic generation (HHG) from solids
has been extensively studied. Under the excitation by an intense laser field,
solid crystals can emit high-order harmonics of the driving frequency, over a
very broad spectrum and with promising properties for applications. While the
role of different mechanisms are still being discussed, the semiconductor
Bloch equations (SBE) serve as a generic versatile approach that yields good
agreement with several experiments and allows for analyzing fine details of
the generation process and the dynamics of the photoexcitations [2, 3, 4, 5,
6, 7, 8].

The interaction between matter and light fields is usually described in either
the length gauge (LG) or the velocity gauge (VG). Early models for HHG in
solids were implemented exclusively in the LG. In this gauge, the HHG spectrum
is usually analyzed in terms of interband and intraband components, which are
believed to originate from optical transition between different bands and the
electrons’ acceleration within each band, respectively. Such an intuitive
picture has been shown to be reliable in the LG, it is, however, questionable
for the VG [9].

However, in the LG the SBE have the disadvantage that they requires the Bloch-
state basis to be smooth in reciprocal space (k-space). When the Bloch states
at each k-point are obtained from numerical diagonalization they are prone to
suffer from random phases and swapping of the band index of degenerate states.
This results in phase jumps of the complex transition dipole
$\bm{\xi}(\mathbf{k})$ between neighboring $\mathbf{k}$ vectors. A primitive
remedy to this problem was to only consider the absolute value
$\absolutevalue{\bm{\xi}(\mathbf{k})}$, which is an often used approximation
when solving the SBE [5]. Recently, the transition-dipole’s phase has been
claimed to carry information on the crystal’s symmetry and to be essential for
the description of even-order harmonics [10]. To incorporate the transition-
dipole phases properly when solving the SBE, it is required to implement a
gauge which results in Bloch functions which vary smoothly in k-space. For
non-degenerate energy bands, methods to fix the phases have been proposed and
applied to investigate HHG [11, 12, 13, 14]. In the presence of degenerate and
crossing bands, this issue has, however, rarely been addressed, since the
Bloch states are further subject to an arbitrary unitary transformation within
each degenerate subspace. Here we apply the method proposed by Virk and Sipe
[15] to construct a smooth parallel-transport gauge of degenerate Bloch
functions which allows us to solve the SBE in the LG and compute HHG emission
spectra.

Several recent publications [16, 17, 14] have demonstrated that HHG
calculations in both the VG and the LG agree if in both gauges the relevant
band are included and no further approximations, in particular, no
approximations that are not gauge invariant are made. Here, we confirm the
gauge independence of our results which verifies the correctness of our gauge
procedure. Even though we obtain in both gauges the same overall results, the
LG converges with a significantly smaller number of bands (less than half)
than required in the VG and thus the numerical effort of the calculations is
strongly reduced. Furthermore, unlike in the VG, in the LG a unique
distinction between inter- and intraband contributions and thus an instructive
physical interpretation is feasible and physically sound.

On the other hand, in HHG experiments it has been reported that there exist
components whose polarizations are perpendicular to that of a linearly-
polarized driving electric field. This emission has been ascribed to originate
from the crystal’s Berry curvature via the anomalous-velocity formula within a
semi-classical approximation [18, 19]. For much weaker THz fields that follow
an optical excitation with circularly-polarized light, similar Berry-
curvature-induced so called anomalous currents that flow in a perpendicular
direction have been proposed and measured in semiconductor quantum wells [20,
21]. In this work, we study that perpendicular HHG component within the
microscopic theory of the SBE. We find that the Berry curvature indeed
contributes predominantly to the perpendicular emission calculated within the
microscopic theory, though only when the laser’s frequency is small compared
to the band-gap and the field strength is strong enough to drive electrons to
regions of large Berry curvature.

The paper is organized as follows. In Section II.1 we present the SBE in both
the LG and the VG. Section II.2 describes the method of constructing a proper
parallel-transport gauge for linearly polarized incident fields including
degenerate Bloch states. Numerical results are presented and discussed in
Section III where we consider GaAs whose band structure and matrix elements
are obtained from the 30-band kp method [22]. In the results section we focus
on three aspects: the equivalence between calculations in the LG and in the
VG, the comparison with a reported HHG experiment performed on GaAs [23] and
the contribution of the Berry curvature to perpendicularly polarized HHG
emission. We close with a brief summary presented in Section IV.

##  II Theoretical methods

In this section we first present the derivation of the SBE in LG and VG and
then describe the parallel-transport gauge which is applied to be able to
numerically solve the SBE in LG.

###  II.1 The SBEs in LG and VG

In second-quantization the Hamiltonian describing Bloch electrons interacting
with a light field is written in LG as

$\displaystyle H^{\mathrm{}}$ $\displaystyle=$ $\displaystyle
H_{0}-e\mathbf{E}(t)\cdot\sum_{\lambda\lambda^{\prime}\mathbf{k}\mathbf{k^{\prime}}}\mathbf{r}_{\lambda\lambda^{\prime}}(\mathbf{k},\mathbf{k^{\prime}})a^{\dagger}_{\lambda\mathbf{k}}a_{\lambda^{\prime}\mathbf{k^{\prime}}},$
(1)

where $H_{0}$ is the Hamiltonian of the crystal,
$a^{\dagger}_{\lambda\mathbf{k}}$ $(a_{\lambda\mathbf{k}})$ is the creation
(annihilation) operator of an electron with wave vector $\mathbf{k}$ in band
$\lambda$, $\mathbf{E}(t)$ is the electric field (whose spatial dependence is
neglected here in the long-wavelength limit), and $\mathbf{r}$ is the position
operator. Although its explicit form is not required for the following
derivations, we would like to mention that in our numerical evaluations we
take $H_{0}$ as the single-particle Hamiltonian of the Bloch electrons, i.e.,
$H_{0}=\sum_{\lambda\mathbf{k}}\epsilon_{\lambda}(\mathbf{k})a^{\dagger}_{\lambda\mathbf{k}}a_{\lambda\mathbf{k}}$,
where $\epsilon_{\lambda}(\mathbf{k})$ is the band structure.

In the Bloch basis, the position operator can be represented as [24]

$\displaystyle\mathbf{r}_{\lambda\lambda^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})=\left(i\delta_{\lambda\lambda^{\prime}}\nabla_{\mathbf{k}}+\bm{\xi}_{\lambda\lambda^{\prime}}(\mathbf{k})\right)\delta({\mathbf{k}-\mathbf{k}^{\prime}}),$
(2)

where
$\bm{\xi}_{\lambda\lambda^{\prime}}(\mathbf{k})=-i\braket{\nabla_{\mathbf{k}}u_{\lambda\mathbf{k}}}{u_{\lambda^{\prime}\mathbf{k}}}$
is the transition dipole matrix element and $\ket{u_{\lambda\mathbf{k}}}$ is
the periodic part of the Bloch function. While one has the gauge freedom in
choosing the phase of $\ket{u_{\lambda\mathbf{k}}}$, this representation of
$\mathbf{r}$ is limited to gauges that smoothen the $k$-dependence of the wave
functions, such that the above derivative with respect to $k$ is well-defined.

The semiconductor Bloch equations [25], that describe the dynamical
optoelectronic response, can be expressed in terms of the reduced density
matrix elements
$\rho_{\lambda\lambda^{\prime}}(\mathbf{k})=\braket{a^{\dagger}_{\lambda^{\prime}\mathbf{k}}a_{\lambda\mathbf{k}}}$
as

$\displaystyle\frac{d\rho(\mathbf{k})}{dt}$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\commutator{H_{0}(\mathbf{k})-e\mathbf{E}(t)\cdot\bm{\xi}(\mathbf{k})}{\rho(\mathbf{k})}$
(3)
$\displaystyle-\frac{e}{\hbar}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}\rho(\mathbf{k}),$

where
$H_{0}(\mathbf{k})=e^{-i\mathbf{k}\cdot\mathbf{r}}H_{0}e^{i\mathbf{k}\cdot\mathbf{r}}$.
One can see that the equations of motion at different $\mathbf{k}$ vectors are
coupled by the derivative term $\nabla_{\mathbf{k}}\rho(\mathbf{k})$, which
stem from the position operator in Eq. (2). It is this term that poses the
requirement of constructing a smooth gauge of the Bloch functions in order to
be able to solve the SBE. In many previous works this problem has been ignored
and it has been assumed that the transition dipoles are independent of $k$ or
have a simple $k$-dependence or often only their absolute values were
considered. The proper smooth gauge is, however, unavoidable if one wants to
work in the preferred LG and consider the full information in the form of
complex matrix elements which arise from band structure calculations like
$\mathbf{k}\cdot\mathbf{p}$ or density functional theory that contain
arbitrary $k$-dependent phases arising from the numerical diagonalization.

In the VG the SBE read

$\displaystyle\frac{d\rho(\mathbf{k})}{dt}=-\frac{i}{\hbar}\commutator{H_{0}(\mathbf{k})-\frac{e}{m}\mathbf{A}(t)\cdot\mathbf{p}(\mathbf{k})}{\rho(\mathbf{k})},$
(4)

where $\mathbf{A}(t)=-\int_{-\infty}^{t}\mathbf{E}(t^{\prime})dt^{\prime}$ is
the vector potential and
$\mathbf{p}(\mathbf{k})=\dfrac{m}{\hbar}\nabla_{\mathbf{k}}H_{0}(\mathbf{k})$
is the canonical momentum matrix. In contrast to Eq. (3), the equations for
the reduced density matrix in VG, Eq. (4), at different $\mathbf{k}$ vectors
are independent of each other, hence a gauge choice ensuring the smoothness of
the wave function in $k$-space is not necessary. As is shown below, whereas we
do obtain very similar final results in both gauges the LG is more efficient
since a smaller number of bands is required to obtain converged results and,
in addition, the interpretation in terms inter- and intraband contributions
seems to be more adequate. The transformation between the LG and the VG can be
performed by the unitary operator
$Q=e^{-ie\mathbf{A}(t)\cdot\mathbf{r}/\hbar}$[16, 14].

It is noted that the common expression for the decoherence term, which
describes the exponential decay of off-diagonal elements of the density
matrix,

$\left(\frac{d\rho(\mathbf{k})}{dt}\Big{\rvert}_{\text{decoh}}\right)_{\lambda\lambda^{\prime}}=-\frac{\rho_{\lambda\lambda^{\prime}}(\mathbf{k})}{T_{2}}(1-\delta_{\lambda\lambda^{\prime}}),$
(5)

is not gauge invariant and the gauge transformation of this term can only be
done by numerical calculations [14]. Instead of Eq. (5), we therefore use its
gauge-covariant form that was proposed in Ref. 17 and allows for a simple
transformation between these two gauges. The decoherence term is given in LG
and VG, respectively, by [17]

$\displaystyle\left.\frac{d\rho(\mathbf{k})}{dt}\right|_{\text{decoh}}$
$\displaystyle=$
$\displaystyle-\frac{1}{T_{2}E_{g}^{2}}\left[H_{0}(\mathbf{k}),\left[H_{0}(\mathbf{k}),\rho(\mathbf{k})\right]\right]$
(6)

and

$\displaystyle\left.\frac{d\rho(\mathbf{k})}{dt}\right|_{\text{decoh}}$
$\displaystyle=$
$\displaystyle-\frac{1}{T_{2}E_{g}^{2}}\left[H_{0}(\mathbf{k})-\frac{e}{m}\mathbf{A}(t)\cdot\mathbf{p}(\mathbf{k}),\right.$
(7)
$\displaystyle\left.\left[H_{0}(\mathbf{k})-\frac{e}{m}\mathbf{A}(t)\cdot\mathbf{p}(\mathbf{k}),\rho(\mathbf{k})\right]\right],$

where $T_{2}$ is the phenomenological decoherence time and $E_{g}$ is the
bandgap energy.

###  II.2 The parallel-transport gauge

In order to solve the the SBE in LG, we follow Virk and Sipe [15] and
construct a local gauge transformation of the Bloch functions. The procedure
is based on the $\mathbf{k}\cdot\mathbf{p}$ perturbation theory for the band
structure. To first order in $\Delta\mathbf{k}$, the perturbation theory
provides the relation between the eigenstates at $\mathbf{k}$ and
$\mathbf{k}+\Delta\mathbf{k}$ [15, 26]

$\ket{u_{\lambda\mathbf{k}+\Delta\mathbf{k}}}=\sum_{\nu\mu}\ket{u_{\mu\mathbf{k}}}(\delta_{\nu\mu}-i\Delta\mathbf{k}\cdot\bm{\xi}^{\text{ter}}_{\mu\nu}(\mathbf{k}))g_{\nu\lambda}({\mathbf{k},\mathbf{k}+\Delta\mathbf{k}}),$
(8)

where
$g(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})=\exp(-i\Delta\mathbf{k}\cdot\bm{\xi}^{\text{tra}}(\mathbf{k}))$
is an unitary matrix accounting for arbitrary phase factors. The transition
dipole matrix $\bm{\xi}(\mathbf{k})$ is separated into
$\bm{\xi}^{\text{tra}}(\mathbf{k})$ and $\bm{\xi}^{\text{ter}}(\mathbf{k})$
which contains matrix elements between connected states and between
disconnected states, respectively. Here, two states
$|u_{\lambda\mathbf{k}}\rangle$ and $|u_{\nu\mathbf{k}}\rangle$ are called
connected if they have the same energy at at least one $k$-point in the
Brillouin zone (BZ), otherwise they are called disconnected.

Introducing the overlap matrices between the two bases at different $k$-points
$S^{\text{tra}}_{\lambda\lambda^{\prime}}({\mathbf{k},\mathbf{k}^{\prime}})=\braket{u_{\lambda\mathbf{k}}}{u_{\lambda^{\prime}\mathbf{k}^{\prime}}}\Delta_{\lambda\lambda^{\prime}}$
and
$S^{\text{ter}}_{\lambda\lambda^{\prime}}({\mathbf{k},\mathbf{k}^{\prime}})=\braket{u_{\lambda\mathbf{k}}}{u_{\lambda^{\prime}\mathbf{k}^{\prime}}}\left(1-\Delta_{\lambda\lambda^{\prime}}\right)$,
where $\Delta_{\lambda\mu}=1$ if $\lambda=\mu$ or if two bands $\lambda$ and
$\mu$ are connected and $\Delta_{\lambda\mu}=0$ otherwise, the above
perturbative relation gives

$g(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})=S^{\text{tra}}({\mathbf{k},\mathbf{k}+\Delta\mathbf{k}})$
(9)

and

$-i\Delta\mathbf{k}\cdot\bm{\xi}^{\text{ter}}(\mathbf{k})g(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})=S^{\text{ter}}(\mathbf{k},\mathbf{k}+\Delta\mathbf{k}).$
(10)

Eqs. (9) and (10) allow one to calculate
$g(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})$ and
$\bm{\xi}^{\text{ter}}(\mathbf{k})$ from the overlap matrices which are
determined by the wave functions obtained from the band structure
calculations. To eliminate the relative phase of eigenstates between
$\mathbf{k}$ and $\mathbf{k}+\Delta\mathbf{k}$ one performs the gauge
transformation

$\displaystyle\ket{u_{\lambda\mathbf{k}+\Delta\mathbf{k}}}\longmapsto\sum_{\mu}g^{*}_{\lambda\mu}(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})\ket{u_{\mu\mathbf{k}+\Delta\mathbf{k}}}.$
(11)

Although in practical calculations $g(\mathbf{k},\mathbf{k}+\Delta\mathbf{k})$
is not exactly unitary, it can be made unitary by a singular value
decomposition (SVD) algorithm.

By extending Eq. (11) to a series of $k$-points on a straight line in
$j$-direction, with unit vector $\hat{\mathbf{e}}_{j}$, and starting from the
first point $\mathbf{k}_{0}$, one constructs a gauge transformation matrix

$\displaystyle
W^{j}(\mathbf{k}_{0},\mathbf{k})=g(\mathbf{k}_{0},\mathbf{k}_{0}+\hat{\bf{e}}_{j}\Delta
k)...g(\mathbf{k}-\hat{\bf{e}}_{j}\Delta k,\mathbf{k})$ (12)

that removes the arbitrary relative phase of the eigenstates between every two
nearest neighboring $k$-points, which makes the Bloch functions smooth and
differentiable with respect to $\mathbf{k}$

$\displaystyle\ket{u_{\lambda\mathbf{k}}}\longmapsto\sum_{\mu}W^{j*}_{\lambda\mu}(\mathbf{k}_{0},\mathbf{k})\ket{u_{\mu\mathbf{k}}}.$
(13)

Mathematically, the above method is equivalent to the optimal alignment
procedure presented in Ref. 27, section 3.6.

In the new basis of Bloch functions, the projection onto the $j$-direction of
$\bm{\xi}^{\text{tra}}(\mathbf{k})$ vanishes and the projection of
$\bm{\xi}^{\text{ter}}(\mathbf{k})$ given by

$\displaystyle\xi^{\text{ter}}_{j}({\mathbf{k}})=\frac{i}{2\Delta
k}W^{j}({\mathbf{k}_{0},\mathbf{k}})$ $\displaystyle\ \ \
\times\Big{[}S^{\text{ter}}(\mathbf{k},\mathbf{k}+\hat{\mathbf{e}}_{j}\Delta
k)S^{\text{tra}\dagger}({\mathbf{k},\mathbf{k}+\hat{\mathbf{e}}_{j}\Delta k})$
(14) $\displaystyle\ \ \
-S^{\text{ter}\dagger}({\mathbf{k}-\hat{\mathbf{e}}_{j}\Delta
k,\mathbf{k}})S^{\text{tra}}({\mathbf{k}-\hat{\mathbf{e}}_{j}\Delta
k,\mathbf{k}})\Big{]}W^{j\dagger}({\mathbf{k}_{0},\mathbf{k}})$

is a smooth function of $\mathbf{k}$. The above gauge of the Bloch functions
is called the parallel-transport gauge [28].

##  III HHG from GaAs: Results and discussion

By integrating the SBE in LG (3) or VG (4) for an one-dimensional $k$-grid
that is parallel to the polarization-direction of the electric field, we
obtain the dynamics of the density matrix $\rho(k,t)$ which determines the
time-dependent electric current density via

$\displaystyle\mathbf{J}(t)=e\sum_{\lambda\lambda^{\prime}k}\mathbf{v}_{\lambda\lambda^{\prime}}(k)\rho_{\lambda^{\prime}\lambda}(k,t),$
(15)

where $\mathbf{v}=\frac{i}{\hbar}\left[H,\mathbf{r}\right]$ is the velocity
operator. The momentum space representations of $\mathbf{v}$ are
$\mathbf{v}(k)=\mathbf{p}(k)/m$ in LG and
$\mathbf{v}(k)=\left(\mathbf{p}(k)-e\mathbf{A}(t)\right)/m$ in VG.

The photoexcited current oscillates rapidly and contains multiples of the
excitation frequency which correspond to high-harmonic radiation. From the
Fourier transform of current density we obtain the spectrum of the emission
intensity

$\displaystyle I_{\text{HHG}}(\omega)\propto\absolutevalue{J(\omega)}^{2}.$
(16)

In the following we carry out numerical calculations for bulk GaAs. The
electronic band structure and wave functions of GaAs for the entire Brillouin
zone (BZ) are obtained from a 30-band $\mathbf{k}\cdot\mathbf{p}$ model [22].
This quite sophisticated model includes spin-orbit coupling and allows to
describe the inversion asymmetry of GaAs crystal.

The electric field of the exciting laser pulse is described by

$\displaystyle\mathbf{E}(t)=E_{0}\hat{\mathbf{e}}\
e^{-2\ln(2)t^{2}/\tau^{2}}\sin(\omega_{0}t),$ (17)

where $\hat{\mathbf{e}}$ denotes the polarization direction, $E_{0}$ is the
maximal amplitude, $\tau$ is the pulse duration (FWHM of the Gaussian
envelope), and $\omega_{0}$ is the central light frequency.

###  III.1 Comparison between LG and VG

The computed time-dependent photocurrent for the two gauges and the resulting
emission spectra are presented in Fig. 1. In the numerical simulations, the
driving laser is linearly polarized in the $[100]$ crystallographic direction
($\Gamma\text{X}$ direction in the BZ) and has an amplitude of
$E_{0}=10\leavevmode\nobreak\ \text{MV/cm}$, $\tau=60\leavevmode\nobreak\
\text{fs}$, and a center frequency of
$\hbar\omega_{0}=0.38\leavevmode\nobreak\ \text{eV}$.

Figure 1: (a) shows the electric field of the time-dependent exciting THz
laser pulse. (b) and (c) show the temporal dynamics of the current density
calculated for a THz field that is linearly polarized in the $[100]$
crystallographic direction (corresponding to the $\Gamma\text{X}$ direction in
the BZ) in the LG and the VG, respectively. The total current (red solid
lines) is given by the sum of the intraband (blue dashed lines) and the
interband (orange dashed lines) currents. (c) Intensity spectra of the emitted
high harmonics in the two gauges.

Fig. 1(b) shows the dynamics of the current density calculated in LG. The red-
solid line represents the total current, while the blue and orange dashed
lines display its two components corresponding to the intraband (diagonal) and
interband (off-diagonal) contributions, respectively. While the interband
current is in phase with the electric field (grey line in Fig. 1(a)), the
intraband current lags behind the field by a phase of $\pi/2$. This phase lag
of the intraband current can be explained by the acceleration theorem for the
electron wave packet motion $\hbar\dot{\mathbf{k}}=e\mathbf{E}(t)$, i.e., the
change of $\mathbf{k}$ of proportional to the time integral over
$\mathbf{E}(t)$.

Fig. 1(c) is the same as Fig. 1(b) but calculated in the VG. Compared to the
LG, the interband and intraband currents in the VG are both much stronger and
are apart from their opposite sign very similar to each other. Obviously,
these currents are gauge dependent and have therefore in the VG no clear
physical significance [16, 9]. Due to the strong cancellations between the
interband and intraband currents in the VG, the resulting gauge-independent
total current is identical to that obtained in the LG.

We would like to emphasize that the calculation in the LG is already converged
when the 6 highest valence and the 8 lowest conduction bands are considered in
the SBE, while the VG requires that all 30 bands included. Since in the VG the
subcurrent components are gauge dependent, a full sum (over all available
bands) is necessary so that their gauge dependencies basically compensate each
other which results in a total current which is independent of the chosen
gauge (if sufficiently many bands can be considered). Since the effort for
numerical solutions of the SBE (without including many-body interactions)
scales quadratically with the number of bands, the evaluations in the LG are
significantly faster than those in the VG and in addition in the VG and also
the required computer memory is reduced similarly.

In Fig. 1(c) we show the spectral intensity of the emitted HHG computed in LG
and VG. Since the incident laser pulse is polarized in the $\Gamma\text{X}$
direction, the spectrum consists of only odd-order harmonics. Just like the
time-dependent total currents, the HHG spectra in the two gauges are almost
identical, i.e., the 30 band in our band-structure model are sufficient to
ensure convergence of the VG results for the considered excitation conditions.

###  III.2 Comparison to experiment

Next, we compare our numerical results with experimental data reported in Ref.
23. In our simulations, we consider the experimental conditions and use the
following parameters for the driving laser: a photon energy of
$\hbar\omega_{0}=E_{g}/4$, a pulse duration (FWHM) of 5 laser cycles, and a
maximal electric field amplitude of $E_{0}=10\leavevmode\nobreak\
\text{MV/cm}$. The only fitting parameter is the decoherence time $T_{2}$,
which is taken to be on the order of a few femtosecond as in several previous
HHG studies [4]. As shown recently, such very small dephasing times can be
justified as they yield similar HHG spectra as arising when including
propagation effects [29]. Both types of the decoherence term described in Sec.
II.1 can be used to match the experiment with nearly the same value of
$T_{2}$. Here, we show the calculated results in the LG using the decoherence
term Eq. (6) with $T_{2}=3$ fs (which corresponds to 3/10 of the laser
period).

Figure 2: Calculated HHG spectra of GaAs with the exciting laser field
polarized linearly along the (a) $\Gamma\text{K}$ and (b) $\Gamma\text{L}$
directions, respectively, in comparison to experimental data which is
extracted from Ref. 23.

For the case that the incident light is polarized linearly in the
$\Gamma\text{K}$ direction, both the calculated and the measured HHG spectra
exhibit only odd-order harmonic peaks, see Fig. 2(a), since the inversion
asymmetry of GaAs is not probed in this excitation configuration. Our
theoretical simulations describe the decrease of the intensity of the odd-
order harmonics with increasing harmonic order in very good agreement with the
experimental observations. An important point is that an optical excitation
polarized in the $\Gamma\text{K}$ direction does also produce harmonics that
are polarized perpendicularly to the direction of the incident laser field.
This perpendicular component of HHG was, however, not measured in [23]. We
will present our theoretical investigation of this component in the next
section.

When the incident light is polarized linearly in the $\Gamma\text{L}$
direction ($[111]$ crystallographic direction), besides the odd-order
harmonics, also even-order harmonics with somewhat weaker intensities arise,
see Fig. 2(b). The even-order harmonics originate from the inversion asymmetry
of GaAs crystal. In the $30$-band $\mathbf{k}\cdot\mathbf{p}$ model [22],
$P^{\prime}$ is the only parameter that models this asymmetry and consequently
the even-order harmonics disappear if one artificially sets $P^{\prime}=0$.
Compared to experiment, the calculated even-order harmonics are somewhat
weaker which is probably the case because the 30-band
$\mathbf{k}\cdot\mathbf{p}$ model underestimates the inversion asymmetry of
GaAs as it uses a smaller $P^{\prime}$ parameter than other
$\mathbf{k}\cdot\mathbf{p}$ models [30, 31], e.g., in the 14-band
$\mathbf{k}\cdot\mathbf{p}$ model this term is almost ten times larger. For
example, if we manually double the value of $P^{\prime}$ in the 30-band model
(not shown in figure), the amplitude of the even-order harmonics increases
accordingly and the calculated even-order peaks agree better with experiment.
This is one example showing that HHG spectra can be used to gain information
on the atomic structure as a detailed comparison between measurements and
calculations can be used to improve the band structure models [32, 33].

###  III.3 Anomalous perpendicular currents arising from Berry curvature

Although we solve the SBE here in just one $k$-space dimension, the current
density obtained from $\eqref{current}$ is still a three-dimensional quantity.
For example, when the laser field is polarized in the $[110]$
($\Gamma\text{K}$) direction, beside the main current that flows parallel to
the field direction, the photoexcited current also has a smaller perpendicular
component that flows in the $[001]$ ($\Gamma\text{X}$) direction. We evaluate
this perpendicular current component in two ways: directly from the
microscopic approach using Eq. (15) and indirectly from the perturbative
equation for the anomalous velocity (Berry curvature approach)

$\mathbf{J}_{\perp}^{a}(t)=-\frac{e^{2}}{\hbar}\mathbf{E}(t)\times\expectationvalue{\mathbf{\Omega}},$
(18)

where
$\expectationvalue{\mathbf{\Omega}}=\sum\limits_{\lambda\lambda^{\prime}k}\mathbf{\Omega}_{\lambda\lambda^{\prime}}(k)\rho_{\lambda^{\prime}\lambda}(k,t)$
is the average of the Berry curvature

$\mathbf{\Omega}_{\lambda\lambda^{\prime}}(k)=i\sum_{\mu}\bm{\xi}^{\text{ter}}_{\lambda\mu}(k)\times\bm{\xi}^{\text{ter}}_{\mu\lambda^{\prime}}(k).$
(19)

Figure 3: Spectra of the HHG radiation with perpendicular polarization
calculated by the microscopic and the Berry curvature approaches for an
excitation photon energy of (a) $\hbar\omega_{0}=E_{g}/10$ and (b)
$\hbar\omega_{0}=E_{g}/20$. The vertical red lines indicate the bandgap energy
$E_{g}$. The amplitudes of the THz fields are $E_{0}=8$ MV/cm in (a) and
$E_{0}=4$ MV/cm in (b) such that the ratio $eE_{0}/\hbar\omega_{0}$ remains
the same.

We solve the SBE in the LG using the decoherence term Eq. (5) with $T_{2}=1/5$
of the laser period. The perpendicularly polarized HHG spectra calculated by
the microscopic (blue line) and Berry curvature (orange line) approaches are
shown in Fig. 3(a) for a photon energy of $\hbar\omega_{0}=E_{g}/10$. We find
that the perpendicularly-polarized HHG component is dominated by even-order
harmonics and its intensity is about two orders smaller than that of the
parallel HHG component. The overall good agreement between two approaches
confirms the Berry curvature as the origin of the perpendicular current
component. Since $\mathbf{\Omega}(k)$ is an odd function of $k$, the anomalous
velocity has a fundamental frequency of $2\omega_{0}$, and hence the
perpendicular HHG basically contains only even-order harmonics [18, 19].

The results shown in Fig. 3(b) are obtained for a more slowly varying laser
with $\hbar\omega_{0}=E_{g}/20$. To ensure that the region which the
electronic wavepacket traverses to is the same as in the previous case, we
adjust the field amplitude so that $eE_{0}/\hbar\omega_{0}$ is unchanged. We
note that the expression for anomalous velocity, Eq. (18), is derived from the
first-order adiabatic perturbation theory. Thus, the difference between two
approaches describes the contributions beyond this approximation. Because the
adiabatic condition holds better for smaller laser frequencies, the agreement
between two approaches in Fig. 3(b) is closer than that in Fig. 3(a).

An intrinsic limitation of one-dimensional models is that they only allow to
compute the interband contribution to the perpendicularly-polarized current
but not the intraband one. Our result therefore do not include the
perpendicular odd-order harmonics which were shown to be produced
predominantly by the intraband current for fields with high amplitudes [34].

##  IV CONCLUSIONS

We present and analyze a microscopic approach to high-harmonic generation in
solids with degenerate bands and crossing band. We confirm that the
calculations in the LG and in the VG produce the same results if sufficiently
many bands are included in the numerical solutions and gauge-dependent
approximation, e.g., when including dephasing, are avoided.

To be able to solve the SBE in the LG we implement a parallel-transport gauge
which is able to properly treat degenerate bands and results in wave functions
and matrix elements which vary smoothly as function of k. To obtain converged
HHG results in the LG a bit less than half the number of bands is required
than in the VG. Thus the numerical effort is strongly reduced in the LG and,
furthermore, the LG allows to distinguish between inter- and intraband
contributions and thus provides an instructive physical interpretation.

Our calculated results for HHG from GaAs are in good agreement with recent
experimental data on the polarization-direction dependence. Furthermore, our
approach is able to describe the perpendicularly-polarized even-order
harmonics. These are caused by the Berry curvature and also have contributions
beyond the adiabatic approximation that are included in our microscopic
approach.

The approach introduced here is very general, can be combined with various
band structure computation methods, and is applicable to other strong field
phenomena also involving resonant optical fields.

###### Acknowledgements.

This work is funded by the Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under the grant No. 103.01-2017.42, the
Deutsche Forschungsgemeinschaft (DFG) through project ME 1916/4, the National
Natural Science Foundation of China under the grant No. 12074240, and by the
Sino-German Mobility Programme (Grant No. M-0031). We thank the PC2 (Paderborn
Center for Parallel Computing) for providing computing time.

## References

  * [1] S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nature Physics 7, 138 (2011).
  * [2] D. Golde, T. Meier, and S. W. Koch, Phys. Rev. B 77, 075330 (2008).
  * [3] D. Golde, M. Kira, T. Meier, and S. W. Koch, phys. stat. sol. (b) 248, 863 (2011).
  * [4] T. T. Luu and H. J. Wörner, Phys. Rev. B 94, 115164 (2016).
  * [5] G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Phys. Rev. Lett. 113, 073901 (2014).
  * [6] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, Nature Photonics 8, 119 (2014).
  * [7] M. Hohenleutner, F. Langer, O. Schubert, M. Knorr, U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature 523, 572 (2015).
  * [8] X. Song, R. Zuo, S. Yang, P. Li, T. Meier, and W. Yang, Optics Express 27, 2225 (2019).
  * [9] G. Ernotte, T. J. Hammond, and M. Taucer, Phys. Rev. B 98, 235202 (2018).
  * [10] S. Jiang, J. Chen, H. Wei, C. Yu, R. Lu, and C. D. Lin, Phys. Rev. Lett. 120, 253201 (2018).
  * [11] U. Lindefelt, H.-E. Nilsson, and M. Hjelm, Semicond. Sci. Technol. 19, 1061 (2004).
  * [12] M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B. Gaarde, Phys. Rev. A 91, 043839 (2015).
  * [13] J. Li, X. Zhang, S. Fu, Y. Feng, B. Hu, and H. Du, Phys. Rev. A 100, 043404 (2019).
  * [14] L. Yue and M. B. Gaarde, Phys. Rev. A 101, 053411 (2020).
  * [15] K. S. Virk and J. E. Sipe, Phys. Rev. B 76, 035213 (2007).
  * [16] P. Földi, Phys. Rev. B 96, 035112 (2017).
  * [17] M. S. Wismer and V. S. Yakovlev, Phys. Rev. B 97, 144302 (2018).
  * [18] H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, Nature Physics 13, 262 (2017).
  * [19] T. T. Luu and H. J. Wörner, Measurement of the Berry curvature of solids using high-harmonic spectroscopy, Nature Communications 9, 1 (2018).
  * [20] K. S. Virk and J. E. Sipe, Phys. Rev. Lett. 107, 120403 (2011).
  * [21] S. Priyadarshi, K. Pierz, and M. Bieler Phys. Rev. Lett. 115, 257401 (2015).
  * [22] S. Richard, F. Aniel, and G. Fishman, Phys. Rev. B 70, 235204 (2004).
  * [23] P. Xia, C. Kim, F. Lu, T. Kanai, H. Akiyama, J. Itatani, and N. Ishii, Optics Express 26, 29393 (2018).
  * [24] E. Blount, Solid State Physics: Advances in Research and Applications (Academic, New York, 1962), Vol. 13, pp. 305–373.
  * [25] H. Haug and S. W. Koch, Quantum theory of the optical and electronic properties of semiconductors, 4th ed. (World Scientific, Singapore, 2004).
  * [26] M. Lax, Symmetry Principles in Solid State and Molecular Physics (Dover, Mineola, 2001).
  * [27] D. Vanderbilt, Berry Phases in Electronic Structure Theory (Cambridge University Press, 2018).
  * [28] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * [29] I. Kilen, M. Kolesik, J. Hader, J. V. Moloney, U. Huttner, M. K. Hagen, and S. W. Koch, Phys. Rev. Lett. 125, 083901 (2020).
  * [30] H. Mayer and U. Rössler, Phys. Rev. B 44 9048 (1991).
  * [31] I. Saidi, S. Ben Radihia, and K. Boujdaria, J. Appl. Phys. 107, 043701 (2010).
  * [32] N. Tancogne-Dejean, O. D. Mücke, F. X. Kärtner, and A. Rubio, Phys. Rev. Lett. 118, 087403 (2017).
  * [33] A. A. Lanin, E. A. Stepanov, A. B. Fedotov, and A. M. Zheltikov, Optica 4, 516 (2017).
  * [34] K. Kaneshima, Y. Shinohara, K. Takeuchi, N. Ishii, K. Imasaka, T. Kaji, S. Ashihara, K. L. Ishikawa, and J. Itatani, Phys. Rev. Lett. 120, 243903 (2018).