# Linking High-Harmonic Generation and Strong-Field Ionization in Bulk
Crystals

Peter Jürgens<EMAIL_ADDRESS>Max-Born-Institute for Nonlinear Optics
and Short Pulse Spectroscopy, Max-Born-Str. 2A, 12489 Berlin, Germany
Advanced Research Center for Nanolithography (ARCNL), Science Park 106, 1098
XG Amsterdam, Netherlands  Sylvianne D.C. Roscam Abbing  Advanced Research
Center for Nanolithography (ARCNL), Science Park 106, 1098 XG Amsterdam,
Netherlands  Mark Mero  Max-Born-Institute for Nonlinear Optics and Short
Pulse Spectroscopy, Max-Born-Str. 2A, 12489 Berlin, Germany  Graham G. Brown
Max-Born-Institute for Nonlinear Optics and Short Pulse Spectroscopy, Max-
Born-Str. 2A, 12489 Berlin, Germany  Marc J.J. Vrakking  Max-Born-Institute
for Nonlinear Optics and Short Pulse Spectroscopy, Max-Born-Str. 2A, 12489
Berlin, Germany  Alexandre Mermillod-Blondin  Max-Born-Institute for Nonlinear
Optics and Short Pulse Spectroscopy, Max-Born-Str. 2A, 12489 Berlin, Germany
Peter M. Kraus  Advanced Research Center for Nanolithography (ARCNL), Science
Park 106, 1098 XG Amsterdam, Netherlands  Department of Physics and Astronomy,
and LaserLaB, Vrije Universiteit, De Boelelaan 1105, 1081 HV Amsterdam,
Netherlands  Anton Husakou  Max-Born-Institute for Nonlinear Optics and Short
Pulse Spectroscopy, Max-Born-Str. 2A, 12489 Berlin, Germany

###### Abstract

The generation of high-order harmonics in bulk solids subjected to intense
ultrashort laser pulses has opened up new avenues for research in extreme
nonlinear optics and light-matter interaction on sub-cycle timescales. Despite
significant advancement over the past decade, a complete understanding of the
involved phenomena is still lacking. High-harmonic generation in solids is
currently understood as arising from nonlinear intraband currents, interband
recollision and ionization-related phenomena. As all of these mechanisms
involve or rely upon laser-driven excitation we combine measurements of the
angular dependence of nonlinear absorption and high-order harmonic generation
in bulk crystals to demonstrate the relation between high-harmonic emission
and nonlinear, laser-induced ionization in solids. An unambiguous correlation
between the emission of harmonics and laser-induced ionization is found
experimentally, that is supported by numerical solutions of the semiconductor
Bloch equations and calculations of orientation-dependent ionization rates
using maximally localized Wannier-functions.

††preprint: APS/123-QED

##  I Introduction

High-order harmonic generation (HHG) in gases marked the birth of attosecond
science allowing to track carrier dynamics on sub-cycle timescales [1, 2, 3,
4]. Transferring the concept of HHG to solid-state systems [5] has opened up
several research avenues ranging from fundamental investigations of strong-
field-driven carrier dynamics [6] towards the development of compact extreme-
ultraviolet (XUV) sources [7] and petahertz electronics [8, 9]. In solids,
four main sources of nonlinearity have been identified and investigated with
regards to harmonic generation. First, traditional perturbative nonlinearities
[10] have been invoked to explain conventional second harmonic generation
(SHG) ever since Franken’s seminal work [11] that marked the advent of
nonlinear optics. It is based on the anharmonic motion of bound electrons in
the valence band and is held responsible for low-order harmonic generation
with multiple applications in frequency conversion as well as in optical-
parametric chirped pulse amplification (OPCPA) [12]. Second, high-order
harmonics extending to frequencies in the XUV [13] have been associated with
interband recollisions (in strong analogy with the three-step model in the
atomic case [14, 1, 2]). Third, high harmonics have also been explained by
intraband currents where the nonlinearity enters through the non-parabolic
shape of the conduction (and valence) bands [15]. Fourth, the nonlinearity
that is inherent to the process of photoionization was proposed as a possible
origin for harmonic generation [16, 17]. Three possible ionization-related
mechanisms have been discussed in the context of HHG [18, 19]. Brunel
harmonics, arising from the acceleration of quasi-free carriers excited by
photoionization have been analyzed in gases, clusters, bulk solids and thin
films [20, 21, 22, 23, 24]. Moreover, after non-resonant photoionization an
excited electron has nonzero velocity which provides sub-cycle contributions
to the polarization and gives rise to the generation of harmonics. However,
the strength of these velocity harmonics is expected to be insignificant under
typical experimental conditions due to the modest photon energies that result
in small excess velocities. Finally, the previously overlooked injection
current, originating from the spatial displacement of an electron wavepacket
during laser-driven tunneling ionization, was identified as the dominant
source of low-order harmonic generation in fused silica [25, 19].
The highly anisotropic angular dependence of the HHG process in crystalline
solids has been correlated to the electronic band structure [13], petahertz
photocurrents [26], real-space trajectories of electronic wavepackets [27] and
van-Hove singularities [28, 29]. Low-order (below-bandgap) harmonic emission
has been associated with intraband dynamics while above-bandgap HHG is
generally attributed to interband recollision [26, 30].
All of the aforementioned HHG mechanisms - except for the Kerr-type
nonlinearity - rely either on the presence or on the excitation of quasi-free
conduction band electrons. In the most common experimental scheme the photon
energy of the driving near-infrared or mid-infrared laser field $\hbar\omega$
is smaller than the bandgap $E_{g}$ of the irradiated crystal. Thus, excited
carriers are generated by nonlinear photoionization or electron-impact
ionization [31]. Despite intense research in the past decade, the connection
between carrier excitation and its impact on the anisotropic angular
dependence of the harmonic emission has not yet been fully understood. While
photoionization is often a key factor in harmonic generation, to date, there
has been no investigation of its angular dependence under HHG conditions.
In this article we experimentally demonstrate a correlation between
orientation-dependent ionization yields and the efficiency of HHG in various
bulk crystals. We combine measurements of the angular dependence of nonlinear
absorption (also referred to as multiphoton crystallography [32, 33]) and
angle-dependent HHG efficiencies to establish a direct link between the
excitation of quasi-free carriers and the emission of high-order harmonics. We
identify an unambiguous connection between intense HHG emission and strong
ionization for distinct orientations in various oxide and fluoride crystals.
Our experimental results are supported by numerical simulations based on the
semiconductor Bloch equations (SBEs) and calculations of the photoionization
rate based on the Wannier-formalism. Our numerical results reproduce the
experimentally observed angular dependence of the harmonic emission and
provide insights into the complex interplay of interband and intraband
dynamics as well as the importance of ionization-related mechanisms for low-
order harmonic generation.

