††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>aachen.de††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>frankfurt.de

# Higher-harmonic generation in boron-doped silicon from band carriers and
bound-dopant photoionization

Fanqi Meng  Frederik Walla  Physikalisches Institut, J. W. Goethe-Universität,
Max-von-Laue-Strasse 1, 60438 Frankfurt am Main, Germany  Sergey Kovalev  Jan-
Christoph Deinert  Igor Ilyakov  Min Chen  Alexey Ponomaryov  Helmholtz-
Zentrum Dresden-Rossendorf, Bautzner Landstrasse 400, 01328 Dresden, Germany
Sergey G. Pavlov  Institute of Optical Sensor Systems, German Aerospace Center
(DLR), 12489 Berlin, Germany  Heinz-Wilhelm Hübers  Institute of Optical
Sensor Systems, German Aerospace Center (DLR), 12489 Berlin, Germany  Institut
für Physik, Humboldt-Universität zu Berlin, 12489 Berlin, Germany  Nikolay V.
Abrosimov  Leibniz-Institut für Kristallzüchtung (IKZ), 12489 Berlin, Germany
Christoph Jungemann  Institut für Theoretische Elektrotechnik, RWTH Aachen,
52062 Aachen, Germany  Hartmut G. Roskos  Physikalisches Institut, J. W.
Goethe-Universität, Max-von-Laue-Strasse 1, 60438 Frankfurt am Main, Germany
Mark D. Thomson  Physikalisches Institut, J. W. Goethe-Universität, Max-von-
Laue-Strasse 1, 60438 Frankfurt am Main, Germany

###### Abstract

We investigate ultrafast harmonic generation (HG) in Si:B, driven by intense
pump pulses with fields reaching $\sim$100$~{}\mathrm{kV}\mathrm{cm}^{-1}$ and
a carrier frequency of 300 GHz, at 4 K and 300 K, both experimentally and
theoretically. We report several novel findings concerning the nonlinear
charge carrier dynamics in intense sub-THz fields. (i) Harmonics of order up
to $n=9$ are observed at room temperature, while at low temperature we can
resolve harmonics reaching even $n=13$. The susceptibility per charge carrier
at moderate field strength is as high as for charge carriers in graphene,
considered to be one of the materials with the strongest sub-THz nonlinear
response. (ii) For $T=300$ K, where the charge carriers bound to acceptors are
fully thermally ionized into the valence subbands, the susceptibility values
decrease with increasing field strength. Simulations incorporating multi-
valence-band Monte-Carlo and finite-difference-time-domain (FDTD) propagation
show that here, the HG process becomes increasingly dominated by energy-
dependent scattering rates over the contribution from band non-parabolicity,
due to the onset of optical-phonon emission, which ultimately leads to the
saturation at high fields. (iii) At $T=4$ K, where the majority of charges are
bound to acceptors, we observe a drastic rise of the HG yields for internal
pump fields of $\sim$30$~{}\mathrm{kV}\mathrm{cm}^{-1}$, as one reaches the
threshold for tunnel ionization. We disentangle the HG contributions in this
case into contributions from the initial “generational”- and subsequent band-
nonlinearities, and show that scattering seriously degrades any coherent
recollision during the subsequent oscillation of the holes.

##  I Introduction

The future of semiconductor electronic devices relies on continuing to test
and refine the physical description of carrier dynamics for increasingly
shorter time scales and higher electric field strengths (i.e., sub-picosecond
times and approaching the $\mathrm{MV}\,\mathrm{cm}^{-1}$-range). Experiments
with ultrashort, high-field pulses, in particular observing nonlinear higher
order harmonic generation (HG) [1, 2, 3, 4, 5], provide a sensitive probe into
the carrier response, as these nonlinear currents are dictated by the precise
carrier distribution and scattering processes, band structure, interactions
between band carriers and dopant/impurity ions, and, at low-temperatures
(where charges are bound to acceptors), field-driven ionization.

In the journey of HG from free-carriers in semiconductors, pioneering studies
using mid- [6] and far-infrared [7, 8] pulses from molecular lasers were
generally restricted to $\sim$nanosecond pulses. Although the estimated peak
electric field reached over $\sim$100$~{}\mathrm{kV}\mathrm{cm}^{-1}$, no
harmonics higher than third order were reported, and the sub-picosecond
dynamics were not resolved. Nevertheless, they fueled the development of
theoretical descriptions, including the role of band non-parobolicity and
carrier relaxation [9, 10]. The availability of high-field (sub-)picosecond
THz pulses (e.g. from free-electron lasers or femtosecond-amplifier-laser
sources) allows one to finally reach the non-perturbative field regime on the
native time scales for the carrier dynamics [1, 5]. Besides fundamental
interest, the HG process can be considered for practical use for frequency
conversion or frequency comb generation, and the efficiency can be enhanced,
e.g., using resonant cavities [11].

In this publication, we investigate HG in weakly p-doped bulk silicon in a
focused 0.3-THz radiation field reaching maximum incident peak fields of
$\sim$100$~{}\mathrm{kV}\mathrm{cm}^{-1}$ (all field values here and in the
following specify the field strength in vacuum, as the field strength within
the specimens varies locally due to the formation of standing waves). These
fields are strong enough to induce tunnel ionization of the impurity atoms at
low temperature and to drive charge carriers out in $k$-space to sufficiently
high kinetic energies to reach the threshold for scattering by optical phonon
emission (energy of the optical phonons at the $\Gamma$ point:
$\mathcal{E}_{op}=63.3$ meV – the TO and LO phonons being degenerate there
[12]).

Figure 1:  Overview of strong-field harmonic generation in Si:B at (a) $T=300$
K – band motion with thermally ionized holes including both band non-
parabolicity and energy-dependent scattering; (b) $T=4$ K, whereby light holes
are first generated by tunnel ionization of acceptors, followed by significant
lh$\rightarrow$hh scattering (amplitude of hole trajectories not to scale).
(c) Example of experimental emission intensity spectra for both $T=300$ K
(peak external incident pump field $E_{0}=115~{}\mathrm{kV}\mathrm{cm}^{-1}$)
and $T=4$ K ($E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$) (linear transmitted
fundamental spectrum with undoped sample included for comparison).

##  II Experimental

The experiments were conducted at the HZDR (Dresden) by employing the TELBE
superradiant undulator source, where we used linearly polarized THz pulses
tuned to a carrier frequency of $\nu_{1}=$ 300 GHz, with a pulse duration of
$T_{1}=$ 14.1 ps (FWHM) and pulse energies up to 1$~{}\mu\mathrm{J}$ at a
50-kHz repetition rate. We performed two series of experiments with two Si:B
samples: (i) at $T=300$ K ($N_{d}=5.75\cdot 10^{15}~{}\mathrm{cm}^{-3}$,
thickness $L=900$ µm) where the dopant acceptors are fully ionized and one
measures the nonlinear response of a constant population of band holes, and
(ii) at $T=4$ K ($N_{d}=5.0\cdot 10^{16}~{}\mathrm{cm}^{-3}$, thickness
$L=272$ µm) where essentially all holes are bound to their parent ions and are
photoionized into the valence bands during the pump pulse. Preceding the
sample, we employed a low-pass filter to suppress residual harmonics from the
source, followed by two polarizers [2] to allow continuous variation of the
strength of the pump field in the sample. Following the sample, a calibrated
high-pass filter [1] was used to reduce the amplitude of the fundamental to
allow a more balanced signal level for the overtones. The emitted field from
the sample was reimaged into an electro-optic (EO) crystal for coherent time-
domain detection. Incident and detected on-axis temporal field strengths were
determined by calibrating the reference signal with the fluence determined
from an additional measurement of the pump beam profile. Details of the
samples and determination of doping concentration, field calibration, and
spectral correction for the filters and EO-response are given in the
Supplementary.

##  III Results

###  III.1 Overview of the findings

Figure 1 shows conceptual aspects of the experiments performed at 300 K (Fig.
1(a)) and 4 K (Fig. 1(b)) together with spectra recorded at the two
temperatures at the maximal field strength available in the experiments (Fig.
1(c)). At room temperature, the impurities are ionized and one observes the
signatures of the nonlinear response of the charge carriers accelerating and
scattering in the incident THz field (linearly polarized in $x$-direction,
propagating in $z$-direction). As indicated in the lower panel of Fig. 1(a),
the nonlinear response is determined by the non-parabolicity of the heavy- and
light-hole valence bands and by the relaxational nonlinearity associated with
the absorption and emission of acoustic and optical phonons, which gives rise
to HG due to the energy dependence of the scattering rates. Charged impurity
scattering turns out not to play a strong role for the nonlinear response. At
low $T$, most of the impurity atoms are initially in the neutral state. When
they are exposed to a radiation field with rising field amplitude, one
observes an increase of the nonlinear susceptibility $\chi^{(n)}$, which
originates from the growth of the number of free holes by impact ionization
and – setting in at fields of $\sim$30$~{}\mathrm{kV}\mathrm{cm}^{-1}$ – by
additional tunnel ionization, as indicated in Fig. 1(b). The tunneling process
mainly populates the light-hole states from which scattering rapidly
redistributes the holes in the valence band. We show that near the ionization
threshold, a significant HG contribution is due to these “generational”
nonlinear currents, while for higher pump fields the subsequent band
nonlinearities dominate.

