Introduction

High harmonic generation (HHG) is a non-perturbative nonlinear optical process during which strong laser fields produce harmonic frequencies of 5th order or above. HHG in atoms has been the primary source of coherent X-ray radiation in tabletop laboratory setups. The physical origin of atomic HHG can be explained within the framework of a three-step model - tunneling ionization of electrons to vacuum, propagation in the electric field of the driving light wave and re-collision with the parent ion1,2. A similar model may also be applied to describe the main features of HHG in periodic solids. The main difference from the atomic case is that the electron (hole) is excited to conduction (valence) bands with nonparabolic dispersion. When an electron-hole wavepacket is coherently driven by harmonically oscillating field, its motion in an anharmonic potential leads to generation of intraband currents with high-frequency spectral components, which are the first source of high frequency photons. The second source of high harmonic radiation comes from the recombination of accelerated electron and hole and it is related to the coherent interband polarization. Besides the electronic dispersion, the solid-state HHG differs from generation in atoms by the fact that the electrons may scatter during the propagation due to electron–phonon or electron–electron interactions. Comparison of experimental HHG spectra with theoretical calculations suggests that the momentum scattering times are very short for high energy electrons and may influence the resulting emission spectra3,4.

From its first observation in a solid-state crystal (ZnO)5, HHG has been investigated in a wide range of materials6,7,8,9,10,11,12 including 2D materials such as monolayers of transition metal dichalcogenides7,13, graphene8,14,15 and other exotic materials, e.g., epsilon-near-zero indium tin oxide films16, AlN17, TiN18, or topological BiSbTeSe219. In general, both the interband polarization and intraband currents contribute to the HHG process and the relative strengths of their contributions depend on the material and the experimental conditions7,15,20.

Solid-state HHG shows a strong dependence on the polarization state and waveform of the driving field. After the tunneling excitation, the excited electron and hole propagate in the crystal far from their equilibrium positions. Their classical trajectories are determined by the superposition of the vector potential of the driving light wave and the potential generated by other electrons and ions in the crystal. When the trajectory crosses neighboring lattice sites, the HHG yield may be strongly enhanced20,21,22. Besides the polarization of the driving light, the trajectories of the electrons can be controlled, e.g., by applying a coherent combination of fundamental and second harmonic frequencies23,24. In combination with polarization control, such two-color fields allow to study the chiral properties of HHG25. The HHG spectra further depend on the exact waveform of ultrashort laser pulses, i.e. the carrier-envelope phase (CEP). Such dependence can serve for tuning the HHG spectra9,13,19,26 or for measuring the CEP of pulses with a limited bandwidth27. The CEP-dependent frequency shift of the harmonics varies depending on the material and can be used to obtain information about the time structure of the emitted harmonics11,28,29.

One of the most prominent directions in HHG spectroscopy of solids is the research of the connection between the spectral and polarization properties of the emitted photons and the electronic structure of the studied materials. As shown recently, HHG can serve for both the reconstruction of the electronic potential in ionic materials30 and for spectroscopy of Van Hove singularities in the joint density of states of wide band gap crystals24. HHG in crystalline silicon has been studied already in20,26,31,32,33, but a clear connection between the spectrally-dependent HHG enhancement for particular polarization directions and the material band structure has not yet been drawn. In this paper we demonstrate that the polarization anisotropy of the HHG spectra generated in silicon depends on the frequency of the driving light. At low frequencies (mid-infrared spectral region) leading to a large displacement of the carriers in both the real and reciprocal spaces, the polarization anisotropy changes as a function of the generated photon energy. We show that this observation is closely related to the critical points in the band structure corresponding to Van Hove singularities in the joint density of states. However, when driven with higher frequency field (near-infrared spectral region), the polarization anisotropy of HHG radiation is independent of the generated photon energy, which points out its different origin. The polarization-dependent generation yield is in this case related to the anisotropy of the reduced mass of electron-hole pair at the minimum band gap position in the Γ-point of the Brillouin zone, which influences the excitation probability. We also show that the HHG spectra change with the CEP of ultrashort driving pulses and that the modulation depends on the polarization direction of the strong infrared wave with respect to the crystallographic orientation of silicon. The experimental results are compared with numerical simulations performed using time-dependent density functional theory (TDDFT).

Results and discussion