##  II Experimental Setup

In the experiments we used the signal beam of a home-built, high-repetition-
rate, dual-beam OPCPA ($100\text{\,}\mathrm{kHz}$) with a central wavelength
of $1500\text{\,}\mathrm{nm}$ and a full-width half-maximum (FWHM) pulse
duration of $\sim$$50\text{\,}\mathrm{fs}$ (architecture based on the one
presented in Ref. [34]) to generate high-order harmonics from bulk crystals.
The linearly polarized beam was focused to an $1/{e^{2}}$ beam diameter of
$\sim$$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (measured by the knife-edge
technique in ambient air) providing a maximum peak intensity in ambient air of
$\sim$$30\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ in the focal plane.

Figure 1: Experimental setup and results obtained in a bulk MgO crystal. (a).
A short-wavelength infrared (SWIR) laser pulse ($\lambda_{0}=\,$
$1500\text{\,}\mathrm{nm}$, $\tau=\,$$50\text{\,}\mathrm{fs}$) is focused into
the bulk of a crystal mounted on a rotation stage. The harmonic radiation is
analyzed using a visible-ultraviolet (VIS/UV) spectrometer. The transmitted
fundamental radiation can be analyzed by inserting a moveable dielectric
mirror directing the SWIR beam onto a photodiode. (b). Illustrative harmonic
spectra containing the third, the fifth and the seventh harmonic of the SWIR
driving field obtained in a bulk MgO crystal at two different orientations
($\theta=\,$$25\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$) and a peak intensity of
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. The large spectral width
of the seventh harmonic is attributed to spectral broadening in the detection
system. (c). Absorption of the fundamental laser field for different peak
intensities below the damage threshold as a function of the orientation angle
$\theta$.

High harmonics emitted from bulk crystals were analyzed in transmission
geometry with the help of a commercial visible-ultraviolet (VIS-UV)
spectrometer (Avantes AvaSpec-HS1024x58/122TEC). The energy remaining in the
short-wave infrared (SWIR) signal beam was directed onto a photodiode [cf.
Fig. 1(a)] with the help of a highly reflective mirror at
$1500\text{\,}\mathrm{nm}$. A possible contamination of the transmitted signal
beam by harmonic radiation was excluded due to the low reflectivity of the
dielectric mirror for the harmonic wavelengths and by the low conversion
efficiencies (typically around $10^{-6}$[35]). The crystalline samples were
mounted on a computer-controlled rotation stage enabling a precise and
reproducible variation of the angle $\theta$ between the main axis of the
crystal and the laser polarization axis.

##  III Experimental Results

Figure 1(b) and (c) show exemplary high-harmonic spectra obtained in a
$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick MgO crystal, cut along the
(100) plane, at two different orientation angles
($\theta=\,$$25\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$, for minimum and maximum signal
strength) and the absorption of the SWIR laser pulse as a function of
$\theta$, for two different peak intensities ($I_{1}=\,$
$11\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, $I_{2}=\,$
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$). The harmonic yield of
all observed orders is higher for
$\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ when compared to
$\theta=\,$$25\text{\,}\mathrm{\SIUnitSymbolDegree}$ (the crystal angles were
calibrated by perturbative third harmonic generation at a moderate excitation
intensity of $\sim$$2\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$). Based
on the inset shown in Fig. 1(b) illustrating the cubic crystal structure of
MgO, such a yield anisotropy can be attributed to high-symmetry directions in
the crystal lattice. At an angle of $45\text{\,}\mathrm{\SIUnitSymbolDegree}$
the real space trajectories of the excited electron wavepackets point towards
the nearest-neighbour of the same species (Mg-Mg or O-O) inside the crystal
structure (monoatomic nearest-neighbour direction).
The SWIR absorption exhibits no pronounced $\theta$-dependence at a low peak
intensities below $\sim\,$$8\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$.
Instead, the absorption stays close to 0 after accounting for the transmission
losses due to Fresnel reflection at the interfaces
($\sim\,$$16\text{\,}\mathrm{\char 37\relax}$ for MgO). At an intermediate
SWIR excitation intensity of
$11\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ the average absorption
increases to $\sim$$4\text{\,}\mathrm{\char 37\relax}$ [yellow circles in Fig.
1(c)]. At the same time distinct extrema begin to emerge. When the excitation
intensity approaches the laser-induced damage threshold of the used samples
(determined to be $\sim$$20\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$),
the modulation depth of these oscillations increases and a clear eight-fold
symmetry becomes apparent in the $\theta-$dependent absorption of the SWIR
pump laser pulse with maxima appearing in the diatomic
($\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$)
and monoatomic
($\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$)
nearest-neighbour directions.