###  III.2 Room-temperature harmonic generation

We begin with the detailed description of the results measured at $T=300$ K,
to study the nonlinear response of a constant density of thermally ionized
band holes. Examples of the fields measured behind the sample are shown in
Fig. 2(a) for an incident peak pump field of
$E_{0}=115~{}\mathrm{kV}\mathrm{cm}^{-1}$, with the corresponding intensity
spectrum and spectrogram in (b) and (c), respectively (the spectrogram
calculated with a temporal gate $w(t)=e^{-2t^{2}/T_{g}^{2}}$ with $T_{g}=10$
ps). Here, the odd overtones from $n=3$ to 9 are clearly resolved. As shown in
Fig. 2(b), the relative spectral intensity vs. harmonic $n$ is quite close to
that expected for an ideal nonlinear HG process ($I_{n}\sim a^{-n}$, dashed
curve), as is the progressive reduction in pulse duration $T_{n}$ for each
harmonic in Fig. 2(c), where some deviations are expected due to
saturation/propagation effects discussed below.

The theoretical modelling of the HG emission is based on Monte-Carlo (MC)
simulations of an ensemble of holes in the time-dependent field. The
simulations include the dynamics of the heavy-hole (hh), light-hole (lh) and
split-off (so) valence bands, realistic 3D band structures, and
acoustic-/optical-phonon scattering (as well as charged impurity scattering,
when expected to contribute significantly), similar to the treatment employed
previously in the context of high-frequency Si devices [13, 14]. As
significant propagation effects also occur, including standing-wave effects
for the pump field [1], one must go beyond a description of the local response
to describe the experimental results. Hence we embedded the MC simulations in
a 1D (plane-wave) finite-difference time-domain (FDTD) propagation scheme with
a full, self-consistent approach for the charges/fields for the simulations at
high fields (as opposed to an earlier perturbative approach [15]), typically
with a step size of $\delta z=1~{}\mu\mathrm{m}$ along the propagation
direction. To investigate the microscopic HG mechanisms, local MC simulations
were performed, i.e., at a single internal $z$-point in the bulk, as in our
previous report [1].

An example of the emitted spectrum/spectrogram from such MC-FDTD simulations
is shown in Fig. 2(b) and (d), respectively, and shows very good agreement
with the experimental data. Note that the noise floor of the MC-FDTD results
is due to the statistical fluctuations for the ensemble, which was restricted
to $\leq 10^{7}$ holes in order to achieve practical calculation times. As
discussed below, these results contain both local-saturation and propagation
effects, involving changes in the carrier distribution and scattering dynamics
during the pump pulse.

Figure 2:  (a) Experimental transmitted temporal field for a Si:B sample
($N_{d}=5.75\cdot 10^{15}~{}\mathrm{cm}^{-3}$, $T=300$ K) for an incident pump
peak field of $E_{0}=115~{}\mathrm{kV}\mathrm{cm}^{-1}$ (blue curve), as well
as the field of the overtones without fundamental (red) and scaled reference
data for an undoped sample (black). Inset shows detail around the pulse peak.
(b) Corresponding intensity spectra, as well as spectrum from MC-FDTD
simulations with $E_{0}=100~{}\mathrm{kV}\mathrm{cm}^{-1}$ (see text for
discussion of pump field scaling between theory and experiment). (c)
Experimental and (d) MC-FDTD spectrograms (common color scale as shown).
Included in (b) and (c) are dashed curves for the relative intensity and
duration, respectively, for ideal nonlinear HG (for the latter, assuming a
Gaussian pulse profile, such that $T_{n}\sim T_{1}/\sqrt{n}$, curves shown for
$T_{n}=\pm 2T_{\mathrm{FWHM}}$).

The dependence of the amplitudes of the transmitted harmonic fields on that of
the incident pump field is shown in Fig. 3(a), for both experimental and MC-
FDTD results. The data are obtained via band-pass filtering the spectral
fields about each $\omega_{n}=2\pi\nu_{n}$ and extracting the peak temporal
field after transforming each back to the time-domain (see Supplementary for
temporal waveforms vs. $n$). In order to aid visual comparison, the
experimental pump field $E_{0}$ is scaled down by a factor of 0.9, while for
$n=5,7,9$ modest correction factors are applied to the field values $E_{n}$
(as listed in the caption). A comparison of the unscaled experimental and
simulated results is included in the Supplementary. Given the experimental
error margin for the field calibration, the agreement is seen to be very
reasonable, especially as no scaling is applied to $E_{1,3}$. Also included
are power-law fits $E_{n}\sim E_{0}^{\eta_{n}}$ using the lowest available
field ranges in each case; the exponent $\eta_{n}$ is plotted in Fig. 3(b).
Several key aspects are evident: Firstly, one sees with increasing field that
the transmitted fundamental wave ($n=1$) grows super-linearly, in both the
experimental and MC-FDTD data. As discussed further below, this is due to the
increase of the scattering rate of the holes at higher fields, which causes a
reduction of the Drude absorption at $\nu_{1}$, and amounts to a bleach factor
of ${\sim}2$ at the highest field (the linear Drude absorption depth of the
sample at $\nu_{1}$ is 230$~{}\mu\mathrm{m}$). For $n=3$, one observes an
initial growth close to the ideal power-law exponent of $\eta_{3}=3$, although
the data visibly saturate by a factor ${\sim}2$ at the highest fields, despite
the bleach of the fundamental pump field. A similar saturation is observed for
$n=5,7,9$ (note that the noise levels in the MC-FDTD preclude fitting the data
for $n=9$ at sufficiently low fields to avoid saturation). These non-ideal
field dependencies are due to an interplay between local (microscopic)
saturation and propagation effects, and will be disentangled below with the
aid of local MC simulations. On the basis of the experimental fields, one can
estimate an effective nonlinear susceptibility $\chi^{(n)}$ for each harmonic,
the results of which are shown in Fig. 3(c) for $n=3,5,7$. These are
calculated following standard approaches [2] for uni-directional propagation
(neglecting absorption and phase-matching) which yields
$E_{n}=\beta_{n}\chi^{(n)}E_{1}^{n}$ with
$\beta_{n}=i\omega_{n}L/(2^{n}n_{r}c)A_{n}$ ($n_{r}$ the refractive index, $c$
the vacuum speed of light). Due to the significant Drude absorption for both
the pump and harmonics (which would lead to a underestimate of $\chi^{(n)}$
from the data), we include the factors
$A_{n}=(1-e^{-a_{n}L}/a_{n}L)e^{-\alpha_{n}L/2}$
($a_{n}=(n\alpha_{1}-\alpha_{n})/2$) where $\alpha_{n}$ is the intensity
absorption coefficient for each harmonic (the experimental bleach factor is
also included to correct $\alpha_{1}$ to avoid overestimating $\chi^{(n)}$,
but only leads to a small correction, see Supplementary). Clearly, the
saturation effects with increasing $E_{0}$ lead to a reduction in the
extracted values. For the lowest fields, we have $\chi^{(3)}=1.0\cdot
10^{-13}~{}\mathrm{m}^{2}\,\mathrm{V}^{-2}$, which is highly consistent with
the value $\chi^{(3)}=0.9\cdot 10^{-13}~{}\mathrm{m}^{2}\,\mathrm{V}^{-2}$ one
can deduce from previous reports for Si:B with nanosecond pulses at
$\nu_{1}=610$ GHz [7] for this dopant density.

We note that besides the pump-field dependence of the emission field strength,
we also extracted the relative phases of the fields. As described in the
Supplementary, these vary by a small fraction of the fundamental cycle
($<0.2\pi$) while the trends vs. $E_{0}$ are also reasonably well reproduced
in the MC-FDTD results. A comparison with the local-MC response shows these
trends comprise both contributions from the local intra-cycle phase as well as
a dominant contribution from the accumulated phase shifts due to propagation.

Figure 3:  (a) Pump-field dependence of transmitted harmonic peak fields
($N_{d}=5.75\cdot 10^{15}~{}\mathrm{cm}^{-3}$, $T=300$ K, $L=900$ µm):
Experimental results (open circles), and MC-FDTD results (filled squares) for
each odd harmonic $n=1-9$. Respective power-law fits included as straight
lines. To aid comparison, experimental peak pump fields ($E_{0}$, incident
external to the sample) scaled by a factor 0.9, while higher-harmonic fields
are scaled by 1.7 ($n=5$), 2.0 ($n=7$), and 2.5 ($n=9$). (b) Power-law
exponents from (a) vs. $n$, from experimental and MC-FDTD data, as well as
single-point (bulk) MC simulations. Note that the resolvable range of
experimental data for $n=7,9$ are already well in the saturation regime. (c)
Nonlinear susceptibilities $\chi^{(n)}$ calculated from experimental peak
fields in (a) after correction for Drude absorption (see main text).