The non-perturbative high harmonic generation in crystalline silicon is studied in reflection geometry34 to avoid propagation effects in the sample volume. The sample with the surface normal in the <001> crystallographic direction is irradiated by linearly polarized ultrashort laser pulses from two distinct frequency regions. The mid-infrared driving pulses (further referred as low-frequency) with the duration of 25 fs (15 fs for the measurements of the dependence of HHG spectra on the CEP), the center wavelength of 2.06 μm (photon energy of 0.6 eV) and the peak electric field inside the material of up to 2.5 GV/m generate harmonics up to the 17th order corresponding to photon energy of about 10.5 eV. In the near-infrared spectral region (further referred as high-frequency), the pulses with the duration of 6 fs, spectra covering the region 650–950 nm and field of up to 2.8 GV/m produce harmonics up to the 7th order with the photon energy of 11 eV.

The properties of HHG spectra are characterized as a function of the direction of linear polarization of the driving field with respect to the crystal orientation and the CEP of the compressed low-frequency pulses (the spectrum and temporal profile of the compressed pulses are shown in Supplementary Fig. 3). The sample is placed on a rotation stage to make sure that the measured polarization dependence of the HHG yield is not influenced by the polarization-dependent efficiency of the collection and detection setups. Due to the fact that silicon is a centrosymmetric crystal and we use close-to-normal incidence of the linearly polarized light with respect to the surface, the polarization of the generated harmonics is linear. The spectra in the spectral region 2.5–6.2 eV are measured using an optical spectrograph with a cooled charged-coupled device camera. The spectral region 5.5–12 eV is measured using vacuum ultraviolet spectrometer and a microchannel plate detector. The spectra are corrected for the efficiency of both setups and normalized for the same height of the peak of 9th harmonics of the low-frequency pulse, which is detected in both spectral regions (for experimental details see Methods and Supplementary Figs. 1 and 2).

The HHG yield of harmonics driven by low-frequency pulses scales approximately with the third power of the peak intensity (see Supplementary Fig. 4 in Supplementary Note 1). This scaling corresponds to the non-perturbative regime of HHG. In contrast, the intensity dependences of the harmonics driven by high-frequency pulses follow the power-law with the corresponding photon order (see Supplementary Fig. 5 in Supplementary Note 1). In Fig. 1a we show the measured HHG spectra generated by the low-frequency pulses (photon energy of 0.6 eV) with linear polarization along [100] (black curve) and [110] (red curve) crystallographic directions, respectively. We observe an intriguing dependence of the HHG yield for different harmonic orders. While the low-order (5th and 7th) and highest order harmonics (15th and 17th) are generated with higher yield for polarization along [110], the harmonics from the central part of the spectrum are most intensive for polarization along [100] direction. The HHG spectra driven by high-frequency pulses (photon energy of 1.65 eV) are shown in Fig. 1b. The anisotropy does not depend on the generated photon energy in this case and all the harmonics are generated with the highest yield with polarization of the driver along [110]. This clearly points out two different contributions to the polarization anisotropy, whose relative strengths depend on the driving frequency. The measured data are compared with the numerical simulations performed using TDDFT in Fig. 1c, d. The TDDFT calculations are performed with the frequency of the incident pulse scaled to the ratio between the numerical and experimental band gaps, which is in this case \({E}_{{{{{{{{\rm{g,num}}}}}}}}}/{E}_{{{{{{{{\rm{g,exp}}}}}}}}}=0.89\). The numerical results are thus shown as a function of photon energy in units corresponding to the width of the first direct band gap of silicon. This procedure allows us to take into account the fact that the band gap energies calculated by density functional theory (DFT) are systematically lower than the experimental value (details of the numerical simulations are described in Methods).

Fig. 1: High harmonic spectra in silicon.
figure 1

a, b Experimental high harmonic spectra (HHG signal) reflected from monocrystalline silicon with linear polarization of the driving pulse along [100] (black curves) and [110] (red curves) crystallographic directions. The central photon energies of the driving pulses are a 0.6 eV (wavelength of 2000 nm) and b 1.65 eV (wavelength of 750 nm). The peak intensity is I0 = 3 TW/cm2 in both a and b. c, d High harmonic spectra for both polarization directions of the driving pulses and the same photon energies as used in a and b calculated using time-dependent density functional theory with Tran-Blaha (TB09) functional.

Full size image

The normalized HHG yield of individual harmonic orders labeled by dotted lines in Fig. 1a are plotted in Fig. 2a for low-frequency excitation and (b) for high-frequency excitation as a function of the angle of linear polarization of the driving pulses θ. For both driving photon energies, the relative change of the yield with the polarization angle is weaker for lower-order harmonics and increases for higher harmonic orders. Similar effect has previously been observed in other materials21,24,35.