Figure 2: Experimentally determined HHG and absorption in MgO. (a) Angular
dependence of the third harmonic in the perturbative (at $I_{0}=\,$
$6\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, orange circles) and in the
non-perturbative regime (at $I_{0}=\,$
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, grey circles). (b)
$\theta$-dependence of the observed odd harmonics and the SWIR absorption in a
$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick MgO crystal obtained at a
peak intensity of $18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. The
cubic crystal structure of MgO is sketched in the background to indicate the
nearest-neighbour directions. The harmonic signals as well as the absorption
are normalized and offset vertically for clarity. (c) Amplitudes extracted
from the fit function (Eq. 1) for the observed harmonic orders as a function
of the SWIR excitation intensity. (d) Average transmission and modulation
amplitude of the fundamental pump laser pulse as a function of intensity.

With the goal to link the nonlinear absorption and the emission of high-order
harmonics, the orientation-dependent harmonic yields are compared to the
nonlinear absorption in Fig. 2.
Similar to the angular distribution of the absorption, the third harmonic
signal exhibits an eight-fold symmetry at a SWIR peak intensity of
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ [see dark green circles in
Fig. 2(a)]. As reported in Ref. [27] the maxima of the third harmonic emission
along the monoatomic nearest-neighbour directions
($\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$)
are suppressed at low intensity where perturbative mechanisms dominate the
nonlinear response [orange circles in Fig. 2(a)]. Hence, a transition from a
four-fold symmetry to an eight-fold symmetry can be observed for the third
harmonic as the excitation intensity increases.
At an excitation intensity of
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ the harmonic yield of all
observed orders is maximized along the diatomic nearest-neighbour directions
(Mg-O bonding directions at
$\theta=\,$$0.0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$)
as shown in Fig. 2(b). Further maxima are formed along the monoatomic nearest-
neighbour directions (at $\theta=$ $45\text{\,}\mathrm{\SIUnitSymbolDegree}$,
$135\text{\,}\mathrm{\SIUnitSymbolDegree}$,
$225\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$315\text{\,}\mathrm{\SIUnitSymbolDegree}$) resulting in an eight-fold
symmetry of all observed harmonic orders whose extrema align exactly with
those of the nonlinear absorption [pink circles in Fig. 2(b)] suggesting that
the harmonic emission is directly correlated with nonlinear photoionization.
Due to the lower signal strength the fifth and seventh harmonic could only be
observed at intensities where all orders exhibited an eight-fold symmetry.
For further analysis the angle-dependent harmonic yields $Y(\theta)$ as well
as the $\theta$-dependent absorption are approximated by a periodic fit
function according to

$Y(\theta)=A_{0}+A_{\rm{dia}}*\cos[2\theta]^{2n}+A_{\rm{mono}}*\sin[2\theta]^{2n}.$
(1)

Here $A_{0}$ is an offset amplitude, $A_{\rm{dia}}$ and $A_{\rm{mono}}$ denote
the amplitudes associated with the oscillations along the mono- and diatomic
nearest-neighbour directions, respectively. The parameter $n$ determines the
width of the maxima and is well approximated by 3 throughout this work.
In Fig. 2(c) the amplitudes that were introduced in Eq. 1 are analyzed as a
function of the SWIR intensity. The eight-fold angular dependence of the fifth
and the seventh harmonic contains equal contributions of mono- and diatomic
amplitudes. However, for the third harmonic $A_{\rm{mono}}$ and $A_{\rm{dia}}$
evolve differently. While $A_{\rm{dia}}$ constantly exhibits values between
0.3 and 0.5 over the full range of intensities, $A_{\rm{mono}}$ sharply
increases up to an SWIR intensity of
$\sim\,$$4.5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ before following
the same trend as $A_{\rm{dia}}$ at a slightly ($\sim 0.15$) lower level. This
intensity-dependent behavior is consistent with the aforementioned transition
from a four-fold to an eight-fold symmetry as shown in Fig. 2(a).
Figure 2(d) demonstrates how the transmission of the SWIR pulse decreases as
the excitation intensity increases towards the laser-induced damage threshold
of the MgO crystal. At the same time, the modulation amplitude (for the
absorption $A_{\rm{mono}}\approx A_{\rm{dia}}$) of the eight-fold symmetry
increases with intensity in a manner that is notably different from the
intensity dependence of the amplitudes that were discussed in Fig. 2(c). This
is consistent with the fact, that although the transmission is related to the
frequency conversion process, these processes only account for a negligible
part of the transmission loss.

Figure 3: Experimentally determined HHG and absorption in Al2O3. (a)
$\theta$-dependence of HHG and SWIR absorption in Al2O3 at an excitation
intensity of $18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. (b) Fitted
amplitudes $A_{\rm{mono}}$ and $A_{\rm{dia}}$ as a function of the SWIR
intensity. (c) Average transmission and modulation amplitude of the 8-fold
symmetric transmission loss. (d) Sketch of the projection of the hexagonal
unit cell of an M-cut Al2O3 sample. (e) Ratio of the mono- and diatomic
amplitudes extracted from the fit function (Eq. 1) showing the transition from
a stronger $A_{\rm{dia}}$ to a stronger $A_{\rm{mono}}$ with increasing
intensity, finally settling at $A_{\rm{dia}}\sim A_{\rm{mono}}$.