In Fig. 4, we plot a selection of results from local (single $z$-point) MC
simulations vs. field, where we stored and analyzed the time-dependent,
ensemble-averaged band velocities $v_{b}(t)$
($b=\mathrm{hh},\mathrm{lh},\mathrm{so}$) and occupation-weighted total
velocity $v(t)$, wavevectors $k_{bj}(t)$ and their rms spreads
$\sigma_{bj}(t)$ ($j=x,y,z$, where the pump field is polarized along the
$x$-direction), band populations $N_{b}(t)$ and scattering rates
$\Gamma_{b}(t)$. As the so-band has only a very minor contribution, it was
omitted in the following analyses. In Fig. 4(a), we plot the relative emission
intensities $V_{n}$ for each harmonic $n$, obtained by integrating the
spectral intensity of $V(\omega)=|\mathcal{F}\\{v(t)\\}|^{2}$ about each
$\omega_{n}$. Also included are the fluences $F_{n}\propto|E_{n}|^{2}$ from
the MC-FDTD results in Fig. 3(a) (scaled for comparison by a common constant
for all $n$), which in the absence of any propagation effects would coincide
with $V_{n}$ for $n>1$. One observes saturation in the local response with
increasing field for all overtones as well as the fundamental. For the latter,
its origin is clearly demonstrated in Fig. 4(b), where we plot $v(t)$ for hh
and lh with a scaled profile of the pump field $E(t)$. The inset displays
three oscillation periods at the center of the pulse. One sees how $v(t)$
undergoes strong clipping at each peak of $E(t)$, which is due to a rapid rise
in momentum scattering, discussed further below. Integrated along the beam
path, the decrease in local emission intensity is, however, overcompensated by
reduced loss as shown by the blue lines (full and dotted) in Fig. 4(a). The
loss reduction is also a consequence of the increase in scattering rate for
the local current which manifests as a weaker Drude absorption at $\omega_{1}$
during propagation, and leads to the bleach effect seen in Fig. 3(a).
Comparing the local and FDTD results in Fig. 4(a) allows one to draw
additional conclusions concerning the role of propagation for the harmonics
$n>1$: The local saturation in the emission intensity for $n=3$ is compensated
by the bleach of the Drude absorption during propagation, weakening the
saturation effect for the transmitted signal. This effect is also present to a
lesser degree for $n=5$, but is essentially absent for $n=7,9$, as here the
increased scattering rate no longer causes a significant bleach of the Drude
absorption. Moreover, with increasing $n$, the effects of a finite phase-
mismatch $\Delta k_{n}\sim\omega_{n}(n_{rn}-n_{r1})/c$ should also lead to
additional saturation effects during propagation. In Fig. 4(c), we compare the
relative contributions to the total emission intensity from hh and lh for each
$n$. As expected, the linear response ($n=1$) is dominated by hh, due to their
much higher occupation, $N_{\mathrm{hh}}/N_{0}\approx 0.85$ in thermal
equilibrium, which does not vary significantly during the pump excitation (in
contrast to our previous simulations at $T=10$ K [1] where $N_{\mathrm{hh}}$
rose to $0.93$ in non-equilibrium, further depleting $N_{\mathrm{lh}}$). For
the overtones, the relative contribution from lh increases with $n$,
indicating that they have a stronger nonlinear response per-hole than the hh
(see below).

Figure 4:  Local MC results ($N_{d}=5.75\cdot 10^{15}~{}\mathrm{cm}^{-3}$,
$T=300$ K). (a) Pump-field dependence of relative emission intensities $V_{n}$
for each odd harmonic $n=1-9$ (solid curves, legend in (d)). Scaled fluence
from MC-FDTD (dotted curves) included for comparison (corresponding to data in
Fig. 3(a)). Black dashed curve corresponds to a linear dependence, for
comparison with the curves for $n=1$. Note that local internal pump fields
used are $E^{\prime}_{0}=\hat{t}_{1}E_{0}$ (where $\hat{t}_{1}=0.46$ is the
Fresnel field transmission entering the sample). (b) Time-domain (ensemble-
average) velocity for $E_{0}=110~{}\mathrm{kV}\mathrm{cm}^{-1}$ for both hh
and lh, as well as scaled pump electric field profile $E(t)$. (c) Respective
fraction of emission intensities from hh/lh for each $n$. (d) Ratio of
nonlinear intensity due to $k_{b}(t)$ only ($K_{bn}$, see text) relative to
emission intensity $V_{bn}$. (e, f) Scattering rates $S(t)$ during the peak
cycle of $E(t)$ for two peak pump fields
($E_{0}=55,110~{}\mathrm{kV}\mathrm{cm}^{-1}$, respectively) plotted as loci
vs. the ensemble averages $k_{bx}(t)$. Also included are the same loci, but
vs. $k_{bx}\cdot(1+|\sigma_{bx}/k_{bx}|)$ to reflect the extent of the high
energy tail of the distributions. Thresholds for opt. phonon emission for each
band (corresponding to an energy of $\mathcal{E}_{op}=63.3$ meV) included as
vertical lines. Orientation of temporal hysteresis indicated by arrow in (f).

As mentioned in Sec. III.1, one can distinguish two main mechanisms for the HG
process: band non-parabolicity (“BNP”) and energy-dependent scattering
(“EDS”). To gauge their relative roles, we employ the same approach we used
previously (Ref. [1], Fig. 4), i.e., to compare the emission intensity
$V_{bn}$ for each band with the spectral intensity calculated from the
occupation-weighted wavevector component $k_{bx}(t)$, i.e., $K_{bn}$ obtained
from integrating $K_{b}(\omega)=|\mathcal{F}\\{k_{bx}(t)\\}|^{2}$ about each
$\omega_{n}$ (in contrast, we determined $K_{b}(\omega)$ in Ref. [1] from a
1D$\langle k\rangle$ calculation for a rigid 1D wave packet instead of Monte
Carlo simulations for an ensemble of holes). While the band non-parabolicity
indeed has a finite influence on the precise scattering processes and hence
$k_{bx}(t)$, the quantity $K_{bn}$ still provides a reasonable measure for the
nonlinear response in the hypothetical absence of BNP. By normalising
$K_{b1}\rightarrow V_{b1}$, for the overtones a ratio of
$K_{bn}/V_{bn}\rightarrow$0 or 1 corresponds to pure BNP or EDS, respectively.
As shown in Fig. 4(d), the $K_{bn}/V_{bn}$ values are all above 0.5 in the
field ranges where harmonics can be extracted, increasing toward unity as the
field amplitude increases. It indicates that, for the highest electric field
used in our experiment, EDS dominates. One could also envisage that, at the
electric field below $50~{}\mathrm{kV}\mathrm{cm}^{-1}$, the BNP also plays an
important role for certain higher order harmonics, i.e $n=7,9$. This is in
stark contrast to our earlier study [1] for $n=3,5$ at $T=10$ K with
$N_{d}\leq 10^{14}~{}\mathrm{cm}^{-3}$ and somewhat lower fields
($E^{\prime}_{0}\lesssim 25~{}\mathrm{kV}\mathrm{cm}^{-1}$, $E_{0}\lesssim
50~{}\mathrm{kV}\mathrm{cm}^{-1}$), where $K_{bn}/V_{bn}\ll 0.01$, and hence
BNP heavily dominated the HG response. The major reason for this is, that the
higher fields in the present study allow to reach the threshold for optical
phonon emission (the TO/LO phonon being degenerate at the $\Gamma$ point with
$\mathcal{E}_{op}=63.3$ meV [12]). This is not immediately obvious in Fig.
4(d) where $K_{bn}/V_{bn}$ for $n=3$ does not fall rapidly as the field
strength decreases to $E_{0}\lesssim 50~{}\mathrm{kV}\mathrm{cm}^{-1}$, the
field range of our previous study [1]. This can be traced to the broadening
effect by the higher value of $T$ (300 K here, 10 K in [1]), whereby the
broadened Fermi-Dirac distribution brings the more energetic carriers already
close to $\mathcal{E}_{op}$ in thermal equilibrium before the pump excitation
(the band-filling at these densities $N_{d}$ is still only a minor effect).
This suppresses any sharp threshold behavior for $K_{bn}/V_{bn}$ vs. $E_{0}$
here, although experiments at somewhat lower $T$ (but sufficiently high to
maintain a thermally ionized hole population) should show a much clearer
transition.

