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

We investigate ultrafast harmonic generation (HG) in Si:B, driven by intense pump pulses with fields reaching ~100 kV/cm 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 30 kV/cm, 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 MV 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 ∼nanosecond pulses.Although the estimated peak electric field reached over ∼100 kVcm −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 pdoped bulk silicon in a focused 0.3-THz radiation field reaching maximum incident peak fields of ∼100 kVcm −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 Γ point: E op = 63.3 meV -the TO and LO phonons being degenerate there [12]).

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 ν 1 = 300 GHz, with a pulse duration of T 1 = 14.1 ps (FWHM) and pulse energies up to 1 µ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 ⋅ 10 15 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⋅10 16 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.

A. 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 xdirection, 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 lighthole 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 χ (n) , which originates from the growth of the number of free holes by impact ionization and -setting in at fields of ∼30 kVcm −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.

B. 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 kVcm −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 2 g 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 ∼ 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 timedomain (FDTD) propagation scheme with a full, selfconsistent 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 δz = 1 µ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 ≤ 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.
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 ω n = 2πν 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 simu- lated 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 powerlaw fits E n ∼ E ηn 0 using the lowest available field ranges in each case; the exponent η 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 ν 1 , and amounts to a bleach factor of ∼2 at the highest field (the linear Drude absorption depth of the sample at ν 1 is 230 µm).For n = 3, one observes an initial growth close to the ideal power-law exponent of η 3 = 3, although the data visibly saturate by a factor ∼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) satu-ration 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 χ (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 = β n χ (n) E n 1 with β n = iω 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 χ (n) from the data), we include the factors where α n is the intensity absorption coefficient for each harmonic (the experimental bleach factor is also included to correct α 1 to avoid overestimating χ (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 χ (3) = 1.0 ⋅ 10 −13 m 2 V −2 , which is highly consistent with the value χ (3) = 0.9 ⋅ 10 −13 m 2 V −2 one can deduce from previous reports for Si:B with nanosecond pulses at ν 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π) 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.
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 = hh, lh, so) and occupation-weighted total velocity v(t), wavevectors k bj (t) and their rms spreads σ bj (t) (j = x, y, z, where the pump field is polarized along the x-direction), band populations N b (t) and scattering rates Γ 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 (ω) = F{v(t)} 2 about each ω n .Also included are the fluences F n ∝ 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 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 (E0, 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 χ (n) calculated from experimental peak fields in (a) after correction for Drude absorption (see main text).
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 ω 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 essen-tially 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 ∆k n ∼ ω 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 hh N 0 ≈ 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 hh rose to 0.93 in non-equilibrium, further depleting N 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).
As mentioned in Sec.III A, one can distinguish two main mechanisms for the HG process: band nonparabolicity ("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 occupationweighted wavevector component k bx (t), i.e., K bn obtained from integrating K b (ω) = F{k bx (t)} 2 about each ω n (in contrast, we determined K b (ω) in Ref. [1] from a 1D⟨k⟩ 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 → V b1 , for the overtones a ratio of K bn V bn →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 kVcm −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 ≤ 10 14 cm −3 and somewhat lower fields (E ′ 0 ≲ 25 kVcm −1 , E 0 ≲ 50 kVcm −1 ), where K bn V bn ≪ 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 Γ point with 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 ≲ 50 kVcm −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 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 kVcm −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) ⋅ (1 + σ x (t) k bx (t) ), i.e., k bx (t) extended by the instantaneous rms width σ 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 kVcm −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 lh,n /K lh,1 compared to K hh,n /K 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→hh intra-cycle energy equilibration takes full effect.

C. 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 ≲ 50 kVcm −1 ) and higher frequency (ν 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 χ (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 ∼80 kVcm −1 (and internal fields E ′ 0 ∼40 kVcm −1 ), and more decisively at a significantly lower frequency of ν 1 = 0.3 THz, such that the ponderomotive energy U p,b = q 2 E 2 0 (4m b ω 2 1 ) [16] will be up to a factor ∼50 higher (U p,lh ∼ 1.3 eV for lh with m 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 ) γ b = I p 2U p,b also 7 times higher (γ lh ∼ 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 kVcm −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 ∼ T 1 √ 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 multiphoton 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.
As in Sec.III B, 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 ∂ 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 ing 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 (γ ≪ 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 ∝ 1 m b dependence for the ponderomotive energy.As a simplified approach, we then take w i,hh → 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,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 (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.

(b) and
As the broadband EO field detection was only sensitive enough to resolve the harmonics for E 0 ≳ 40 kVcm −1 , we augmented the measurements of the 3HG-emission with a more sensitive low-bandwidth EO sensor for pump fields down to E 0 ∼10 kVcm −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 = 0 c ∫ 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 ∝ E 2η3 0 with η 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 η 3 = 4.2 was determined, which exceeds the value η 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 η 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 = ∆N N r ∼ 1 − E 2 m E 2 , with the characteristic field constant E m increasing with ω 1 .Hence ∂ E f ∼ +2E 2 m E 3 , which should be larger for the previous experiments with ν 1 = 1.29 THz and hence produce a larger value of η 3 , which is at least qualitatively consistent with the two results.
Turning now to the higher field range E 0 ≳ 25 kVcm −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 pumpfield 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) = Σ 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 fieldenhancement peaks (at each surface and in the center of the Si:B sample), which exceed the field E ′ 0 one would have if only accounting for the Fresnel transmission coefficient t1 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 ′ 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.
One discrepancy remains.It concerns the pumpinduced Drude absorption of the fundamental predicted in the MC-FDTD results (see data for n = 1 in Fig. 6 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 ω 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 pos-sible 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, 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 ∼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 B, 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 ⋅ 10 16 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 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 kVcm −1 (E ′ 0 = 37 kVcm −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.
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 b comprises both the contributions from photoionization ("generational nonlinearity") and from hh ↔ lh scattering.Although our physical situation differs from the wellstudied 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 v(IB) b vanishes for Ṅ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 ω 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 kVcm −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, v(IB) is comparable or even somewhat larger than v(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 kVcm −1 ), the hh contribution is significantly larger for all harmonics, as is the contribution from v(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).
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 σ x (t) (along k x ) and σ y,z (t) (along k y,z ), for both the (a) room-temperature (T = 300 K, Sec.III B) 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 σ j < 0.01(2π 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 with holes generated by photoionization, one sees that the spread σ 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 σ y,z (t) pursuing σ 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 B), we show here in Fig. 9(c) results for the distance r = ⃗ r h − ⃗ r 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 ∼ 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 ∼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 ioniza-tion 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 ∂ 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 Γ (∂ t v = qE(t) m − Γv), that the holes indeed re-approach their parent ions after a half-cycle for Γ ≪ ω 1 (see Supplementary).

A. 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 B, Fig. 3(c)).Here it is instructive to compare these values of χ (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 χ (n) (in respective SI units, m (n−1) V −(n−1) ) of 1.7 ⋅ 10 −9 (n = 3), 1.2 ⋅ 10 −22 (n = 5) and 1.7 ⋅ 10 −38 (n = 7).While the value of χ (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 pdoping) of N 2D = 2.1 ⋅ 10 12 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 χ (3) N d = 1.7 ⋅ 10 −35 m 5 V −2 and χ (3) N g = 2.4 ⋅ 10 −35 m 5 V −2 , i.e., the nonlinear response per carrier is almost equal.One also notes that the ratios between successive χ (n) values (n = 3 ∶ 5 ∶ 7) are also loosely correlated, in SI units ∼ 1 ∶ 10 −14 ∶ 5⋅10 −14 for Si:B, and ∼ 1 ∶ 7⋅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.

B. Contribution of recollisions to the nonlinear response?
We turn now to the low-temperature HG process with photoionization.As mentioned in Sec.III C, one can compare our scenario with that of the more wellestablished HG process observed for electron-hole pairs in solids with an energy bandgap, following interband excitation via either tunnel ionization [23,24,29,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 "interband" 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 (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ν em corresponds to complete e-h recombination, i.e., hν em = I p + E r , where the kinetic energy upon recollision 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 C, 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 E r during any residual recollision events lower than 0.1 ⋅ U p which is compared to 3.17 ⋅ U p for the ballistic case (U p being again the ponderomotive energy, see 2nd paragraph of Sec.III C).We note in passing, that this should not result in a highfrequency cut-off of HG in the frequency range covered in this study (4 THz), as that cut-off in the emission photon energy is given by I p + E r , where the ionization energy I p = 45 meV ≙ 10.9 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 [31,[33][34][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 kVcm −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.

FrequencyFIG. 1 .
FIG. 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→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 E0 = 115 kVcm −1 ) and T = 4 K (E0 = 81 kVcm −1 ) (linear transmitted fundamental spectrum with undoped sample included for comparison).

FIG. 2 .
FIG. 2. (a) Experimental transmitted temporal field for a Si:B sample (N d = 5.75 ⋅ 10 15 cm −3 , T = 300 K) for an incident pump peak field of E0 = 115 kVcm −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 E0 = 100 kVcm −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 Tn ∼ T1 √ n, curves shown for Tn = ±2T FWHM ).

FIG. 3 .
FIG. 3. (a) Pump-field dependence of transmitted harmonic peak fields (N d = 5.75 ⋅ 10 15 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 (E0, 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 χ (n) calculated from experimental peak fields in (a) after correction for Drude absorption (see main text).

FIG. 4 .
FIG. 4. Local MC results (N d = 5.75 ⋅ 10 15 cm −3 , T = 300 K).(a) Pump-field dependence of relative emission intensities Vn 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 ′ 0 = t1E0 (where t1 = 0.46 is the Fresnel field transmission entering the sample).(b) Time-domain (ensemble-average) velocity for E0 = 110 kVcm −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 (E0 = 55, 110 kVcm −1 , respectively) plotted as loci vs. the ensemble averages k bx (t).Also included are the same loci, but vs. k bx ⋅ (1 + σ 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 Eop = 63.3 meV) included as vertical lines.Orientation of temporal hysteresis indicated by arrow in (f).

FIG. 5 .
FIG. 5. (a) Experimental transmitted temporal field for Si:B sample (N d = 5.0 ⋅ 10 16 cm −3 , T = 4 K) for an incident pump peak field of E0 = 81 kVcm −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 E0 = 80 kVcm −1 (N d = 5.9 ⋅ 10 16 cm −3 ).Both spectra plotted with the same absolute intensity scale.(c) Experimental and (d) MC-FDTD spectrograms.Additional annotations per Fig. 2.

FIG. 7 .
FIG. 7. Local MC results (N d = 5.9 ⋅ 10 16 cm −3 , T = 4 K).(a) Pump-field dependence of relative emission intensities Vn 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 E0 correspond to the external incident fields (see Fig. 6(b)).(b) Fraction of emission intensities r hh,n = V hh,n (V hh,n + V lh,n ) from hh for each n vs. E0 (r lh,n = 1 − r hh,n omitted for visual clarity).(c,d) Timedependent hole band populations (relative to N d ) for two pump fields E0 = 49 kVcm −1 and 81 kVcm −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.

FIG. 9 .
FIG. 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 kx(t) (parallel to pump field), and rms spreads σx(t) (along kx) and σy(t) = σz(t) (along ky,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 vg(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