We repeated the same experiments in bulk sapphire samples
($150\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick) cut along the M-plane
($10\overline{1}0$) for which the projection of the hexagonal unit cell onto
the M-plane exhibits a similar structure as MgO [see Fig. 3(d)]. For an
excitation intensity of $18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$
(damage threshold at $\sim\,$
$21\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$) we observe an eight-fold
symmetry for all detected harmonic orders as well as for the nonlinear
absorption of the fundamental SWIR laser pulse [see Fig. 3(a)]. Again, we
observe maxima in the angular distribution of the SWIR absorption when the
high-harmonic emission is maximized. Even though the orientation angles of the
HHG and absorption maxima do not correspond to mono- and diatomic nearest-
neighbour directions we keep the terminology that was introduced in the MgO
case. Figure 3(b) shows the modulation amplitudes extracted from Eq. 1 as a
function of the SWIR peak intensity. Even at the lowest intensity of
$\sim\,$$4\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ the third harmonic
signal exhibits an eight-fold symmetry resulting in comparable magnitudes of
$A_{\rm{mono}}$ and $A_{\rm{dia}}$ [see blue circles and blue triangles in
Fig. 3(b)]. For all observed harmonic orders $A_{\rm{mono}}$ and
$A_{\rm{dia}}$ display a similar intensity dependence. We find a transition
from a stronger diatomic amplitude for intensities below
$12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ to a stronger monoatomic
amplitude, before, at the highest measured intensities, the two amplitudes
become equal again. This observation is confirmed by plotting the ratio
$A_{\rm{dia}}/A_{\rm{mono}}$ in Fig. 3(e). The ratio changes from values
$\geq\,$1 for excitation intensities up to
$12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ to ratios $\leq\,$1 for
intensities between $12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ and
$18\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. This indicates that the
emission from real-space trajectories towards monoatomic nearest-neighbours at
$\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$
becomes more important than the emission from shorter diatomic nearest-
neighbour trajectories at
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$
for $I_{0}\geq\;$$12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. The
transmission results obtained from the sapphire sample display the same
characteristics as in the MgO case with the exception that $A_{\rm{mono}}$ is
constantly slightly larger than $A_{\rm{dia}}$ [see Fig. 3(c)].
In order to verify the generality of our experimental approach we furthermore
analyzed the angular dependence of the HHG process and the nonlinear
absorption in wide-bandgap fluoride crystals. As a prominent example we used
LiF crystals ($200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thickness) with a
cubic crystal structure similar to the one of MgO [see inset in Fig. 4(a)].
Apart from the fact that no transition from a four-fold to an eight-fold
symmetry was observed for the third harmonic, the results shown in Fig. 4
exhibit comparable characteristic features as in MgO. Both, the observed odd
harmonic orders as well as the nonlinear absorption display an eight-fold
symmetry in which the maxima of the HHG signal coincide with those of the SWIR
absorption [Fig. 4(a)]. The angular contrast of the odd harmonics and of the
absorption reaches comparable values as in the MgO case [see Fig. 4(b)] and
the average transmission shows the characteristic drop from the Fresnel-
related losses at low intensity to lower values at higher intensity towards
the damage threshold [see Fig. 4(c)].

Figure 4: Experimentally determined HHG and absorption in LiF. (a)
$\theta$-dependence of HHG and SWIR absorption in LiF. (b) Fitted amplitudes
$A_{\rm{mono}}$ and $A_{\rm{dia}}$ as a function of the SWIR intensity. (c)
Average transmission and modulation amplitude of the 8-fold symmetric
transmission loss.

The experimental results presented in this section provide consistent evidence
for a correlation between the emission of high-order harmonics and nonlinear
absorption of the SWIR pump laser pulse. In all investigated materials the
harmonic emission was the strongest at orientation angles where also the
nonlinear absorption was at its peak. We interpret these findings as an
indication that ionization, which results from the nonlinear absorption,
drives the maximization of the HHG process. In MgO a transition from a four-
fold to an eight-fold symmetry was observed for the angular dependence of the
third harmonic as the SWIR intensity increased. In Al2O3 all detected harmonic
orders exhibited an eight-fold symmetry irrespective of the pump laser
intensity. However, a change of the strongest emission angles from
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$
(called diatomic nearest-neighbour directions in the MgO case) to
$\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$
(monoatomic nearest-neighbour directions) was identified.

Based on our experimental results we identified an unambiguous link between
HHG and laser-induced nonlinear ionization. However, all widely discussed
mechanisms for HHG in solids (interband recollision, intraband dynamics and
ionization-related mechanisms) require the presence of quasi-free, excited
electrons in the conduction bands. Hence, an unambiguous identification of the
dominating HHG mechanism solely based on the presented experimental results is
out of reach.

##  IV Numerical Results and Discussion

To improve our understanding of the dominant HHG mechanism and of the observed
$\theta$-dependent features we performed numerical simulations based on the
SBEs [36, 37] for the case of MgO. In our simulations we implement the band
structure of MgO (taken from Ref. [27], using one valence and one conduction
band) to estimate the $\theta$-dependence of HHG by projecting the electric
field of the SWIR pump laser pulse onto the $\Gamma-K$ and $\Gamma-X$
directions. We then calculate the resulting high-harmonic emission due to the
interband polarization $P(\omega)$ and the intraband current $J(\omega)$
(details of the numerical model can be found in Ref. [38]). The relative
contribution of the two competing mechanisms is displayed in Fig. 5(a).
Strikingly, the intraband current (orange triangles) dominates the emission of
H5, H7 and H9 for SWIR peak intensities
$\geq\,$$1\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ while for the third
harmonic the emission due to the interband polarization (dark circles) remains
stronger. Thus, contrary to previous publications [26, 6, 30, 39, 15] we do
not observe a dominant intraband mechanism for all observed below-bandgap
harmonics. Instead our simulations indicate a transition from a prevailing
intraband mechanism at $0.5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ to
a dominating interband mechanism at a peak intensity of
$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, while both contributions
lead to comparable harmonic yields in between.