In Fig. 4(e,f) we plot the loci of the instantaneous, ensemble-averaged
scattering rates $S(t)$ vs $k_{bx}(t)$ (for a cycle at the peak of the pump
pulse) for peak pump fields $E_{0}=55$ and
$110~{}\mathrm{kV}\mathrm{cm}^{-1}$, respectively. To better reflect the
magnitude of $k_{bx}$ reached by the high-energy tail of the carrier
distribution, we also plot the data vs.
$k_{bx}(t)\cdot(1+|\sigma_{x}(t)/k_{bx}(t)|)$, i.e., $k_{bx}(t)$ extended by
the instantaneous rms width $\sigma_{x}(t)$ of the distribution. Here one
clearly sees that the $k$-space extent of the carrier distributions is
significantly broader than the ensemble average. At $k_{x}=0$, the hh
distribution (unlike the lh distribution) already touches the threshold for
LO/TO phonon emission, although the hh and lh scattering rates are still
dominated by acoustic phonon scattering. When driven out in $k$-space, the hh
and lh populations experience a strong increase of the scattering rate, much
of which appears to be due to reaching and crossing the threshold for LO/TO
phonon emission. At $E_{0}=110~{}\mathrm{kV}\mathrm{cm}^{-1}$, the increase in
scattering is steeper for lh than hh. Moreover, a comparison of $K_{bn}$ vs.
$E_{0}$ (not shown) also shows a higher ratio
$K_{\mathrm{lh},n}$/$K_{\mathrm{lh},1}$ compared to
$K_{\mathrm{hh},n}$/$K_{\mathrm{hh},1}$, i.e., the lh nonlinear response is
dominantly from the higher EDS, and not a higher BNP. These effects arise due
to the fact that the lh initially absorb more kinetic energy from the pump
field, allowing the highest-energy lh to experience stronger optical phonon
emission before the lh$\rightarrow$hh intra-cycle energy equilibration takes
full effect.

###  III.3 Low-temperature harmonic generation

We turn now to the low-temperature case, where the vast majority of acceptors
are in the neutral state and only a small residual density of holes reside in
the valence bands before the pump pulse excitation. In our previous study of
Si:B at low temperature [1], due to both the lower peak pump fields
($E_{0}\lesssim 50~{}\mathrm{kV}\mathrm{cm}^{-1}$) and higher frequency
($\nu_{1}=1.29$ THz), we did not reach the threshold for significant tunnel
ionization. Nevertheless, we could resolve both 3HG and 5HG in the emitted
fields (measured only as intensity spectra obtained with far-infrared Fourier
transform spectrometry), which arose from the residual band-hole density
(which was seen to grow with increasing field due to carrier multiplication
during the pump pulse, resulting in a moderately field-dependent value of
$\chi^{(3,5)}$). Moreover, as mentioned above, the nonlinearity giving rise to
the HG was found to be dominated by BNP, in particular as most holes did not
reach the threshold energy for optical phonon emission whereby EDS should
become very strong.

Here at low temperature, we can reach peak external incident pump fields of
$E_{0}{\sim}80~{}\mathrm{kV}\mathrm{cm}^{-1}$ (and internal fields
$E^{\prime}_{0}{\sim}40~{}\mathrm{kV}\mathrm{cm}^{-1}$), and more decisively
at a significantly lower frequency of $\nu_{1}=0.3$ THz, such that the
ponderomotive energy $U_{p,b}=q^{2}E_{0}^{2}/(4m_{b}\omega_{1}^{2})$ [16] will
be up to a factor ${\sim}50$ higher ($U_{p,\mathrm{lh}}\sim 1.3$ eV for lh
with $m_{\mathrm{lh}}/m_{e}=0.15$), and the Keldysh parameter (an inverse
measure of the onset of tunnel ionization for a required ionization energy
$I_{p}$) $\gamma_{b}=\sqrt{I_{p}/2U_{p,b}}$ also 7 times higher
($\gamma_{\mathrm{lh}}\sim 0.13$, well into the tunnel regime). Note we
tentatively employ these expressions from the gas phase, where scattering
(dephasing) is neglected and one assumes a ballistic trajectory for the
ionized charges, as addressed below.

As per the results in Fig. 2 for $T=300$ K, we first show a summary of the
experimental emitted fields in Fig. 5 for the highest pump field
$E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$. Compared to the results for $T=300$
K, the bandwidths of the overtones here are larger (see also Fig. 1(c)),
increasingly so for higher $n$, and the intensity ratios between successive
overtones are significantly smaller. The latter is also evident in the time-
domain overtone field (Fig. 5(a), red curve) which shows a more complex
interference of the $n=3$ wave with those for $n>3$. The larger bandwidths
could in principle indicate either successively shorter harmonic pulses
(decreasing faster than the $T_{n}\sim T_{1}/\sqrt{n}$ rate for an ideal
nonlinearity with a Gaussian pulse profile) or a frequency modulation of each
HG field: An inspection of the experimental spectrogram (Fig. 5(c)) shows that
the latter effect is significant, which we attribute to the fact that the hole
populations vary during the pulse due to the tunnel ionization (addressed
further below), such that the wave-mixing is more complex than just a multi-
photon convolution of the components of the pump pulse spectrum. That the
broadening is mostly due to local effects (rather than nonlinear refraction
during propagation) is confirmed by comparing the spectra from FDTD-MC and
local-MC simulations (presented below), the latter also showing such increased
spectral widths.

Figure 5:  (a) Experimental transmitted temporal field for Si:B sample
($N_{d}=5.0\cdot 10^{16}~{}\mathrm{cm}^{-3}$, $T=4$ K) for an incident pump
peak field of $E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$ (blue curve), as well
as the field of the overtones without fundamental (red). Inset shows detail
around the pulse peak. Note that a temporal window was applied to the wings of
the pulse to suppress noise and reflections in the sample. (b) Corresponding
intensity spectrum, as well as that from MC-FDTD simulations including tunnel
ionization for $E_{0}=80~{}\mathrm{kV}\mathrm{cm}^{-1}$ ($N_{d}=5.9\cdot
10^{16}~{}\mathrm{cm}^{-3}$). Both spectra plotted with the same absolute
intensity scale. (c) Experimental and (d) MC-FDTD spectrograms. Additional
annotations per Fig. 2.

As in Sec. III.2, we performed MC simulations of the HG experiments. This
required extending the MC scheme to incorporate the tunnel ionization process
injecting holes into the bands, with a time-dependent density $N_{b}(t)$. To
calculate the ionization rate $\partial_{t}N_{b}(t)$, we employed the
ionization probability rates $w_{i,b}(E)$ established in the literature for a
static electric field [17, 18] and integrate
$\partial_{t}N_{b}=w_{i,b}(t)(N_{d}-N)$ ($N=\Sigma_{b}N_{b}$) during each time
step with the field $E(t)$. The use of a static-field model for the
instantaneous tunnel ionization rate has been shown to hold reasonably well in
the context of ionization of gas atoms/molecules, at least in the strong-field
regime ($\gamma\ll 1$) [19, 20, 21]. As asserted in [18], the tunneling vs.
field should be dominated by ionization into the lh band, which seems
reasonable due to the $U_{p,b}\propto 1/m_{b}$ dependence for the
ponderomotive energy. As a simplified approach, we then take
$w_{i,\mathrm{hh}}\rightarrow 0$ and inject holes only into the lh band, which
then can rapidly scatter into the hh band, as modeled by the MC treatment.
Details of the ionization rate (formula, parameters for Si:B and plot of
$w_{i,\mathrm{lh}}(E))$ are given in the Supplementary.

The corresponding intensity spectrum and spectrogram from MC-FDTD calculations
for conditions close to the experimental ones are shown in Fig. 5(b) and (d),
respectively. The simulated results do exhibit a reasonable qualitative
correspondence to the experimental data, although certain systematic
deviations are present, as discussed in the following. We first compare the
dependence of the HG emission vs. pump field in Fig. 6, which includes both
experimental and MC-FDTD data. As the broadband EO field detection was only
sensitive enough to resolve the harmonics for $E_{0}\gtrsim
40~{}\mathrm{kV}\mathrm{cm}^{-1}$, we augmented the measurements of the 3HG-
emission with a more sensitive low-bandwidth EO sensor for pump fields down to
$E_{0}{\sim}10~{}\mathrm{kV}\mathrm{cm}^{-1}$, as shown. As significant
distortions are seen in the envelopes of the filtered time-domain harmonics
(see Supplementary), we choose here to plot the results in terms of the
emitted fluence ($F_{n}=\epsilon_{0}c\int{dtE^{2}_{n}(t)}$) to avoid any
artifacts. One sees that at low pump fields, the 3HG emission closely follows
a power-law dependence, $F_{n}\propto E_{0}^{2\eta_{3}}$ with $\eta_{3}=3.54$.
This behavior is comparable to that seen in our previous study at $T=10$ K
(with the same B dopant concentration as here) [1], where a value of
$\eta_{3}=4.2$ was determined, which exceeds the value $\eta_{3}=3$ (for an
ideal nonlinear process in a static medium) due to field-driven multiplication
[22] of the residual band carriers (density $N_{r}$) during the pump pulse.
The different value of $\eta_{3}=3.54$ here can be attributed to the different
pump frequency. In [1, 11], we employed the impact ionization model from [22]
in its high-field limit, where the carrier multiplication factor $f=\Delta
N/N_{r}\sim 1-E_{m}^{2}/E^{2}$, with the characteristic field constant $E_{m}$
increasing with $\omega_{1}$. Hence $\partial_{E}f\sim+2E_{m}^{2}/E^{3}$,
which should be larger for the previous experiments with $\nu_{1}=1.29$ THz
and hence produce a larger value of $\eta_{3}$, which is at least
qualitatively consistent with the two results.