Fig. 2: Dependence of high harmonic generation in silicon on the orientation of linear polarization and the electric field amplitude of the driving pulses.
figure 2

a, b Measured dependence of the high harmonic generation (HHG) yield of different harmonic orders on the angle of linear polarization with respect to [100] crystallographic direction (see inset) for a low-frequency (wavelength of 2000 nm, peak vacuum intensity of 3 TW/cm2) and b high-frequency (750 nm, 3 TW/cm2) driving frequencies. c Calculated dispersion of energy gaps of silicon \({E}_{c}^{i}({{{{{{{\bf{k}}}}}}}})-{E}_{v}^{j}({{{{{{{\bf{k}}}}}}}})\) for different combinations of the first 8 conduction and 4 valence bands. Color coding is used to relate the bands with critical points corresponding to Van Hove singularities in the joint density of states at specific energies for [100] (Γ-X) and [110] (Γ-K) directions to the angular dependences of harmonics shown in a and b. d HHG yield for the four highest harmonic orders as a function of polarization angle and peak electric field inside the sample for mid-infrared driving frequency (low-frequency pulse).

Full size image

For qualitative understanding the polarization anisotropy of harmonic yield in silicon, two sources of the enhancement of the HHG for polarization along specific crystallographic directions can be identified. The first is related to the probability of excitation of an electron-hole wavepacket by nonresonant light, which can be approximated using Keldysh formula36,37 for both the perturbative and the non-perturbative interaction regimes. The theory predicts a dependence of the excitation rate on the reduced mass of the electron-hole pair. When the band structure is anisotropic in the vicinity of the lowest direct band gap, the effective mass changes as a function of the direction of the applied electric field. The dependence of highly-nonresonant carrier excitation rate on the orientation of polarization vector of the driving light was experimentally observed, e.g., in diamond, sapphire or quartz38,39. The polarization anisotropy of electron excitation rate is directly transferred to the HHG signal due to the linear scaling of the emitted power with the number of excited carriers. This mechanism is thus expected to generate the anisotropy of the yield of all the harmonic frequencies, which is in agreement with our observation with high-frequency driving pulses. For the low-frequency pulse, this mechanism is clearly not the main source of the observed polarization anisotropy of the HHG yield, which changes as a function of the generated photon energy. The reduced mass anisotropy of the light hole-electron transition in the Γ-point of the silicon Brillouin zone obtained from DFT band structure calculations is \({m}_{[100]}^{* }/{m}_{[110]}^{* }=1.3\). Application of the Keldysh formula36 for the nonresonant excitation rate with the field amplitude of 2 GV/m gives the ratio between the electron populations excited by polarizations along [100] and [110] directions of W[100]/W[110] = 0.6 for the low-frequency driving pulse and W[100]/W[110] = 0.7 for the high-frequency pulse, which is close to the observed HHG anisotropy at low photon energies represented by blue curves in Fig. 2a, b.

The second source of the angular anisotropy of the HHG yield, which is dominant at low driving frequencies allowing the electrons and holes to reach larger momenta, is related to the presence of Van Hove singularities in the joint density of states. These singularities are connected to the critical points of the band structure and correspond to the regions of the Brillouin zone, in which the derivative of the energy difference between two bands with respect to quasimomentum is zero40,41. This condition can be written as \({\nabla }_{{{{{{{{\bf{k}}}}}}}}}({E}_{c}^{i}({{{{{{{\bf{k}}}}}}}})-{E}_{v}^{j}({{{{{{{\bf{k}}}}}}}}))={\nabla }_{{{{{{{{\bf{k}}}}}}}}}{\epsilon }_{g}^{ij}({{{{{{{\bf{k}}}}}}}})=0\), where k is the quasi-momentum and \({E}_{c}^{i}({{{{{{{\bf{k}}}}}}}})\) and \({E}_{v}^{j}({{{{{{{\bf{k}}}}}}}})\) are the i-th conduction and j-th valence band energies. As shown in3,24, the semi-classical description of the electron dynamics using semiconductor Bloch equations (SBE) combined with application of stationary phase approximation allows to relate the enhancement of interband HHG emission to the presence of Van Hove singularities at specific points of the Brillouin zone (details are discussed in Methods section).