Figure 5: Numerical simulations of HHG in MgO using SBEs. (a) Relative
importance of the interband polarization and the intraband current for the
total HHG signal obtained by integrating the harmonic yields for
$\theta\in[0,2\pi]$. (b) - (e) Orientation-dependence of the odd harmonics
(H3-H9) for various SWIR peak intensities.

Figures 5(b)-(e) show the orientation-dependence of H3-H9 for four different
SWIR peak intensities
($0.5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}1\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$).
The angular distribution of H3 changes from a four-fold symmetry with maxima
at
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$
(diatomic nearest-neighbour directions) at
$I_{0}=\,$$0.5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ to an eight-
fold symmetry at $I_{0}=\,$$1\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$.
At even higher SWIR peak intensities a four-fold symmetry with maxima at
$\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$
(monoatomic nearest-neighbour directions) develops. Comparing these
observations to Fig. 5(a) indicates that the $\theta-$dependence associated
with a dominating intraband mechanism (at
$I_{0}=\,$$0.5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, four-fold
symmetry with maxima along diatomic nearest-neighbour directions) differs from
the angular dependence due to the interband polarization (at
$I_{0}=\,$$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, four-fold
symmetry with maxima along monoatomic nearest-neighbour directions). Thus, the
transition from a four-fold to an eight-fold symmetry corresponds to a
transition from a purely intraband case to an intermediate situation where the
combination of the interband polarization and the intraband current generates
an eight-fold symmetry of the total H3 signal. The experimentally observed
transition from a four-fold to an eight-fold symmetry of H3 (compare Fig. 2)
could hence be interpreted as a result of such a combined emission.
The $\theta$-dependence of H5 exhibits an eight-fold symmetry over the full
range of SWIR intensities with one exception at
$I_{0}=\,$$5\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ where the maxima
along the diatomic nearest-neighbour directions break up into two peaks
leading to a twelve-fold symmetry that was experimentally not observed.
However, as the excitation intensity further increases up to
$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, an eight-fold symmetry in
good agreement with our measurements is observed again. In the intensity range
where we could experimentally access the seventh harmonic, the numerical
simulations based on the SBEs reproduce the constant eight-fold symmetry of H7
while H9 tends to display a four-fold symmetric angular dependence whose
preferred orientation for efficient HHG varies between the diatomic and the
monoatomic nearest-neighbour directions [see Fig. 5(e)].
Generally, several microscopic physical mechanisms contribute to the two
observables [$P(\omega)$ and $J(\omega)$] in the SBE-calculations. While
interband recollisions are expected to dominate the above-bandgap response of
the interband polarization $P(\omega)$, high-harmonic generation due to the
intraband current is usually associated with quasi-free carriers being
accelerated to non-parabolic regions of the energy bands. However, other HHG
mechanisms that cannot easily be isolated are naturally included in these two
observables. In particular, ionization-related HHG mechanisms are expected to
severely contaminate the low-order harmonics [16, 31]. We expect the Brunel
mechanism that is traditionally associated with harmonic generation due to
photoionization in gases, clusters and solids to contribute to the intraband
current as it relies upon carrier motion within the excited state. The
injection current, in contrast, is directly connected to the spatial
displacement of an electron during the interband excitation itself and is thus
expected to contribute to the interband polarization. To investigate a
potential contribution of these ionization-related HHG mechanisms that cannot
directly be extracted from the SBE calculations we performed further numerical
simulations of the angular dependence of the photoionization rate. We have
conceived a simple model for the photoionization rate that can support
physical insights and provide realistic qualitative estimates of the angular
dependence of solid-state HHG. Conventional computations of the
$\theta$-dependence of the photoionization rate only consider carrier dynamics
within distinct energy bands, i.e. transitions between Bloch-states that are
based on single-atom wavefunctions localized at the same atom. Here we work in
the Wannier-basis [40, 41, 42] which allows to explicitly resolve transitions
from one atomic site to another. Since such nearest-neighbour transitions are
not included in the dipole momenta used for the SBE calculations we treat them
separately and qualitatively [see Fig. 6(a)]. The photoionization rate (i.e.
the rate of transitions from the valence to the conduction band) is
approximated by the multiphoton formalism [43] as

$\Gamma\approx\sum_{s}|\vec{E}\cdot\vec{s}|^{2N}d_{s}^{2N}$ (2)

where the sum is taken over the neighbouring atomic sites,
$\vec{E}\cdot\vec{s}$ is the projection of the electric field on the direction
towards the atomic site, and $d_{s}$ is the transition dipole moment between a
state in the valence zone and a state in the conduction zone at atomic site
$s$. The electric field modifies the energy difference beetween these states
and correspondingly the number of photons $N$ that is needed for a multiphoton
transition.

Figure 6: Numerical simulations of HHG in MgO using maximally-localized
Wannier-functions. (a) Depiction of electronic excitation at a single atom (as
captured by the SBEs), to interatomic space and to neighbouring atoms
(Wannier-Jumping). (b) Orientation-dependence of the photoionization rate
calculated according to Eq. 2 for different SWIR laser intensities. (c) % (d)
Angular dependence of the photoionization rate of the various nearest-
neighbour and next-nearest neighbour transitions at a SWIR peak intensity of
$0.1\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ (c) and
$20\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ (d).