Turning now to the higher field range $E_{0}\gtrsim
25~{}\mathrm{kV}\mathrm{cm}^{-1}$, one sees that the experimental 3HG emission
grows rapidly (by more than 4 orders of magnitude in fluence) with a
reasonably sharp onset, accompanied by the higher harmonics $n=5,7,9$ almost
in proportion. Indeed, the MC-FDTD simulations also show a similar pump-field
dependence, which is due to the onset of tunnel ionization and resultant
nonlinear response of the ionized band holes, which then begins to saturate
due to similar effects seen above for $T=300$ K. Considering that no rescaling
is applied to the results, the agreement is quite remarkable.

The importance of employing a full FDTD treatment for modeling the
experimental data is demonstrated in Fig. 6(b,c), where we plot the spatial
distributions of both the peak field $E(z)$ (i.e., the maximum amplitude
during FDTD propagation) and photoionized hole density (after propagation,
$N(z)=\Sigma_{b}N_{b}(z)$) for a value of $E_{0}$ at the onset of tunnel
ionization (Fig. 6(b)) and above (Fig. 6(c)). In Fig. 6(b), one clearly sees
the standing-wave profile in $E(z)$ due to multiple reflections of the field
in the sample, which manifests as three field-enhancement peaks (at each
surface and in the center of the Si:B sample), which exceed the field
$E^{\prime}_{0}$ one would have if only accounting for the Fresnel
transmission coefficient $\hat{t}_{1}$ of the incident field $E_{0}$ – this
will be addressed again below in assessing the local-MC results. (Note that
these effects are not significant for the $T=300$ K case above, due to the
strong, pre-existing Drude absorption from the thermally ionized carriers.)
This standing-wave effect also gives rise to a strong dependence of the
photoionized hole profile $N(z)$, and for the case in (b) where one is still
close to the initial, exponential onset of ionization, the small asymmetry in
$E(z)$ manifests as a significant asymmetry in $N(z)$. For the higher-field
case in Fig. 6(c), one sees that the field profile becomes more complex,
affected by depletion and propagation effects, leading to a reshaping of the
standing-wave pattern in $N(z)$. Moreover, the field enhancement is actually
suppressed and one sees peak fields closer to the nominal value
$E^{\prime}_{0}$. Clearly a neglect of these effects (e.g. using a simplified,
single-pass uni-directional propagation model, or a bulk layer – see below)
would have a significant impact on the predicted HG emission.

Figure 6:  (a) Pump-field dependence of transmitted harmonic fluence
($N_{d}=5.0\cdot 10^{16}~{}\mathrm{cm}^{-3}$, $T=4$ K, $L=272$ µm):
Experimental (open circles, with additional higher sensitivity measurements
for $n=3$ as open squares), and MC-FDTD results (filled squares, using
$N_{d}=5.9\cdot 10^{16}~{}\mathrm{cm}^{-3}$) for each odd harmonic $n=1-9$.
Power-law fits for $n=1$ and $n=3$ included as straight dashed lines
(exponents $\eta$ as indicated). (b,c) Spatial profile of peak electric field
$E(t)$ and ionized hole density $N(z)$ for two values of pump field $E_{0}=30$
and $80~{}\mathrm{kV}\mathrm{cm}^{-1}$, respectively. Only incident field
amplitude $E_{0}$ (without reflected field) shown for $z<0$. Also included are
the nominal internal pump fields $E^{\prime}_{0}=\hat{t_{1}}E_{0}$ (used for
the local MC simulations when specifying $E_{0}$).

One discrepancy remains. It concerns the pump-induced Drude absorption of the
fundamental predicted in the MC-FDTD results (see data for $n=1$ in Fig.
6(a)). While the experimental transmitted fluence continues to grow linearly
in the pump fluence (even showing a superlinear growth at the highest fields),
for the MC-FDTD results one sees that this undergoes a significant saturation
due to the Drude absorption by the ionized carriers. This manifests as a
suppression of the trailing half of the pump pulse exiting the sample, as seen
in the spectrogram in Fig. 5(d). We performed several test simulations to look
for any possible resolution to this discrepancy: (i) reduction of the dopant
density, (ii) allowing for rapid recombination of ionized holes with their
parent ions (both of which would reduce the Drude absorption at $\omega_{1}$)
– however, in all tests the agreement for the overtones was significantly
degraded, with the predicted fluences $F_{n}$ falling further below the
experimental levels. Moreover, we tested that injection of the holes rather
into the hh band does not lead to any significant changes. While the tunnel-
ionization-rate model has not been tested in this regime in the literature,
evidently the onset vs. pump field in the simulations is close to
quantitative. One possible hypothesis is that fewer holes would be generated
than the employed ionization model predicts, and that there is an additional
contribution to the harmonics e.g. from recollision with their parent ions
(not included in the $k$-space MC simulations). As presented in Sec. IV, this
seems unlikely due to the strong intra-cycle scattering, which should strongly
suppress any (coherent) recollision.

Proceeding on the basis that important aspects of the HG process are described
by the MC simulations, we inspect the local response, again using MC
simulations at a single $z$-point, as summarized in Fig. 7. In Fig. 7(a), we
plot the local emission intensity vs. pump field for each harmonic. Also
included for comparison is the scaled fluence from the MC-FDTD simulations
above. One notices immediately that the threshold fields for the latter
(governed in both cases by the onset of tunnel ionization) are significantly
lower (by a factor ${\sim}2$). This is due to the standing-wave field
enhancement effect presented above (Fig. 6(b,c)). With increasing pump fields,
this field enhancement becomes increasingly suppressed, giving rise to a
saturation behavior in the MC-FDTD (and, according to Fig. 6(a), also the
experimental) fluence which is not due to the inherent local response. One
sees also that for the local response, the highest harmonics ($n=7,9$) emerge
rather close to the onset of local saturation effects, which are due to
similar effects for the $T=300$ K case in Sec. III.2, i.e., primarily the
onset of the optical phonon emission. However, as shown before in Fig. 6(a),
the FDTD propagation effects (including phase mismatch which becomes more
severe with increasing $n$) result in relative harmonic yields more consistent
with the experimental data.

In Fig. 7(b) we show the relative emission intensity contribution from hh (the
remaining fraction again being dominated by lh, with a negligible contribution
from the so-band holes). Compared to the $T=300$ K results above (Fig. 4(c),
where the lh contribution was below 20% for all harmonics) here for low pump
fields, one sees that the lh contribution can actually dominate the HG
emission for certain $n$, with a complex dependence on $E_{0}$, although the
hh contribution becomes strongest for all $n$ at the highest fields. To assess
this intriguing result, in Fig. 7(c,d) we plot the local time-dependent hole
populations $N_{b}(t)$ (relative to the dopant density taken as
$N_{d}=5.9\cdot 10^{16}~{}\mathrm{cm}^{-3}$), for both low and higher fields,
respectively. As expected from the implementation of the photoionization
model, at each intra-cycle peak of the pump field one sees a rapid growth in
$N_{\mathrm{lh}}$. Within the subsequent half-cycle, interband scattering
drives the majority of the newly ionized holes into the hh band, such that
during the pulse the relative densities approach those dictated by the density
of states in the respective band. This intra-cycle interband scattering effect
is specific to the case with photoionization, while the modulation of the band
populations with the thermally populated bands ($T=300$ K) is much weaker. As
a control simulation, we also modified the MC treatment to rather inject all
holes into the hh band, which results in the time-dependent band densities
shown in Fig. 7(e). In this case, one sees that the intra-cycle inter-band
scattering is weaker, and (as expected) inverted between lh and hh.
Interestingly, at this high field of $E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$
($E^{\prime}_{0}=37~{}\mathrm{kV}\mathrm{cm}^{-1}$), the predicted harmonic
emission intensities are almost identical to the values for lh ionization,
although the relative lh contribution decreases. Given the highly non-linear
intra-cycle kinetics of $N_{b}(t)$, in the following we inspect more closely
how this might contribute to the HG emission.