The energies of direct transitions between conduction and valence bands as a function of quasimomentum for monocrystalline silicon are shown in Fig. 2c (results of DFT calculations using TB09 functional), where we see significant differences between the Γ-X [100] and Γ-K [110] directions. The bands with critical points are labeled with the same color coding which is used for the angular dependence of different harmonic orders shown in Fig. 2a, b. In the case of low driving frequency, we clearly see the correspondence between the spectrally-dependent HHG anisotropy and the momentum distribution of the critical points. While there is one critical point at low energy (blue) and several critical points at high energies (red) in the Γ-K direction ([110]), the critical points in Γ-X direction ([100]) are associated with the central region of the HHG spectra shown in Fig. 1a. The two gray curves in Fig. 2c label the bands which also show critical points in [110] direction. The enhancement caused by these singularities at moderate intensity of the driving pulse is weaker compared to the enhancement caused by the critical points in [100] direction in this energy region. However, as we show in Fig. 2d, these singularities start to play role at higher field amplitudes.

The main reason why the Van Hove singularities do not influence the HHG driven by high-frequency pulses is the smaller displacement of the electron-hole wavepacket in the reciprocal space compared to the low-frequency driving field42. The maximum momentum which the electron can reach can be estimated using the Bloch acceleration theorem for the time-dependent wavevector, which can be written as k(t) = k0 − A(t)e/ℏ. Here k0 is the initial electron momentum and A(t) is vector potential of the driving wave. When assuming k0 = 0 (tunneling occurs at minimum band gap in Γ-point), the momentum-space displacement of an electron interacting with oscillating electric field with the amplitude of 2 GV/m can be estimated as \({k}_{\max }(t)\approx 0.6{k}_{{{{{{{{\rm{X}}}}}}}}}\) for the low-frequency driving field and only \({k}_{\max }(t)\approx 0.2{k}_{{{{{{{{\rm{X}}}}}}}}}\) for the high-frequency field, where kX is the momentum in the X-point at the edge of the first Brillouin zone. This means that the Van Hove singularities located in the region \(\left\vert k\right\vert =0.25{k}_{{{{{{{{\rm{X}}}}}}}}}\) and \(\left\vert k\right\vert =0.85{k}_{{{{{{{{\rm{X}}}}}}}}}\) can only be reached in the case of the low-frequency driving pulse, while they do not play any significant role for the high-frequency excitation. We note that the electrons generated by the high-frequency pulse via single-photon indirect transition do not contribute to HHG. The reason is their incoherent nature due to the fact that participation of phonon is required during the indirect transition, which leads to population of many different electronic states which contribute to the emitted radiation with random phases.

To verify the role of the maximum electron-hole momentum for the observed HHG polarization anisotropy we measured the yield of the four highest harmonic orders as a function of the polarization angle and peak electric field of the low-frequency pulse (see Fig. 2d). While at moderate field of 1.8 GV/m we observe the anisotropy shown in Fig. 2a, it changes with the increasing field. At the maximum of about 2.5 GV/m, we see new maxima in the angular dependence of 11th and 13th harmonics for polarization along [110] direction, which are related to the two gray bands in 2a. The enhancement of the 15th harmonic for polarization along [100] at high field amplitude might be related to the fact that four energy degenerate direct transitions are possible in the X-point of the Brillouin zone, which strongly increase the joint density of states at the energy of about 8.6 eV. This point can be reached by the electrons driven by the strongest applied field.

Besides the anisotropy of HHG driven by pulses with random CEP we studied the role of the CEP of the compressed low-frequency pulses (2060 nm, 15 fs) on the HHG anisotropy and spectral shape. The resulting HHG spectra modulated as a function of the CEP for driving field polarization along [110] direction are shown in Fig. 3. Only the low order harmonics are measured due to the limitations of dispersion compensation of compressed pulses in the vacuum setup, where the window and the polarizer introduce additional prolongation of the pulses. We can see clear modulation of the spectra, which is caused by the interference between electron trajectories influenced by the exact waveform of the driving pulse. Similar modulation has been observed for the HHG spectra measured in the transmission geometry in silicon driven by pulses at longer central wavelength26. When we consider the relation between the time evolution of high energy photon emission43 and the spectral properties of the HHG radiation, the coherent emission of photons during each half-cycle of the strong driving field leads to harmonic orders in the frequency domain separated by 2ℏω. When the duration of the pulse envelope is comparable to the oscillation period, the field amplitude in subsequent cycles varies significantly. Because the phase of the emitted photons depends on the dynamics of the electron-hole wavepacket during each half-cycle, the resulting harmonic spectrum is modulated as a function of the CEP. In the frequency domain, the phase of each harmonic order is locked to the phase of the driving field due to the coherence of the HHG process. The HHG phase changes as φn = nφ, where n is the harmonic order and φ is the CEP of the fundamental pulse. When the spectra of two neighboring harmonic orders overlap, the photons in the overlapping region interfere leading to two oscillations of the intensity per 2π CEP shift of the fundamental pulse. An interesting feature relating the CEP dependence of harmonic spectra to the band structure anisotropy is shown in Fig. 3c, d, where we plot the difference between the measured CEP-dependent and CEP-averaged spectra \(\left({{{{{{{{\rm{Sig}}}}}}}}}_{{{{{{{{\rm{CEP}}}}}}}}}-{{{{{{{{\rm{Sig}}}}}}}}}_{ < {{{{{{{\rm{CEP}}}}}}}} > }\right)/{{{{{{{{\rm{Sig}}}}}}}}}_{ < {{{{{{{\rm{CEP}}}}}}}} > }\) (signal is normalized by the CEP-averaged spectra to show the relative modulation depth) for polarization of the driving field along [100] and [110]. We observe clear differences between the two signals. There is a phase shift of the modulation for each harmonics and also the slopes of the modulation fringes, which are visible in Fig. 3c, d, differ for the two polarization directions. The polarization anisotropy of the CEP-dependent HHG modulation can be related to different electron dynamics in the two distinct directions in the reciprocal space42.