For MgO we calculated the transition dipole moment based on the known energy
structure of the isolated atoms, the bandgap value and the interatomic
distances. For the MgO crystal the valence-zone electrons are localized at
oxygen atoms, while the localized electrons in the conduction zone are
predominantly positioned at oxygen atoms and partially at magnesium atoms
[44]. For the calculations of the angle-dependent photoionization rate we have
taken nearest-neighbour O-O transitions (corresponding to
$\theta=\,$$45\text{\,}\mathrm{\SIUnitSymbolDegree}135\text{\,}\mathrm{\SIUnitSymbolDegree}225\text{\,}\mathrm{\SIUnitSymbolDegree}315\text{\,}\mathrm{\SIUnitSymbolDegree}$),
second-nearest-neighbour O-O transitions (corresponding to
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$)
and nearest-neighbour O-Mg transitions (corresponding to
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$)
into account [see Fig. 6(b)-(d)]. We only analyze the $\theta$-dependence of
the photoionization rate itself since all ionization-related harmonic orders
are proportional to this rate and will inherit the same angular distribution.
Numerical results for the $\theta$-dependence of the ionization rate are
presented in Fig. 6(b). At a low SWIR intensity of
$2\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ a four-fold symmetry with
maxima along the monoatomic nearest-neighbour directions is observed. This
corresponds to the same directions for which the interband polarization
dominates H3 [see brown line for
$I_{0}=\,$$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ in Fig. 5(b)].
As the excitation intensity increases, a transition from a four-fold to an
eight-fold symmetry is found in the $\theta$-dependence of the photoionization
rate [see orange line for a peak intensity of
$13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ in Fig. 6(b)]. The
additional maxima appearing at
$\theta=\,$$0\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}180\text{\,}\mathrm{\SIUnitSymbolDegree}270\text{\,}\mathrm{\SIUnitSymbolDegree}$
can be associated with ionization along the diatomic nearest-neighbour
directions (excitation from an O ion to a Mg ion or vice versa). Since this is
a direct signature of Wannier-transitions from states localized at one ion to
states localized at another ion, the SBE results do not include these maxima
in the interband polarization [see Fig. 6(a)]. The combination of both
numerical approaches indicates that the interband polarization becomes
dominant for H3 at sufficiently high intensities
($\geq\,$$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$), while the
experimentally observed eight-fold symmetry is linked to Wannier-transitions
between adjacent atoms. H3 (photon energy: $2.48\text{\,}\mathrm{eV}$) is
situated well below the bandgap of bulk MgO ($\sim$$7.8\text{\,}\mathrm{eV}$),
hence interband recollision can not be responsible for the detected emission.
As the injection current is the only mechanism that generates below-bandgap
harmonics in a non-perturbative manner and contributes to the interbant
polarization, we attribute the emission of H3 at SWIR intensities
$\geq\,$$10\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ to the injection
current.

##  V Conclusion

In summary, we have presented experimental and numerical results of the
angular dependence of HHG in periodic crystals. To correlate the nonlinear
frequency conversion process to laser-induced ionization we performed
simultaneous transmission measurements that provide information on the
ionization yield at different crystallographic orientations. Our
investigations were focused on bulk MgO crystals with a cubic crystal
structure and were complemented by experimental investigations in Al2O3 and
LiF. In all materials a distinct correlation between HHG and nonlinear
ionization was observed which manifested itself in a well-defined orientation
dependence of both signals. In detail, we observed maxima of the harmonic
emission coinciding with maxima of the nonlinear absorption. We interpreted
this as an enhanced HHG conversion efficiency at angles where laser-induced
ionization is also maximized. A symmetry analysis of the angular dependence of
the experimentally observed signals connected the angles of maximum HHG
emission and SWIR absorption with monoatomic and diatomic nearest-neighbour
directions in the crystal lattice.
To further substantiate our interpretation and investigate the dominant
mechanisms responsible for HHG we performed two sets of numerical simulations
using MgO as a model system. First, we numerically solved the SBEs using one
valence and one conduction band. Our results reproduced the experimentally
observed transition from a four-fold to an eight-fold symmetry of H3 at
intermediate intensity while at high intensity the interband polarization was
found to dominate, resulting in a four-fold symmetry with maxima along the
monoatomic nearest-neighbour directions. For H5 and H7 an eight-fold symmetry
was predicted for the experimentally used SWIR intensities in agreement with
our experimental findings. Within the SBE calculations the transition in the
angular distribution was caused by a change from a dominating intraband
current at low intensities to a prevailing interband polarization at higher
intensities. The observed shift in dominant amplitude from $A_{\rm{dia}}$ to
$A_{\rm{mono}}$ in Al2O3 may also be attributed to the same phenomenon.
Second, we numerically investigated the possible influence of ionization-
related HHG mechanisms on the experimentally detected harmonic spectra. Our
results unveiled a transition from a four-fold to an eight-fold symmetry of
the photoionization rate for increasing intensity where the strongest signal
is predicted along the monoatomic nearest neighbour directions equivalent to
the dominant interband polarization contribution predicted by the SBEs. As
Wannier-jumping (i.e. the transition from one atom to a neighbouring atom) is
not included in the dipole momenta used in the SBE calculations the
additionally emerging peaks leading to the eight-fold symmetry cannot be seen
in the SBE results. Since the injection current is directly associated with
interband excitations and thus contributes to $P(\omega)$ and since the
recollision mechanism can be excluded due to the below-bandgap photon energy,
we attribute - like in our previous work [25, 19] \- the emission of H3 in the
strong-field regime to the injection mechanism, the emission of harmonics due
to the spatial displacement of electrons during the excitation process.

## Funding.