Figure 7:  Local MC results ($N_{d}=5.9\cdot 10^{16}~{}\mathrm{cm}^{-3}$,
$T=4$ K). (a) Pump-field dependence of relative emission intensities $V_{n}$
for each odd harmonic $n=1-9$ (solid curves, legend in (b)). Scaled fluence
from MC-FDTD (data from Fig. 6(a), here dotted curves) included for
comparison. As per Fig. 4, the field values $E_{0}$ given correspond to the
external incident fields (see Fig. 6(b)). (b) Fraction of emission intensities
$r_{\mathrm{hh},n}=V_{\mathrm{hh},n}/(V_{\mathrm{hh},n}+V_{\mathrm{lh},n})$
from hh for each $n$ vs. $E_{0}$ ($r_{\mathrm{lh},n}=1-r_{\mathrm{hh},n}$
omitted for visual clarity). (c,d) Time-dependent hole band populations
(relative to $N_{d}$) for two pump fields
$E_{0}=49~{}\mathrm{kV}\mathrm{cm}^{-1}$ and
$81~{}\mathrm{kV}\mathrm{cm}^{-1}$, respectively. Initial photoionized holes
are taken to enter exclusively into the lh band – see text. (e) As per (d),
only using a simulation where initial photoionized holes enter exclusively
into the hh band for comparison. Vertical scale for (c) at left, for (d,e) at
right. Insets show magnified ranges.

One approach to analyze the relative contributions to HG emission from motion
in a given band (“intraband” contribution) with that due to band-population
changes (“interband” contribution), is to calculate the (occupation-weighted)
acceleration components, i.e., with $v(t)=\Sigma_{b}n_{b}(t)v_{b}(t)$ (where
$n_{b}=N_{b}/N_{d}$ is the density of ionized holes in band $b$ relative to
the dopant density) one obtains
$\dot{v}=\partial_{t}v=\Sigma_{b}(\dot{v}_{b}^{\mathrm{(B)}}+\dot{v}_{b}^{\mathrm{(IB)}})$,
where $\dot{v}_{b}^{\mathrm{(B)}}=n_{b}\dot{v}_{b}$ and
$\dot{v}_{b}^{\mathrm{(IB)}}=\dot{n}_{b}v_{b}$. In this case,
$\dot{v}_{b}^{\mathrm{(IB)}}$ comprises both the contributions from
photoionization (“generational nonlinearity”) and from
$\mathrm{hh}\leftrightarrow\mathrm{lh}$ scattering. Although our physical
situation differs from the well-studied HG process from dissociated e-h pairs
[23, 24, 25], this is reasonably consistent with the terminology introduced
there. While this definition tacitly assumes that carriers entering a band
adopt the ensemble average velocity in that band $v_{b}(t)$, in the absence of
a more rigorous treatment one at least observes that
$\dot{v}_{b}^{\mathrm{(IB)}}$ vanishes for $\dot{N}_{b}(t)=0$, while the
occupation weighting provides a direct measure of the relative contribution
from each band. By applying a band-pass filter about each harmonic frequency
$\omega_{n}$, we can inspect their respective contributions to each harmonic.
This is presented in Fig. 8 for both a low (onset of photoionization) and high
(saturation regime) pump fields for both hh and lh bands. Beginning with the
low-field case (Fig. 8(a), $E_{0}=49~{}\mathrm{kV}\mathrm{cm}^{-1}$), several
features can be observed in the results. Firstly, one sees that the lh
contribution is significant, and even dominates the HG emission for $n=5$ (as
seen in Fig. 7(b) for this field value, noting that the emission intensity
scales with $[v(t)]^{2}$). Moreover, for the lh band,
$\dot{v}^{\mathrm{(IB)}}$ is comparable or even somewhat larger than
$\dot{v}^{\mathrm{(B)}}$, depending on $n$. Hence the interband contribution
at low pump fields plays in important role in the HG emission, especially
around the temporal peak ($t=0$ ps) of the pump pulse. In contrast, for the
higher field case (Fig. 8(b), $E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$), the
hh contribution is significantly larger for all harmonics, as is the
contribution from $\dot{v}^{\mathrm{(B)}}$, i.e., the HG emission is dominated
by the EDS of the hh after they have entered the band, similar to the $T=300$
K case (although the lower temperature here also affects the acoustic phonon
scattering rates and equilibrium Fermi-Dirac distribution).

Figure 8:  Local MC results ($N_{d}=5.9\cdot 10^{16}~{}\mathrm{cm}^{-3}$,
$T=4$ K): Decomposition of (occupation-weighted) hole acceleration into
“intraband” ($\dot{v}_{b}^{\mathrm{(B)}}$) and “interband”
($\dot{v}_{b}^{\mathrm{(IB)}}$) contributions (see text for definitions) for
(a) $E_{0}=49~{}\mathrm{kV}\mathrm{cm}^{-1}$ and (b)
$E_{0}=81~{}\mathrm{kV}\mathrm{cm}^{-1}$, for each band
($b=\mathrm{hh},\mathrm{lh}$), filtered for each harmonic $n$. Note different
vertical scaling in (a) and (b).

Finally, we address the evolution of the carrier distribution (both along
$k_{x}$ parallel to the pump field, and transversely along $k_{y,z}$) and the
real-space trajectories of photoionized holes, as predicted by the local-MC
simulations. We consider only the high-field case. In Fig. 9(a,b), we plot the
time evolution of the ensemble-averaged wavevector $k_{x}(t)$ (separately for
hh and lh), along with the rms spreads $\sigma_{x}(t)$ (along $k_{x}$) and
$\sigma_{y,z}(t)$ (along $k_{y,z}$), for both the (a) room-temperature
($T=300$ K, Sec. III.2) and (b) low-temperature ($T=4$ K) situations. For
$T=300$ K, one sees again the relatively broad distribution before the pump
pulse, due to the thermal Fermi-Dirac distribution, as discussed above in
connection with Fig. 4(e,f), whose displacement/spread is only moderately
perturbed during the pump excitation. This is in contrast to the case we
simulated in Ref. [1] for $T=10$ K (assuming a small density of thermally
ionized carriers) where $\sigma_{j}<0.01(2\pi/a)$ (for all $j=x,y,z$) before
the pulse, and grew significantly during the pulse (although remaining
somewhat smaller than the amplitude of $k_{x}$), as here the relative heating
of the hole ensemble is small compared to $k_{B}T$. For the case of $T=4$ K in
Fig. 9(b), with holes generated by photoionization, one sees that the spread
$\sigma_{j}$ rapidly acquires a comparable magnitude as for $T=300$ K, which
is caused by rapid scattering processes even during the first half-cycle of
the pump wave following respective ionization of each hole, with the spread
$\sigma_{y,z}(t)$ pursuing $\sigma_{x}(t)$ closely due to transverse
scattering events. Although we defer the discussion of subsequent recollision
of ionized holes with their parent ions to the next section (see Sec. IV.2),
we show here in Fig. 9(c) results for the distance
$r=|\vec{r}_{h}-\vec{r}_{\mathrm{B}^{-}}|$ between holes and their parent ions
following ionization (608 holes tracked in total). One can see clearly the
coherent bursts of new holes about the field-peak of each pump half-cycle, and
the ballistic initial acceleration to distances $r\sim 50$ nm, which then
becomes diffuse due to momentum scattering, although one can still perceive a
wave-like trend in each half-cycle as the holes accelerate back in the
direction of their parent ions. To inspect for possible recollision, in Fig.
9(d) we plot a vertical zoom of the data in (c). A close inspection (in
particular for the half-cycles in the range of a few ps about $t=0$ where one
has $\sim$100 ionized holes/burst) shows that a small fraction do return to
distances $r<30$ nm, although very few $r<10$ nm, as the transverse momentum
acquired from scattering causes them to pass $x=0$ displaced from $y=z=0$ (the
same applies for holes returning after two half-cycles from the opposite
$x$-direction). Hence scattering is seen to seriously degrade the probability
of a close recollision for the ionized holes. The fact that the first near-
recollisions occur close to a half-cycle after ionization might at first seem
unexpected, as it is well known from the semi-classical treatment of gas-phase
recollision [16] that ionized charges created close to the field peaks,
following the ballistic equation $\partial_{t}v=qE(t)/m$, return after a full
cycle. However, one can show for the case with scattering, e.g. taking a
simple constant damping rate $\Gamma$ ($\partial_{t}v=qE(t)/m-\Gamma v$), that
the holes indeed re-approach their parent ions after a half-cycle for
$\Gamma\ll\omega_{1}$ (see Supplementary).

Figure 9:  Time-dependence of $k$-space distributions from local-MC
simulations for (a) $T=300$ K and (b) $T=4$ K, for both hh (left panels) and
lh (right panels), including ensemble average $k_{x}(t)$ (parallel to pump
field), and rms spreads $\sigma_{x}(t)$ (along $k_{x}$) and
$\sigma_{y}(t)=\sigma_{z}(t)$ (along $k_{y,z}$). In (b), data truncated at
early times before sufficient photoionized holes exist to perform statistics.
(c) Time-dependent distance of hole from parent ion for the case in (b) (hh
and lh combined), obtained by integrating $v_{g}(t)$ for each particle (608
holes in total, data only for the first 3 ps after respective ionization of
each hole, downsampled to a time step of 50 fs for the first 500 fs, and 100
fs thereafter). Color coding for each burst of photoionized holes during each
pump half-cycle. (d) Vertical zoom of (c) to allow inspection of small number
of holes returning to parent ion during subsequent half-cycles