Fig. 3: Carrier-envelope phase (CEP) dependence of reflective high harmonic generation (HHG) in silicon.
figure 3

a Comparison of high harmonic spectra generated in crystalline silicon for linear polarization of the driving pulses along [100] direction for two values of CEP shifted by π. b CEP-dependent HHG spectra (color bar shows the spectral power density) with driving polarization along [100] direction. c, d Measured and e, f calculated HHG spectra for two directions of linear polarization of the driving pulse along c, e [100] and d, f [110] directions. The photon energy of the pulse was normalized in the time-dependent density functional theory (local density approximation) calculations to account for the underrated width of the band gap of \({E}_{{{{{{{{\rm{g}}}}}}}}}^{{{{{{{{\rm{LDA}}}}}}}}}=2.56\) eV. The color scale represents the relative difference between spectral power densities of the CEP-dependent and CEP-averaged high harmonic radiation.

Full size image

To correlate the observed polarization anisotropy and the CEP dependence of HHG spectra in silicon with the theory we numerically simulate the light-matter interaction using TDDFT. The non-relativistic Kohn-Sham equation is solved using a real-time and real-space implementation using the Hamiltonian that can be written as:

$$\begin{array}{c}\left[{\left(-\frac{i\hslash }{2{m}_{e}}{{{{{{{{\boldsymbol{\nabla }}}}}}}}}_{{{{{{{{\boldsymbol{r}}}}}}}}}+\frac{\left\vert e\right\vert }{c}{{{{{{{\boldsymbol{A}}}}}}}}\left(t\right)\right)}^{2}+{\hat{v}}_{{{{{{{{\rm{ion}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}}\right)+{\hat{v}}_{{{{{{{{\rm{H}}}}}}}}}\left[n\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)\right]\left({{{{{{{\boldsymbol{r}}}}}}}}\right)+\right.\\ \left.+{\hat{v}}_{{{{{{{{\rm{xc}}}}}}}}}\left[n\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)\right]\left({{{{{{{\boldsymbol{r}}}}}}}}\right)\right]\times {\psi }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)=i\hslash \frac{\partial }{\partial t}{\psi }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}},t\right),\end{array}$$
(1)

where ∇r is the real-space gradient operator, ℏ is the reduced Planck constant, c is the light velocity in vacuum, me is the electron mass at rest, r is real-space coordinate, t is time, \(n\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)\) is the density distribution of electrons and \({\psi }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)\) are the Kohn-Shan orbitals in n-th electronic band and k is the electron wave-vector. In Hartree atomic units and Coulomb gauge, the vector potential \({{{{{{{\boldsymbol{A}}}}}}}}\left(t\right)\) is related to the incident electric field \({{{{{{{\boldsymbol{E}}}}}}}}\left(t\right)\) by the relation

$${{{{{{{\boldsymbol{A}}}}}}}}\left(t\right)=-c\int\nolimits_{-\infty }^{t}{{{{{{{\boldsymbol{E}}}}}}}}\left(t^{\prime} \right)\,dt^{\prime} .$$
(2)

The fields generated by the charge motion are not taken into account in our calculations.