Funding by the German Research Foundation - SFB1477 ”Light-Matter Interaction
at Interfaces,” Project No. 441234705, is gratefully acknowledged. The work of
S.D.C.R.A. and P.M.K. has been carried out at the Advanced Research Center for
Nanolithography (ARCNL), a public-private partnership of the University of
Amsterdam (UvA), the Vrije Universiteit Amsterdam (VU), the Dutch Research
Council (NWO), and the semiconductor equipment manufacturer ASML, and was
partly financed by Toeslag voor Topconsortia voor Kennis en Innovatie (TKI)
from the Dutch Ministry of Economic Affairs and Climate Policy. P.M.K.
acknowledges support from ERC Starting Grant ANACONDA (grant no. 101041819).

## Disclosures.

The authors declare no conflicts of interest.

## Acknowledgements

It is our pleasant duty to thank M. Ivanov and T. Fennel for discussion of our
experimental results and the numerical methods.

## Data availablity.

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

  * Corkum [1993] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Physical Review Letters 71, 1994 (1993).
  * Lewenstein  _et al._ [1994] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Theory of high-harmonic generation by low-frequency laser fields, Physical Review A 49, 2117 (1994).
  * Krausz and Ivanov [2009] F. Krausz and M. Ivanov, Attosecond physics, Reviews of Modern Physics 81, 163 (2009).
  * Sansone  _et al._ [2010] G. Sansone, F. Kelkensberg, J. Pérez-Torres, F. Morales, M. F. Kling, W. Siu, O. Ghafur, P. Johnsson, M. Swoboda, E. Benedetti, _et al._ , Electron localization following attosecond molecular photoionization, Nature 465, 763 (2010).
  * Ghimire  _et al._ [2011] 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, Nature Physics 7, 138 (2011).
  * Schubert  _et al._ [2014] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, _et al._ , Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations, Nature Photonics 8, 119 (2014).
  * Garg  _et al._ [2018] M. Garg, H.-Y. Kim, and E. Goulielmakis, Ultimate waveform reproducibility of extreme-ultraviolet pulses by high-harmonic generation in quartz, Nature Photonics 12, 291 (2018).
  * Garg  _et al._ [2016] M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Klostermann, A. Guggenmos, and E. Goulielmakis, Multi-petahertz electronic metrology, Nature 538, 359 (2016).
  * Ossiander  _et al._ [2022] M. Ossiander, K. Golyari, K. Scharl, L. Lehnert, F. Siegrist, J. Bürger, D. Zimin, J. Gessner, M. Weidman, I. Floss, _et al._ , The speed limit of optoelectronics, Nature Communications 13, 1 (2022).
  * Stolen and Ashkin [1973] R. Stolen and A. Ashkin, Optical Kerr effect in glass waveguide, Applied Physics Letters 22, 294 (1973).
  * Franken  _et al._ [1961] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Generation of optical harmonics, Physical Review Letters 7, 118 (1961).
  * Elu  _et al._ [2017] U. Elu, M. Baudisch, H. Pires, F. Tani, M. H. Frosz, F. Köttig, A. Ermolov, P. S. J. Russell, and J. Biegert, High average power and single-cycle pulses from a mid-IR optical parametric chirped pulse amplifier, Optica 4, 1024 (2017).
  * Luu  _et al._ [2015] 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 (2015).
  * Krause  _et al._ [1992] J. L. Krause, K. J. Schafer, and K. C. Kulander, High-order harmonic generation from atoms and ions in the high intensity regime, Physical Review Letters 68, 3535 (1992).
  * Goulielmakis and Brabec [2022] E. Goulielmakis and T. Brabec, High harmonic generation in condensed matter, Nature Photonics 16, 411 (2022).
  * Brunel [1990] F. Brunel, Harmonic generation due to plasma effects in a gas undergoing multiphoton ionization in the high-intensity limit, JOSA B 7, 521 (1990).
  * Bauer  _et al._ [1998] D. Bauer, R. Salomaa, and P. Mulser, Generation of ultrashort light pulses by a rapidly ionizing thin foil, Physical Review E 58, 2436 (1998).
  * Geissler, Michael and Tempea, Gabriel and Scrinzi, Armin and Schnürer, Matthias and Krausz, Ferenc and Brabec, Thomas [1999] Geissler, Michael and Tempea, Gabriel and Scrinzi, Armin and Schnürer, Matthias and Krausz, Ferenc and Brabec, Thomas, Light propagation in field-ionizing media: extreme nonlinear optics, Physical Review Letters 83, 2930 (1999).
  * Jürgens  _et al._ [2022] P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, T. Witting, A. Husakou, A. Rouzee, M. Ivanov, T. Fennel, M. J. Vrakking, _et al._ , Characterization of laser-induced ionization dynamics in solid dielectrics, ACS Photonics 9, 233 (2022).
  * Verhoef  _et al._ [2010] A. J. Verhoef, A. V. Mitrofanov, E. E. Serebryannikov, D. V. Kartashov, A. M. Zheltikov, and A. Baltuška, Optical detection of tunneling ionization, Physical Review Letters 104, 163904 (2010).
  * Babushkin  _et al._ [2022] I. Babushkin, Á. J. Galán, J. R. C. de Andrade, A. Husakou, F. Morales, M. Kretschmar, T. Nagy, V. Vaičaitis, L. Shi, D. Zuber, _et al._ , All-optical attoclock for imaging tunnelling wavepackets, Nature Physics 18, 417 (2022).
  * Gao  _et al._ [2019] X. Gao, B. Shim, and M. C. Downer, Brunel harmonics generated from ionizing clusters by few-cycle laser pulses, Optics Letters 44, 779 (2019).
  * Mitrofanov  _et al._ [2011] A. V. Mitrofanov, A. J. Verhoef, E. E. Serebryannikov, J. Lumeau, L. Glebov, A. M. Zheltikov, and A. Baltuška, Optical detection of attosecond ionization induced by a few-cycle laser field in a transparent dielectric material, Physical Review Letters 106, 147401 (2011).
  * Li  _et al._ [2022] W. Li, A. Saleh, M. Sharma, M. Sierka, C. Hünecke, M. Neuhaus, L. Hedewig, B. Bergues, M. Alharbi, A. M. Azzeer, _et al._ , Resonance effects in Brunel harmonic generation in thin film organic semiconductors, arXiv preprint arXiv:2211.07062 (2022).
  * Jürgens  _et al._ [2020] P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, D. Engel, A. Husakou, T. Witting, M. Ivanov, M. Vrakking, T. Fennel, _et al._ , Origin of strong-field-induced low-order harmonic generation in amorphous quartz, Nature Physics 16, 1035 (2020).
  * Lanin  _et al._ [2019] A. A. Lanin, E. A. Stepanov, A. V. Mitrofanov, D. A. Sidorov-Biryukov, A. B. Fedotov, and A. M. Zheltikov, High-order harmonic analysis of anisotropic petahertz photocurrents in solids, Optics Letters 44, 1888 (2019).
  * You  _et al._ [2017] Y. S. You, D. A. Reis, and S. Ghimire, Anisotropic high-harmonic generation in bulk crystals, Nature Physics 13, 345 (2017).
  * Uzan  _et al._ [2020] A. J. Uzan, G. Orenstein, Á. Jiménez-Galán, C. McDonald, R. E. Silva, B. D. Bruner, N. D. Klimkin, V. Blanchet, T. Arusi-Parpar, M. Krüger, _et al._ , Attosecond spectral singularities in solid-state high-harmonic generation, Nature Photonics 14, 183 (2020).
  * Suthar  _et al._ [2022] P. Suthar, F. Trojánek, P. Malỳ, T. J.-Y. Derrien, and M. Kozák, Role of Van Hove singularities and effective mass anisotropy in polarization-resolved high harmonic spectroscopy of silicon, Communications Physics 5, 1 (2022).
  * Vampa  _et al._ [2015a] G. Vampa, T. Hammond, N. Thiré, B. Schmidt, F. Légaré, C. McDonald, T. Brabec, and P. Corkum, Linking high harmonics from gases and solids, Nature 522, 462 (2015a).
  * Vampa and Brabec [2017] G. Vampa and T. Brabec, Merge of high harmonic generation from gases and solids and its implications for attosecond science, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 083001 (2017).
  * Gertsvolf  _et al._ [2008] M. Gertsvolf, H. Jean-Ruel, P. Rajeev, D. Klug, D. Rayner, and P. Corkum, Orientation-dependent multiphoton ionization in wide band gap crystals, Physical Review Letters 101, 243001 (2008).
  * Zhang  _et al._ [2021] Y. Zhang, T. Huang, L. Li, P. Lan, and P. Lu, Intensity and wavelength dependence of anisotropic nonlinear absorption inside MgO, Optical and Quantum Electronics 53, 1 (2021).
  * Heiner  _et al._ [2018] Z. Heiner, V. Petrov, G. Steinmeyer, M. J. Vrakking, and M. Mero, 100-kHz, dual-beam OPA delivering high-quality, 5-cycle angular-dispersion-compensated mid-infrared idler pulses at 3.1 $\mu$m, Optics Express 26, 25793 (2018).
  * Gholam-Mirzaei  _et al._ [2017] S. Gholam-Mirzaei, J. Beetar, and M. Chini, High harmonic generation in ZnO with a high-power mid-IR OPA, Applied Physics Letters 110, 061101 (2017).
  * Lindberg and Koch [1988] M. Lindberg and S. W. Koch, Effective Bloch equations for semiconductors, Physical Review B 38, 3342 (1988).
  * Yu  _et al._ [2016] C. Yu, X. Zhang, S. Jiang, X. Cao, G. Yuan, T. Wu, L. Bai, and R. Lu, Dependence of high-order-harmonic generation on dipole moment in SiO2 crystals, Physical Review A 94, 013846 (2016).
  * Abbing  _et al._ [2022] S. D. Abbing, F. Campi, B. de Keijzer, C. Morice, Z.-Y. Zhang, M. L. van der Geest, and P. M. Kraus, Efficient extreme-ultraviolet high-order wave mixing from laser-dressed silica, arXiv preprint arXiv:2209.15561 (2022).
  * Vampa  _et al._ [2015b] G. Vampa, T. Hammond, N. Thiré, B. Schmidt, F. Légaré, C. McDonald, T. Brabec, D. Klug, and P. Corkum, All-optical reconstruction of crystal band structure, Physical review letters 115, 193603 (2015b).
  * Osika  _et al._ [2017] E. N. Osika, A. Chacón, L. Ortmann, N. Suárez, J. A. Pérez-Hernández, B. Szafran, M. F. Ciappina, F. Sols, A. S. Landsman, and M. Lewenstein, Wannier-Bloch approach to localization in high-harmonics generation in solids, Physical Review X 7, 021017 (2017).
  * Parks  _et al._ [2020] 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 (2020).
  * Brown  _et al._ [2022] G. G. Brown, Á. Jiménez-Galán, R. E. Silva, and M. Ivanov, A Real-Space Perspective on Dephasing in Solid-State High Harmonic Generation, arXiv preprint arXiv:2210.16889 (2022).
  * Ambrosek  _et al._ [2003] D. Ambrosek, M. Oppel, L. González, and V. May, Theory of ultrafast non-resonant multi-photon transitions: basics and application to CpMn (CO) 3, Chemical Physics Letters 380, 536 (2003).
  * De Boer and De Groot [1998] P. De Boer and R. De Groot, The conduction bands of MgO, MgS and HfO2, Journal of Physics: Condensed Matter 10, 10241 (1998).