##  IV Discussion

###  IV.1 Comparison with harmonic generation in graphene

We first address the magnitude of the band nonlinearities determined for the
experiments at $T=300$ K (Sec. III.2, Fig. 3(c)). Here it is instructive to
compare these values of $\chi^{(n)}$ with those determined recently for
graphene [2, 26], where the HG mechanism was attributed to an intra-cycle
instability in the Drude heating/absorption of the carrier distribution.
Taking the effective thickness of the monolayer as $L_{g}=0.3$ nm, the authors
deduced values $\chi^{(n)}$ (in respective SI units,
$~{}\mathrm{m}^{(n-1)}\,\mathrm{V}^{-(n-1)}$) of $1.7\cdot 10^{-9}$ ($n=3$),
$1.2\cdot 10^{-22}$ ($n=5$) and $1.7\cdot 10^{-38}$ ($n=7$). While the value
of $\chi^{(3)}$ is 4 orders of magnitude higher than our present case, this is
essentially due to the fact that graphene is an extremely potent electronic
system, where the 2D carrier density in [2] (due to substrate-induced
electrostatic p-doping) of $N_{2D}=2.1\cdot 10^{12}~{}~{}\mathrm{cm}^{-2}$ is
concentrated in a single monolayer. If we normalize the nonlinear coefficients
to the 3D carrier densities $N_{d}$ (Si:B) and $N_{g}=N_{2D}/L_{g}$
(graphene), one arrives at $\chi^{(3)}/N_{d}=1.7\cdot
10^{-35}~{}\mathrm{m}^{5}\,\mathrm{V}^{-2}$ and $\chi^{(3)}/N_{g}=2.4\cdot
10^{-35}~{}\mathrm{m}^{5}\,\mathrm{V}^{-2}$, i.e., the nonlinear response per
carrier is almost equal. One also notes that the ratios between successive
$\chi^{(n)}$ values ($n=3:5:7$) are also loosely correlated, in SI units $\sim
1:10^{-14}:5\cdot 10^{-14}$ for Si:B, and $\sim 1:7\cdot 10^{-14}:10^{-16}$
for graphene. It is hence not correct to conclude that the unique band
structure of graphene would lead directly or indirectly to a stronger
nonlinear response per charge carrier as compared to conventional
semiconductors. While beyond the scope of the current work, this comparison
raises the interesting question whether a more universal sum-rule [27, 28] may
apply for the (odd-order) THz nonlinear susceptibilities of charges in solid-
state bands.

###  IV.2 Contribution of recollisions to the nonlinear response?

We turn now to the low-temperature HG process with photoionization. As
mentioned in Sec. III.3, one can compare our scenario with that of the more
well-established HG process observed for electron-hole pairs in solids with an
energy bandgap, following interband excitation via either tunnel ionization
[29, 23, 24, 30] or optical pre-excitation [25, 31]. These studies provide
strong support for recollision as a dominant mechanism for HG emission (in
that context referred to as “inter-band” generation), despite the presence of
scattering processes, although the role of the subsequent anharmonic “intra-
band” currents are also proposed to dominate, at least in certain regimes
[32].

In the case with pump fields at mid-infrared or higher frequencies (and hence
carrier periods $<50$ fs), it seems reasonable that such a coherent
recollision can occur for a significant fraction of the e-h pairs.
Nevertheless, even for recollision of e-h pairs (arising from tunnel
ionization of optically excited excitons) in GaAs/AlGaAs quantum wells [25]
with a 570-GHz driving field, the authors deduce that LO phonon emission there
does not destroy the recollision process as the required kinetic energy
($\mathcal{E}_{op}{}=36$ meV) is only acquired during a short time directly
before recollision. Note that this is based on the assumption that the
dominant trajectory for a given emitted photon energy $h\nu_{em}$ corresponds
to complete e-h recombination, i.e., $h\nu_{em}=I_{p}+\mathcal{E}_{r}$, where
the kinetic energy upon recollision $\mathcal{E}_{r}$ is governed by their
birth time in the driving field. This assumption is deduced from the condition
for coherent emission in the purely ballistic case, as in dilute gases [16],
and may well need revision when stochastic scattering processes occur such
that the trajectories are no longer deterministic.

Our photoionization+MC results (Sec. III.3, Fig. 9(c)) suggest that the
scattering processes occurring during the hole trajectories severely disrupt
the return paths to their parent ions (both in terms of proximity and
coherence). A simple treatment of the scattering (see Supplementary) also
indicates a decisive influence of scattering, as it results in a kinetic
energy $\mathcal{E}_{r}$ during any residual recollision events lower than
$0.1\cdot U_{p}$ which is compared to $3.17\cdot U_{p}$ for the ballistic case
($U_{p}$ being again the ponderomotive energy, see 2nd paragraph of Sec.
III.3). We note in passing, that this should not result in a high-frequency
cut-off of HG in the frequency range covered in this study
($4~{}\mathrm{THz}$), as that cut-off in the emission photon energy is given
by $I_{p}+\mathcal{E}_{r}$, where the ionization energy
$I_{p}=45~{}\mathrm{meV}\hateq 10.9~{}\mathrm{THz}$ is already above the
covered photon energy range.

Nevertheless, given that the MC simulations here underestimate the HG emission
strength (Fig. 6(a)), further theoretical studies should aim at quantifying
any residual HG from recollision. More generally, the results here strongly
motivate future efforts to refine the description of the photoionization and
any interactions with the parent ions. While we have employed the literature
static-field ionization rate and initially populated only the lh band [17,
18], this approach should be compared to time-dependent quantum-mechanical
treatments, e.g. propagating the time-dependent density-matrix (or
semiconductor-Bloch) equations [33, 34, 31, 35]. Such simulations may require
a detailed treatment of both (i) the excited bound acceptor states [36] if
these are involved as intermediate states during ionization [33, 37], and (ii)
the subsequent scattering processes occurring in the bands (where our MC
treatment here does include explicitly effects such as the angular dependence
and Pauli blocking).

##  V Conclusion

In conclusion, we have presented a combined experimental and theoretical study
of harmonic generation in Si:B with multi-cycle high-field THz pump pulses,
for both temperature regimes where the band holes are either initially
thermally ionized ($T=300$ K), or photoionized during the pulse ($T=4$ K).
Pumping at 300 GHz, we observed up to the 9th harmonic order at room
temperature and up to the 13th order at cryogenic temperature. Near
quantitative agreement with experiment was achieved with Monte-Carlo
simulations, but only when one includes both lh and hh bands and treats the
propagation effects rigorously. This agreement allowed us elucidate the
microscopic harmonic generation processes on the basis of simulations of the
local dynamics, which showed that the nonlinear response of pre-existing holes
are dominated by energy-dependent scattering (and not band non-parabolicity)
in this excitation regime, whereas for photoionized holes, the initial inter-
band scattering processes also play an important role, especially close to the
ionization threshold. We find that scattering during the first sub-cycle
should strongly influence the trajectories of photoionized holes (and any
coherent recollision with their parent ions), which strongly motivates further
studies of the photoionization process and subsequent interactions with the
parent ions for the case of harmonic generation in solids. Here, THz pump THz-
probe (and/or photocurrent-probe) experiments would also be invaluable.

The nonlinear susceptibility per charge carrier was shown to be comparable to
that for harmonic generation in graphene. The latter’s Dirac-type band
structure does apparently not lead to a higher per-carrier nonlinear response
than that found in our doped Si. For practical applications, if one creates a
(micrometer) thin layer of dopants on Si surface by using ion implantation,
one would expect a high harmonic generation from this doped layer, which might
open a new platform for future Si-based free-space or on-chip frequency
multipliers or frequency-comb generators. In addition, if one further employs
the waveguides or resonators, one might even generate high harmonics with an
electric pump field of less than $10~{}\mathrm{kV}\mathrm{cm}^{-1}$.

We finally note that this work serves to validate the application of MC
simulations at this level of description for even higher frequencies and field
strengths than it was previously employed, which will become increasingly
relevant for the description of future (opto-)electronic devices.

## Acknowledgments

We acknowledge funding by the German Research Foundation (DFG) under the
contract RO 770/41 and via the Collaborative Research Center TRR 288
(422213477, project B08). Parts of this research were carried out at the ELBE
Center for High-Power Radiation Sources Sources at the Helmholtz-Zentrum
Dresden - Rossendorf e. V., a member of the Helmholtz Association; for more
information about the facility, see DOI: 10.17815/jlsrf-2-58. The authors
thank Alexej Pashkin for his help for the experiments. We thank the ELBE team
for the operation of the TELBE facility.