In Eq. (1), \({\hat{v}}_{{{{{{{{\rm{ion}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}}\right)\) describes the potential induced by the presence of atoms of the crystal. The primitive cell of Si was constructed using the method of ab-initio norm-conserving pseudo-potentials44. The non-local part of the pseudo-potential is not shown for simplicity. \({\hat{v}}_{{{{{{{{\rm{H}}}}}}}}}\) is the Hartree term, corrected by the exchange-correlation potential \({\hat{v}}_{{{{{{{{\rm{xc}}}}}}}}}\). In this work, both the local density approximation (LDA) and Tran-Blaha (TB09) functionals were employed. We note that the predictor-corrector method is necessary upon usage of TB09 functional45.

From the time-evolved Kohn-Sham orbitals \({\psi }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\left({{{{{{{\boldsymbol{r}}}}}}}},t\right)\), the time-dependent electrical current \({{{{{{{\boldsymbol{J}}}}}}}}\left(t\right)\) was computed as indicated in ref. 46. The current is Fourier transformed to obtain the HHG spectra. The TDDFT spectra calculated using TB09 functional for orientations along [100] and [110] directions are shown in Fig. 1c, d. The main features of the experimental spectra shown in Fig. 1a, b, namely the HHG polarization anisotropy, its spectral dependence and the cut-off photon energies, are well reproduced by the simulations. In the calculations, the pulse envelope was described by \({\sin }^{4}\) function, which improves the quality of HHG spectra47 and pulse duration was the same as in the corresponding experiments. The photon energy of the driving field was normalized by the ratio between the experimental and theory values of the energies of the first direct band gap in silicon (experimental: 3.4 eV48, numerical value using DFT-TB09: 3.06 eV, using LDA: 2.56 eV) to take into account the underestimation of the band gap energy by DFT. To include the effect of dephasing in the simulations, which suppresses the coherent oscillations present in the calculated current after the driving pulse, we apply the procedure described in Methods section.

TDDFT with LDA is applied to simulate the dependence of the spectra on the CEP of short pulses. The comparison of the calculated spectra with the experimental results is shown in Fig. 3c–f for both directions of polarization of the driving field along [100] and [110]. We can see that the results of the numerical simulations reproduce the oscillations of the HHG spectra with changing CEP and also the spectral shift of interference fringes agrees well with the experimental data. The calculated CEP-dependent spectra also differ for the two polarization directions, but the agreement with the experimental data is limited. This is probably caused by the fact that the simulations were performed with idealized pulse with \({\sin }^{4}\) envelope while the real shape of pulses used in experiments was more complex (see Supplementary Fig. 3).

Conclusions

In summary, we present the observation of reflective high harmonic generation in silicon in the vacuum ultraviolet spectral region up to the photon energy of 10.5 eV driven by ultrashort pulses in two distinct spectral regions. We demonstrate that the polarization anisotropy of the HHG yield has two sources whose relative strengths are related to the experimental conditions. The reduced mass anisotropy is the main source at higher driving frequencies corresponding to a small maximum displacement of the electrons in the reciprocal space. In contrast, the enhancement of HHG yield due to the presence of Van Hove singularities in the joint density of states plays significant role at low frequencies of the driving field, for which the electron reaches higher momenta. Further we observe the polarization anisotropy of the CEP-dependent modulation of high harmonic spectra, which demonstrates the importance of sub-femtosecond electron dynamics during the strong-field nonresonant light-matter interactions. The good qualitative agreement of our results with the TDDFT simulations proves the variability of this method and its applicability for description of light-matter interactions in a wide range of experimental settings.

Methods

Experimental setup

HHG in monocrystalline silicon is driven by mid-infrared (low-frequency) CEP-stabilized laser pulses and near-infrared (high-frequency) ultrashort pulses with random CEP. The low-frequency pulses are produced in a setup based on non-collinear optical parametric amplifier combined with difference frequency generation (NOPA-DFG). The output from NOPA-DFG is further amplified by an optical parametric amplifier (OPA). The setup is pumped by an ytterbium femtosecond laser (Pharos SP, Light Conversion), for more details see49. The resulting pulse parameters are the full-width at half maximum (FWHM) pulse duration of 25 fs (15 fs in CEP-dependent measurements, measured using interferometric third harmonic generation frequency-resolved optical gating, ITH FROG), center wavelength of 2.06 μm and the peak intensity in focus of up to 5 × 1012 W/cm2. For the high-frequency excitation, the driving pulses with the FWHM duration of 6 fs and spectra covering the region 650–950 nm are generated using a broadband NOPA and compressed using the combination of a prism compressor and broadband chirped mirrors (Venteon DCM-7).

Crystalline silicon samples with the surface normal in the <001> crystallographic direction are irradiated by linearly polarized light. The laser beam is tightly focused using off-axis parabolic mirror and the harmonics are collected in the reflection geometry to avoid propagation effects. The angle of incidence is of about 10∘. The sample is placed on a rotation stage to make sure that the polarization dependence is not influenced by the focusing and collection setups (Supplementary Fig. 1). The crystallographic orientation of the sample was measured using Laue X-ray diffraction.

Two separate setups are used for detection of HHG spectra in ultraviolet (UV, λ > 200 nm) and vacuum-ultraviolet (VUV, 80 nm < λ < 300 nm) spectral regions. In the VUV setup, the experiments occur under high vacuum conditions. The HHG radiation is dispersed by a curved diffraction grating (VUV spectrometer EasyLight, H+P Spectroscopy) and detected by a microchannel-plate detector (MCP) with phosphor screen (Photonis VID240/P43/GL), which is imaged by a CMOS camera. In the UV setup, the visible radiation is filtered out by using a UV fused silica prism pair placed in the HHG beam in front of the detection setup. The spectra are measured using a grating spectrograph (Shamrock 163, Andor) with cryo-cooled charge-coupled device (CCD) detector (Andor iDUS 420). Spectra measured by both setups are corrected for their spectral sensitivities using the efficiencies of gratings and CCD detector supplied by the manufacturers. The MCP efficiency curve provided by the manufacturer was used in the VUV to UV spectral region. The two spectra are merged together by normalizing the spectra for the same height of the 9th harmonics peak (photon energy of 5.4 eV), which is detected in both spectral windows. The absolute spectral power density calibration is only approximate but it does not influence the comparison of the HHG yield measured with different orientation of the sample or different values of CEP of the driving pulses.

For measurements of the dependence of HHG spectra on the CEP of few-cycle pulses we broaden the spectrum of the driving pulses by focusing the beam into a gadolinium gallium garnet crystal, where self-phase modulation takes place (see Supplementary Fig. 2). Due to negative group velocity dispersion in the infrared spectral range of the covered by the pulse spectrum, the pulse self-compresses in the crystal and the final pulse duration on the sample is 15 fs. The CEP of the laser beam is controlled using a 4f setup, in which the fundamental beam is dispersed by a BaF2 prism. The prism is shifted with respect to the laser beam by a computer-controlled translation stage. Due to the difference of group and phase velocities in BaF2, variation of the optical path in the material leads to the change of the pulse CEP. The pulses are characterized using third order interferometric frequency-resolved optical gating. The resulting spectrum and temporal profile of instantaneous intensity of the pulses used for CEP-dependent measurements are shown in Supplementary Fig. 2.

Semi-classical description of the enhancement of high harmonic generation due to Van Hove singularities

The reason for an enhancement of the yield of high harmonic generation in a particular range of photon energies due to the presence of Van Hove singularities in the band structure can be qualitatively understood when considering the description of the electron dynamics using semiconductor Bloch equations (SBE) presented e.g. in refs. 3,24. The enhancement is related to an increase of interband current at the frequency corresponding to the singularity. The frequency-dependent interband current in the two-band approximation can be written as:

$${j}_{er}(\omega )=\omega {\int}_{BZ}{{{{{{{\bf{kd}}}}}}}}({{{{{{{\bf{k}}}}}}}}){d}^{3}{{{{{{{\bf{k}}}}}}}}\int\nolimits_{-\infty }^{\infty }dt{e}^{-i\omega t}\int\nolimits_{-\infty }^{t}d{t}^{{\prime} }F({t}^{{\prime} })d({{{{{{{{\boldsymbol{\kappa }}}}}}}}}_{{t}^{{\prime} }}){e}^{-iS({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)-(t-{t}^{{\prime} })/{T}_{2}},$$
(3)

where \({{{{{{{{\boldsymbol{\kappa }}}}}}}}}_{{t}^{{\prime} }}={{{{{{{\boldsymbol{k}}}}}}}}+{{{{{{{\boldsymbol{A}}}}}}}}({t}^{{\prime} })-{{{{{{{\boldsymbol{A}}}}}}}}(t)\), A(t) is vector potential of the driving pulse and k is the crystal momentum. \(S({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)\) is the semiclassical action defined as:

$$S({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)=\int\nolimits_{t{\prime} }^{t}{\epsilon }_{g}\left[{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{A}}}}}}}}(t)+{{{{{{{\bf{A}}}}}}}}(\tau )\right]\,d\tau ,$$
(4)

Eq. (3) can be rewritten to a form which contains a slowly varying term \(g({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)\) and a part with classical action multiplied by the complex exponential oscillating at frequency ω (the frequency of high harmonic photons):

$${j}_{er}(\omega )=\omega {\int}_{BZ}{d}^{3}{{{{{{{\bf{k}}}}}}}}\int\int{{{{{{{\bf{g}}}}}}}}({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t){e}^{-S({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)+i\omega t}dt\,dt^{\prime} .$$
(5)

In the vicinity of Van Hove singularities, the classical action becomes stationary leading to an enhancement of HHG emission due to constructive interference of different trajectories of the re-colliding electron and hole. The largest contribution to HHG therefore comes from points where \(S({{{{{{{\bf{k}}}}}}}},t^{\prime} ,t)\) is stationary24. Using stationary phase approximation (SPA) around these stationary points (kst) one can link the harmonic spectra to the band structure:

$$I(\omega )\propto {\left\vert {j}_{er}\right\vert }^{2}\propto {\omega }^{2}{\left\vert \mathop{\sum}\limits_{{{{{{{{{\bf{k}}}}}}}}}_{st}}\frac{\tilde{{{{{{{{\bf{g}}}}}}}}}({{{{{{{{\bf{k}}}}}}}}}_{st},t^{{\prime} }_{st},{t}_{st}){e}^{-iS({{{{{{{{\bf{k}}}}}}}}}_{st},t^{{\prime} }_{st},{t}_{st})}}{\sqrt{\left\vert {\nabla }_{{{{{{{{\bf{k}}}}}}}}}{\epsilon }_{g}^{ij}({{{{{{{{\bf{k}}}}}}}}}_{st})\right\vert }}\right\vert }^{2},$$
(6)

where \({\epsilon }_{g}^{ij}({{{{{{{\bf{k}}}}}}}})=({E}_{c}^{i}({{{{{{{\bf{k}}}}}}}})-{E}_{v}^{j}({{{{{{{\bf{k}}}}}}}}))\) is the momentum-dependent band gap.

Details of time-dependent density functional theory simulations

The primitive cell of silicon was constructed using two atoms along with non-orthogonal periodic boundaries. Interatomic distance was taken equal to 5.431 angstroms. At initial instant t = 0, ground state was prepared by solving the stationary Kohn–Sham equation. The k-space grid was meshed using 4 shifted grids of 40 × 40 × 40 points distributed uniformly in momentum space. The integration time-step was taken equal to 0.25 atomic unit. The numerical band gap energies obtained by DFT are 2.56 eV with LDA and 3.04 eV with TB09 functionals, respectively50. The experimental value for the first direct band gap of silicon at room temperature is 3.4 eV41. The central frequency of pulses used in the calculations is scaled down by the factor corresponding to the ratio between experimental and numerical band gaps to compensate this difference.

Since our TDDFT simulations do not include energy dissipation and decoherence effects, the time-dependent electrical current obtained from the Kohn-Sham equation51 was damped using the decaying half-cycle of a \({\sin }^{4}\) function52 by writing

$${{{{{{{{\boldsymbol{J}}}}}}}}}_{{{{{{{{\rm{damped}}}}}}}}}\left(t\right)={{{{{{{\boldsymbol{J}}}}}}}}\left(t\right)\times \eta \left(t\right)$$

where

$$\eta \left(t\right)= \,1-H\left(t-{t}_{{{{{{{{\rm{final}}}}}}}}}+{t}_{{{{{{{{\rm{damping}}}}}}}}}\right)+{\sin }^{4}\left[\frac{\pi }{2}\times \frac{t-{t}_{{{{{{{{\rm{final}}}}}}}}}}{{t}_{{{{{{{{\rm{damping}}}}}}}}}}\right] \\ \times H\left(t-{t}_{{{{{{{{\rm{final}}}}}}}}}+{t}_{{{{{{{{\rm{damping}}}}}}}}}\right)\times \left[1-H\left(t-{t}_{{{{{{{{\rm{final}}}}}}}}}\right)\right]$$

and \(H\left(t\right)\) is the Heaviside function. This method avoids adding spurious oscillations to the Fourier spectrum upon damping the electrical current. By setting the time-window analysis (using tfinal) to 1.5 − 1.75 × τ, where τ is the laser pulse duration, the simulated spectra are directly comparable to the acquired experimental ones. The obtained agreement suggests that the temporal window selected for the Fourier analysis may be related to the time for electron decoherence beyond which dissipation should be taken into account53.