## References

  * Meng  _et al._ [2022] F. Meng, F. Walla, Q. ul Islam, A. Pashkin, H. Schneider, C. Jungemann, M. D. Thomson, and H. G. Roskos, Importance of valence-band anharmonicity and carrier distribution for third- and fifth-harmonic generation in Si:B pumped with intense terahertz pulses, Phys. Rev. B 106, 075203 (2022).
  * Hafez  _et al._ [2018] H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, Extremely efficient terahertz high-harmonic generation in graphene by hot dirac fermions, Nature 561, 507 (2018).
  * Deinert  _et al._ [2020] J.-C. Deinert, D. A. Iranzo, R. Pérez, X. Jia, H. A. Hafez, I. Ilyakov, N. Awari, M. Chen, M. Bawatna, A. N. Ponomaryov, S. Germanskiy, M. Bonn, F. H. Koppens, D. Turchinovich, M. Gensch, S. Kovalev, and K.-J. Tielrooij, Grating-graphene metamaterial as a platform for terahertz nonlinear photonics, ACS Nano 15, 1145 (2020).
  * Kovalev  _et al._ [2020] S. Kovalev, R. M. A. Dantas, S. Germanskiy, J.-C. Deinert, B. Green, I. Ilyakov, N. Awari, M. Chen, M. Bawatna, J. Ling, F. Xiu, P. H. M. van Loosdrecht, P. Surówka, T. Oka, and Z. Wang, Non-perturbative terahertz high-harmonic generation in the three-dimensional dirac semimetal cd3as2, Nat. Commun. 11 (2020).
  * Dessmann  _et al._ [2021] N. Dessmann, N. H. Le, V. Eless, S. Chick, K. Saeedi, A. Perez-Delgado, S. G. Pavlov, A. F. G. van der Meer, K. L. Litvinenko, I. Galbraith, N. V. Abrosimov, H. Riemann, C. R. Pidgeon, G. Aeppli, B. Redlich, and B. N. Murdin, Highly efficient THz four-wave mixing in doped silicon, Light Sci. Appl. 10 (2021).
  * Patel  _et al._ [1966] C. K. N. Patel, R. E. Slusher, and P. A. Fleury, Optical nonlinearities due to mobile carriers in semiconductors, Phys. Rev. Lett. 17, 1011 (1966).
  * Mayer and Keilmann [1986] A. Mayer and F. Keilmann, Far-infrared nonlinear optics. ii. $\chi^{(3)}$ contributions from the dynamics of free carriers in semiconductors, Phys. Rev. B 33, 6962 (1986).
  * Urban  _et al._ [1995] M. Urban, C. Nieswand, M. R. Siegrist, and F. Keilmann, Intensity dependence of the third-harmonic-generation efficiency for high-power far-infrared radiation in n-silicon, J. Appl. Phys. 77, 981 (1995).
  * Wang and Ressler [1969] C. C. Wang and N. W. Ressler, Nonlinear optical effects of conduction electrons in semiconductors, Phys. Rev. 188, 1291 (1969).
  * Brazis  _et al._ [1998] R. Brazis, R. Raguotis, and M. R. Siegrist, Suitability of drift nonlinearity in Si, GaAs, and InP for high-power frequency converters with a 1 THz radiation output, J. Appl. Phys. 84, 3474 (1998).
  * Meng  _et al._ [2020] F. Meng, M. D. Thomson, Q. ul Islam, B. Klug, A. Pashkin, H. Schneider, and H. G. Roskos, Intracavity third-harmonic generation in Si:B pumped by intense terahertz pulses, Phys. Rev. B 102, 075205 (2020).
  * Dargys and Kundrotas [1996] A. Dargys and J. Kundrotas,  _Handbook on the physical properties of Ge, Si, GaAs and InP_ , Vol. 39 (Lithuania: Science and Encyclopedia Publishers, 1996).
  * Jungemann  _et al._ [1999] C. Jungemann, S. Keith, M. Bartels, and B. Meinerzhagen, Efficient full band monte carlo hot carrier simulation for silicon devices, IEICE Trans. Electron. E82-C, 870 (1999).
  * Jungemann and Meinerzhagen [2003] C. Jungemann and B. Meinerzhagen,  _Hierarchical Device Simulation: The Monte-Carlo Perspective_ , Computational Microelectronics (Wien, New York: Springer, 2003).
  * Jungemann  _et al._ [2022] C. Jungemann, F. Meng, M. D. Thomson, and H. G. Roskos, Massively parallel FDTD full-band Monte Carlo simulations of electromagnetic THz pulses in p-doped silicon at cryogenic temperatures, Solid State Electron. 197, 108439 (2022).
  * 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, Phys. Rev. A 49, 2117 (1994).
  * Dargys and Žurauskas [1995] A. Dargys and S. Žurauskas, Tunnel ionization of shallow acceptors and donors in GaAs, J. Phys.: Condens. Matter 7, 2133 (1995).
  * Nie and Coon [1984] H. Q. Nie and D. Coon, Tunneling of holes from acceptor levels in an applied field, Solid State Electron. 27, 53 (1984).
  * Yudin and Ivanov [2001] G. L. Yudin and M. Y. Ivanov, Nonadiabatic tunnel ionization: Looking inside a laser cycle, Phys. Rev. A 64, 013409 (2001).
  * Reiss [2008] H. R. Reiss, Limits on tunneling theories of strong-field ionization, Phys. Rev. Lett. 101, 043002 (2008).
  * Boroumand  _et al._ [2022] N. Boroumand, A. Thorpe, A. M. Parks, and T. Brabec, Keldysh ionization theory of atoms: mathematical details, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 213001 (2022).
  * Ganichev  _et al._ [1986] S. D. Ganichev, A. P. Dmitriev, S. A. Emel’yanov, Y. V. Terent’ev, I. D. Yaroshetskii, and I. N. Yassievich, Impact ionization in semiconductors under the influence of the electric field of an optical wave, Sov. Phys. JETP 63, 256 (1986).
  * Vampa  _et al._ [2014] G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Theoretical analysis of high-harmonic generation in solids, Phys. Rev. Lett. 113, 073901 (2014).
  * Vampa  _et al._ [2015] G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, Linking high harmonics from gases and solids, Nature 522, 462 (2015).
  * Banks  _et al._ [2013] H. Banks, B. Zaks, F. Yang, S. Mack, A. C. Gossard, R. Liu, and M. S. Sherwin, Terahertz electron-hole recollisions in $\mathrm{GaAs}/\mathrm{AlGaAs}$ quantum wells: Robustness to scattering by optical phonons and thermal fluctuations, Phys. Rev. Lett. 111, 267402 (2013).
  * Hafez  _et al._ [2019] H. A. Hafez, S. Kovalev, K.-J. Tielrooij, M. Bonn, M. Gensch, and D. Turchinovich, Terahertz nonlinear optics of graphene: From saturable absorption to high-harmonics generation, Adv. Opt. Mater. 8, 1900771 (2019).
  * Kuzyk  _et al._ [2013] M. G. Kuzyk, J. Pérez-Moreno, and S. Shafei, Sum rules and scaling in nonlinear optics, Physics Reports 529, 297 (2013), sum Rules and Scaling in Nonlinear Optics.
  * Passos  _et al._ [2021] D. J. Passos, G. B. Ventura, J. M. B. L. dos Santos, and J. M. V. P. Lopes, Nonlinear optical conductivity of a two-band crystal i, Journal of Physics: Condensed Matter 33, 465701 (2021).
  * 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).
  * You  _et al._ [2017] Y. S. You, Y. Yin, Y. Wu, A. Chew, X. Ren, F. Zhuang, S. Gholam-Mirzaei, M. Chini, Z. Chang, and S. Ghimire, High-harmonic generation in amorphous solids, Nature Communications 8, 724 (2017).
  * Wang  _et al._ [2017] Z. Wang, H. Park, Y. H. Lai, J. Xu, C. I. Blaga, F. Yang, P. Agostini, and L. F. DiMauro, The roles of photo-carrier doping and driving wavelength in high harmonic generation from a semiconductor, Nature Communications 8, 1686 (2017).
  * Kemper  _et al._ [2013] A. F. Kemper, B. Moritz, J. K. Freericks, and T. P. Devereaux, Theoretical description of high-order harmonic generation in solids, New Journal of Physics 15, 023003 (2013).
  * Kitamura [2015] H. Kitamura, Spherical-harmonics expansion method for density-matrix simulations of quantum electron dynamics in continuum states, International Journal of Quantum Chemistry 115, 1587 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24988 .
  * Luu and Wörner [2016] T. T. Luu and H. J. Wörner, High-order harmonic generation in solids: A unifying approach, Phys. Rev. B 94, 115164 (2016).
  * Yue and Gaarde [2020] L. Yue and M. B. Gaarde, Imperfect recollisions in high-harmonic generation in solids, Phys. Rev. Lett. 124, 153204 (2020).
  * Buczko and Bassani [1992] R. Buczko and F. Bassani, Shallow acceptor resonant states in si and ge, Phys. Rev. B 45, 5838 (1992).
  * Serebryannikov and Zheltikov [2016] E. E. Serebryannikov and A. M. Zheltikov, Strong-field photoionization as excited-state tunneling, Phys. Rev. Lett. 116, 123901 (2016).