1. Introduction

Since the observation of nonperturbative high harmonic generation (HHG) from a solid crystal in 2011 [1], intensive study has been performed to understand the mechanism and to demonstrate its potential applications [2]. Because solids are characterized by high atomic density, crystal symmetry, and a band structure unique to each material, solid HHG is potentially promising for compact short-wavelength light sources [3,4], all-optical band-structure reconstruction [4–6], and studies on nonperturbative physics in condensed phases [7,8]. Engineered nanostructures could form novel platforms for high repetition-rate HHG with a possibility of manipulating wavefront and polarization properties [9–12].

Recent studies on solid HHG have indicated two sources as a radiation mechanism that follows electron-hole generation via Zener tunneling: intraband current leading to bremsstrahlung in Bloch oscillation and interband polarization leading to electron-hole recombination [2,6,8,13,14]. It is naturally expected that the spectral, temporal, and polarization characteristics of HHG differ according to the radiation mechanisms. At the same time, information on the band structure and on the attosecond electron dynamics would be imprinted onto HHG in a different manner depending on the radiation mechanism. Therefore, it is important to gain a deeper understanding of HHG for each radiation mechanism for applications in short-wavelength light sources and high-harmonic (HH) spectroscopy.

There have been several studies on the HH properties for each of the radiation mechanisms. Numerical calculations have shown that HHG originating from the intraband current is free of an attosecond chirp (or atto-chirp) and that HHG originating from the interband polarization is accompanied by an atto-chirp [6,15]. The cutoff and plateau of the HH spectrum have been reported to reflect the maximum bandgap at the Brillouin zone edge for the interband process [5,15] and the shape of the conduction bands for the intraband process [16]. The HH photon energy shift that is sensitive to the carrier envelope phase of the input pulse is explained by the interband model [17], while the pure intraband model does not predict such a shift [18].

To gain a deeper understanding of how the HH properties depend on the generation mechanism, it is essential to investigate them for a material in a broad spectral range, including both below and above the bandgap. A gallium selenide (GaSe) crystal is suited for this purpose because its two-dimensional band structure simplifies the trajectories of charge carriers and because its limited energy spread of the conduction band suppresses the possible impact ionization, which may otherwise hamper the imprinting of the band structure onto HHG. While HHG from GaSe has been intensively studied in the range below the bandgap [19–22], much less is known about HHG in the range above the bandgap.

In this paper, we study HHG from a GaSe crystal in a deep-UV range well above the bandgap, aiming to elucidate the properties of high harmonics emitted by the interband polarization mechanism. By use of an input pulse of a 2.1 µm center wavelength, we observe harmonics extending up to the tenth order (210 nm in wavelength), which corresponds to approximately three times the bandgap energy. The harmonic spectrum represents a clear dependence on the crystal orientation with respect to the input electric field, reflecting the $\bar{6}2m$ point group. The 3D time-dependent density matrix (TDDM) simulations reveal that the interband polarization mechanism dominates over the intraband current mechanism above the material’s bandgap. Most importantly, it is revealed that the dominating interband polarization mechanism can create attosecond bursts with characteristic sine-like waveforms. The successful demonstration of the deep-UV HHG and the numerical findings indicate a great promise for attosecond pulse generation through solid-state HHG as well as form the basis for the potential applications of short-wavelength light sources, HH waveform control, and HH spectroscopy at frequencies above the bandgap.

2. Method

A schematic of the optical setup is shown in Fig. 1(a). A 30-µm-thick GaSe single crystal (c-cut, ɛ-type, EKSMA Optics) is irradiated at normal incidence by infrared (IR) pulses with a center wavelength of 2.1 µm, a temporal duration of 100 fs, and a repetition rate of 1 kHz [12]. Crystalline GaSe has a bandgap energy of 1.98 eV (625 nm in wavelength), and its transparency range extends from 0.65 to 18 µm [23]. The ɛ-type GaSe belongs to the $\bar{6}2\textit{m}$ (D_3h) crystallographic point group, having a honeycomb lattice structure, as shown in Fig. 1(c), with a fourfold (Se-Ga-Ga-Se) layer. Here, we define θ as the angle between the a-axis of GaSe and the electric field of the linearly polarized input pulse. The GaSe crystal is mounted on a rotation holder at normal incidence to allow rotation around the c-axis. The IR pulses are focused by a parabolic mirror to a spot with a diameter of 210 µm at the GaSe crystal, which results in a peak intensity of 0.48 TW/cm² (1.1 V/nm) inside GaSe unless otherwise mentioned. High harmonics are spectrally dispersed by a diffraction grating (2400 line/mm) and measured by a photomultiplier tube (PMT, H10722-113, Hamamatsu Photonics). The Rochon prism is used f polarization-resolved detection.

 

Fig. 1. (a) Schematic of the optical setup. (b) Power spectrum of the input pulse with a center frequency of 140 THz (a center wavelength of 2.1 µm). (c) Top view of the GaSe lattice structure and an input electric field E_in with an orientation angle θ.

Download Full Size | PDF

To study the HH properties above the bandgap, we perform numerical simulations based on a 3D TDDM model combined with the band information obtained by density functional theory (DFT). The electronic structures, band energies, and matrix elements of the momentum operator of GaSe are obtained by DFT calculations using experimental atomic position performed with the Elk code [24]. Here, we employ Perdew-Wang local-density approximation [25] for the exchange-correlation potential. We take into account initially occupied 38 valence bands and 65 unoccupied conduction bands, spanning the energy range of 39 eV, as the dynamical degree of freedom. The first Brillouin zone is sampled at 64×64×12 points in 3D. We perform time-evolution of the one-body density matrix for the Kohn-Sham orbitals (see Supplemental Material of Ref. [21] for detailed information on our TDDM model). We employ velocity gauge coupling to the A-field with a temporal duration of 72 fs and a photon energy of 0.62 eV. A good correspondence with the experiments in the orientation angle dependence, shown in the next section, confirms that the phases of the momentum operator matrix elements have proper values without additional consideration typically made in the length-gauge coupling [26]. Note that the 3D treatment of the Brillouin zone is crucial because electron-hole excitation occurs over the whole Brillouin zone to be driven in the 3D momentum space; it is hard to define a dominant 1D path for an electron-hole pair. The HH spectra are obtained by the Fourier transform of the electric current induced in GaSe. The contribution from each of the intraband and interband mechanisms is separately calculated by using the projection onto instantaneous eigenfunctions. We include the transverse relaxation of crystal electrons with a characteristic time constant of 10 fs (see Appendix A for details). Note that our simulations do not include any propagation effects, such as self-phase modulation, nonlinear wave-mixings, polarization change due to linear/nonlinear birefringence, etc.

3. High harmonic spectra and orientation dependence

Figure 2 shows the HH spectrum of the polarization component parallel to the input electric field measured for an orientation angle of θ = 30 deg. The HH spectrum extends up to the tenth order (a deep-UV wavelength of 210 nm), corresponding to the photon energy three times the bandgap. Note that the input-intensity dependence of the HH yield (shown in Appendix B) indicates that the excitation is in the nonperturbative regime.

 

Fig. 2. The HH spectrum of the polarization component parallel to the input electric field measured for an orientation angle of θ = 30 deg.

Download Full Size | PDF

We investigate how the HH spectrum depends on the orientation of the crystal axis with respect to the input electric field. Figure 3 shows the HH spectrum measured at various orientation angles θ for the HH components (a) parallel and (b) perpendicular to the input electric field. Here, the HH spectral intensity is normalized for each order. For the parallel components shown in Fig. 3(a), the even-order harmonics exhibit a 60° period with nearly zero minima at θ = 60°× n (n: integer). The odd-order harmonics also exhibit a 60° period with only small intensity contrast. For the perpendicular components shown in Fig. 3(b), both the even- and odd-order harmonics exhibit a 60° period with nearly zero minima at θ = 30° + 60°×n for the even-order and at θ = 30°×n for the odd-order (the odd-order harmonics take two nodes within each 60° period). The observed orientation angle dependences are similar to those reported in previous studies on GaSe HHG [20,21], where the harmonics below the bandgap driven by a longer-wavelength pulse were ascribed to be dominated by the intraband current mechanism.

 

Fig. 3. High harmonic spectra at various orientation angles of the input electric field. (a) Parallel and (b) perpendicular components of the experimental HH spectra. (c) Parallel and (d) perpendicular components of the HH spectra simulated with an input electric field of 1.0 V/nm. The HH spectral intensity is normalized for each order.

Download Full Size | PDF

The simulated orientation angle dependences are shown in Fig. 3 for (c) parallel and (d) perpendicular components, where the HH spectral intensity is normalized for each order. We see that the simulated angle dependences agree well with the experimental results in both the periodicities and the node angles. A discrepancy is found for the parallel components of odd-order harmonics. Although much smaller modulation contrast is commonly observed for both the experiments and simulations, the peak angles and valley angles are exchanged between experiments and simulations. Such discrepancy would originate from intrinsically small orientation dependence and the possible difference in the intensity dependence between θ = 0° and 30°. In fact, our simulations show that the peak angles and valley angles can be exchanged depending on the center wavelength and the input field strength. Regardless of the minor disagreement, the overall coincidence between experiments and simulations indicate that the intrinsic orientation dependence is experimentally observed without severe distortion by the linear/nonlinear optical birefringence [26–28]. Here, we note that GaSe belonging to the $\bar{6}2\textit{m}$ point group is optically isotropic within the ab-plane with respect to not only the linear refractive index but also the nonlinear refractive index originating from the third-order nonlinearity (see, for example, Ref. [29]). The orientation angle dependence is important because it provides deeper insight into carrier motions in momentum space [21] and enables waveform control of HH pulses using crystal symmetry [20].

Figure 4(a) shows the simulated HH spectra for the parallel components at θ = 30° with separated contributions from the intraband and interband mechanisms. Here, we note that the odd-order harmonics appear as parallel components at θ = 30°, as shown in Fig. 3(c) and (d). Both the intraband and interband components exhibit an exponentially decreasing tendency with increasing order. Importantly, the interband components decrease less steeply than the intraband components and are dominant at frequencies above the bandgap of GaSe (∼480 THz), (similar discussion is found for HHG from GaAs in Ref. [30]). Figure 4(b) shows the simulated HH spectra for the perpendicular components at θ = 0° with separated intraband/interband mechanisms. We see that the even-order harmonics appearing as perpendicular components at θ = 0° and those appearing as parallel components at θ = 30° are dominated by the interband mechanism. They decrease only moderately with increasing order, exhibiting a plateau-like spectral distribution. The simulated HH spectra represent reasonable agreement with the experimentally measured ones (shown as black dashed lines, plotted against the horizontal axis of harmonic order) with respect to the spectral linewidth of each harmonic order and the decreasing tendency with increasing order. From a more quantitative point of view, however, we find a discrepancy between the measured (dashed line) and simulated (shaded area) spectral profiles in each of Figs. 4(a)(b). For the parallel components shown in Fig. 4(a), the intensity contrast between even- and odd-order harmonics is less pronounced in experiments. The same trend has been experimentally observed for GaAs: even-order harmonics were much weaker than the neighboring odd harmonics in the reflection geometry and the intensity contrast between even- and odd-order harmonics was less pronounced in the transmission geometry [30]. The less-pronounced intensity contrast was ascribed to one of the propagation effects: even-order harmonics are enhanced by cascaded sum and/or difference frequency mixing processes between odd-order harmonics and fundamental pulses [30]. Considering that the propagation effect is present (absent) in our experiments (simulations), the less pronounced intensity contrast in our measured HH spectra would be attributed to the same scenario that even-order harmonics are enhanced by the cascaded nonlinear processes. For the perpendicular components shown in Fig. 4(b), we find that the fifth- and seventh-order harmonics are enhanced in the measured spectrum, which we ascribe to the cascaded sum/difference frequency mixings between the even-order harmonics and fundamental pulses.

 

Fig. 4. HH spectra for the intraband (red) and interband (blue) components, simulated with an input electric field of 1.0 V/nm for (a) parallel components for θ = 30 deg. and (b) perpendicular components for θ = 0 deg. The shaded area represents the sum of the intraband and interband components. The dashed lines represent experimentally measured HH spectra. The dash-dotted lines indicate the bandgap of GaSe.

Download Full Size | PDF

Figure 3(a) and (b) show that the orientation angle dependence is common for different orders within each of the even- and odd-orders throughout the observed frequency range. Moreover, the observed angle dependence is very similar to those reported for intraband HHG at frequencies below the bandgap of GaSe [21]. This means that the major trend of the orientation angle dependence is universal regardless of the frequency range or of the generation mechanism (although there may be differences in detailed profiles of the angle dependence). In fact, the major trend of the observed angle dependence is reasonably explained by symmetry considerations (see Appendix C), which take into account the spatial symmetries of GaSe and the time translation symmetry of the HH fields, understood as the dynamical symmetry [31]. This explanation, based solely on symmetry considerations, without specifying microscopic processes, should be applied universally to both the intraband and interband mechanisms. Note that similar correlation between the orientation dependence and the crystal symmetry has been observed also for ZnO of the 6mm point group [26]. Our experimental observation, together with the symmetry consideration, indicates that the crystal’s spatial symmetry would be clearly imprinted onto the spectral/temporal properties of HHG even at frequencies well above the bandgap, as demonstrated with GaSe at frequencies below the bandgap [20].

4. Temporal waveforms of high harmonics

We study temporal waveforms of the HH fields by numerical simulations. For simplicity, we investigate the waveforms of the parallel components generated at θ = 0°, which contain only odd-order harmonics. The calculated HH spectra for input electric fields of E = 1, 2, 3, and 4 V/nm are shown in Fig. 5(a). Similar to the case shown in Fig. 4(a), the interband components decrease more slowly with increasing harmonic order to dominate over the intraband components above the bandgap. The temporal waveforms of the harmonics simulated with varied input fields are displayed in Fig. 5 for the (b) intraband and (c) interband components. Note that the fundamental frequency component is not included in the displayed HH waveforms. Both the intraband and interband components exhibit temporal bursts that are more pronounced at higher input fields.

 

Fig. 5. Simulated HH properties of the parallel components at θ = 0 deg. for varied input electric field amplitudes of 1, 2, 3, and 4 V/nm. (a) Power spectra of HH components generated via the intraband (red) and interband (blue) mechanisms. The shaded area represents the sum of the intraband and interband components. The dash-dotted lines indicate the bandgap of GaSe. The simulated temporal waveforms of HH electric fields generated via the (b) intraband and (c) interband mechanisms. Note that the fundamental frequency component is not included in the displayed HH temporal waveforms.

Download Full Size | PDF

To see in more detail the temporal waveforms, we replot the HH waveforms shown in Fig. 5(b) and (c) into Fig. 6(a) and (b), respectively, in the fundamental sub-cycle range. For the intraband components shown in Fig. 6(a), sub-cycle spikes appear as the input field increases. Similarly, for the interband components shown in Fig. 6(b), subfemtosecond or attosecond bursts become more pronounced at higher input fields. Here, we find a striking difference in the HH temporal waveform between the intraband and interband mechanisms. The intraband components exhibit a cosine-like waveform with the amplitude maximum coinciding in time with that of the input field, whereas the interband components exhibit a sine-like waveform with its center delayed from the time of the maximum input electric field. Here, we note that simulations with different transverse relaxation times give almost identical HH waveforms for the intraband current and slightly modified waveforms for the interband polarization as presented in Appendix A. Therefore, the essence of the findings on HH waveforms discussed in the following paragraphs do not change regardless of the accuracy of the relaxation time.

 

Fig. 6. The temporal waveforms of HH components generated via the (a) intraband and (b) interband mechanisms simulated for varied input field amplitudes. The data plotted in Fig. 5(b)(c) are replotted in the fundamental sub-cycle range. The fundamental frequency component (gray shaded) is not included in the HH temporal waveforms. The characteristic time scale of the subcycle burst is denoted as $\tau $. Input electric-field dependence of the characteristic time $\tau $ of HH bursts for (c) the intraband component and (d) the interband component. In (c), the circle and square markers represent the time of the strongest HH electric field (the negative peak in Fig. 6(a)) and those of its surrounding two positive peaks, respectively. The difference between two square markers characterizes the temporal duration of the cosine-like HH burst. In (d), the circle markers and the cross markers represent the time of zero electric field at the center of the sine-like burst and its surrounding two zero-crossings of the sine-like burst, respectively. The difference between the two cross markers characterizes the temporal duration of the sine-like HH burst. The dashed lines are fitted to circles markers.

Download Full Size | PDF

We define a characteristic time $\tau $ that gives a measure of the temporal duration of a single burst. It is defined as a time period between two positive peaks of the HH field for the cosine-like burst for the intraband component and as a time period between two zero-crossings surrounding the center of the sine-like burst for the interband component. As the input field increases from 1 to 4 V/nm, $\tau $ decreases slightly from 2.15 fs to 1.86 fs for the intraband component and obviously from 1.95 fs to 0.83 fs for the interband component. Figures 6 (c)(d) display how the characteristic times of the HH bursts vary with the input field. For the intraband components shown in Fig. 6(c), the circle and square markers represent the time of the strongest HH electric field (the negative peak in Fig. 6(a)) and those of its surrounding two positive peaks, respectively. The time period between two square markers corresponds to the temporal duration $\tau $ of the cosine-like HH burst. For the interband components shown in Fig. 6(d), the circle markers and cross markers represent the time of zero electric field at the center of the sine-like burst and its surrounding two zero-crossings, respectively. The time period between the two cross markers cresponds to the temporal duration $\tau $ of the sine-like HH burst.

Here, we discuss the possible physics behind the calculated HH burst for each of the intraband and interband components. The cosine-like waveform of the intraband components is composed of odd-order harmonics, which are superposed with each other all in phase. It means that the harmonic burst is free of atto-chirp and group delay. This feature is in line with the picture of bremsstrahlung, which radiates a burst at the time of the maximum acceleration or the time of the strongest driving electric field. The sine-like waveform of the interband components is composed of odd-order harmonics, superposed with each other possessing some phase differences among them. It indicates the existence of atto-chirp and group delay. This feature is in line with the picture of interband polarization, where the photon energy radiated upon recollision of an electron hole pair depends on its excursion time in a momentum space. It can, in principle, be possible to correlate the HH waveforms with carrier dynamics in the crystal and with the resonant interband polarization responses, but further study is necessary to confirm it. Since the HH fields generated via the intraband mechanism are composed mainly by the lower-order harmonics in the investigated intensity range, the characteristic time τ decreases only moderately with the input field strength. In contrast, since the contribution from the higher-order components is noticeable and enhanced with increasing intensity for the interband mechanism (as seen from the plateau-like structures displayed in Fig. 5(a)), the characteristic time τ decreases more clearly with the input field strength. The fact that the interband components represent shorter bursts indicates the possibility of attosecond pulse generation at frequencies well above the bandgap where the interband mechanism dominates.

5. Conclusions

We studied HHG from GaSe in the deep-UV range well above the bandgap. We successfully observed HHG up to the tenth order (210 nm in wavelength, ∼3 times the bandgap energy), which is the shortest HH wavelength ever reported for GaSe HHG. The experimentally observed orientation-angle dependence of HHG, reproduced by the 3D TDDM simulations, reflects the $\bar{6}2m$ crystallographic point group, which indicates the possibilities of HH spectroscopy and HH waveform control at frequencies above the bandgap. The numerical simulations reveal that the interband polarization mechanism dominates over the intraband current mechanism above the bandgap energy. Most importantly, it is revealed that the dominating interband polarization mechanism can create attosecond bursts with characteristic sine-like waveforms. The successful demonstration of HHG in the deep-UV range and the numerical findings indicate a great promise for attosecond pulse generation through solid-state HHG as well as form the basis for the potential applications of short-wavelength light sources, HH waveform control, and HH spectroscopy at frequencies above the bandgap.

Appendix A. Impact of transverse relaxation on the harmonic spectra

In this section, we describe how we include the transverse relaxation of crystal electrons in our numerical simulations and argue for the impact of the relaxation on the harmonic spectra. We take into account only transverse relaxation because the time period of light-matter interaction, set by the input pulse duration of 100 fs, is much shorter than the longitudinal relaxation time.

We employ velocity gauge coupling to numerically solve the von Neumann equation, or semiconductor Bloch equation without the electron-hole interaction, for the one-body density matrix $\rho $.

$$i\hbar \frac{\textrm{d}}{{\textrm{d}t}}\rho = [{h(t ),\rho } ]+ s(\rho ),$$

where $h(t )$ is the one-body Hamiltonian for a crystal electron derived from density-functional theory and $s(\rho )$ is a phenomenological relaxation term. This relaxation term is introduced to imitate the electron many-body effect and interaction with different subsystems, such as atomic motion. Similar relaxation time approximations have been frequently employed in quantum dynamics simulations for crystalline solids [32].

It is important to choose a proper representation that leads to physically reasonable conclusions for the relaxation term. In the absence of any external field, the common choice is the eigenstates of the field-free Hamiltonian. In the presence of an external field, however, the representation based on a field-free Hamiltonian leads to unreasonable results: the system is excited even for an infinitely sloexternal field. Here, we employ a representation based on the instantaneous eigenfunction of the time-dependent Hamiltonian [33]. With this representation, the system remains at the ground state (of the time-dependent Hamiltonian) under an external field that changes infinitely slowly. This result is consistent with the adiabatic theorem conclusion.

Matrix elements of the density matrix, $\left\langle {{\phi _{b{\mathbf k}}}\textrm{|}\rho \left( t \right)\textrm{|}{\phi _{c{\mathbf k}}}} \right\rangle ,$ represented with the eigenfunctions of the field-free Hamiltonian, ${\phi _{b{\mathbf k}}}$ and ${\phi _{c{\mathbf k}}}$, (b and c identify the bands, and k identifies the crystal momentum) are the dynamical degrees of freedom in our theoretical framework. In simulating the relaxation, we perform a conversion from the eigenfunctions of the field-free Hamiltonian to the instantaneous eigenfunctions of the Hamiltonian with a time-varying field. This conversion is realized by the unitary transformation as

$$ s(\rho)=-\frac{i \hbar}{\tau_{R}} \sum_{\beta \gamma(\beta \neq \gamma) \boldsymbol{k}}\left|\varphi_{\beta \boldsymbol{k}}^{(t)}\right\rangle\left\langle\varphi_{\beta \boldsymbol{k}}^{(t)} \bigg| \rho(t) \bigg| \varphi_{\gamma \boldsymbol{k}}^{(t)}\right\rangle\left\langle\varphi_{\gamma \boldsymbol{k}}^{(t)}\right| $$

where ${\tau _R}$ and $\varphi _{\beta {\mathbf k}}^{(t )}$ are the relaxation time constant and an instantaneous eigenfunction of the time-dependent Hamiltonian, respectively [33]. To evaluate the matrix element $\left\langle {{\phi _{b{\mathbf k}}}\textrm{|}s\left( \rho \right)\textrm{|}{\phi _{c{\mathbf k}}}} \right\rangle$, we perform diagonalization of the Hamiltonian at each time step in the given Hilbert space spanned by the field-free basis ${\phi _{b{\mathbf k}}}$.

In Fig. 7, we show the simulated HH spectra above the bandgap of GaSe as a function of the angle $\theta $ for different relaxation times: (a-d) parallel components, (e-h) perpendicular components, (a, e) ${\tau _R} = \infty $, (b, f) $100\; \textrm{fs}$, (c, g) $10\; \textrm{fs}$, and (d, h) $1\; \textrm{fs}$. For ${\tau _R} = \infty $ and 100 fs, the simulated HH spectrum is not a series of discrete harmonic orders at any angle θ, and the angle dependence is distorted to be nonperiodic at any harmonic order. These results are hardly comparable to the experimental observations. In contrast, spectral shapes that are very similar to the experimental results, characterized by features I∼IV in Appendix C, are obtained for the simulations with ${\tau _R} = 10\; \textrm{fs}$ nd 1 fs. This fact suggests that the electrons and holes accelerated by such an intense electric field experience transverse relaxation with a time constant comparable to or shorter than 10 fs.

 

Fig. 7. Orientation angle dependence of the simulated HH spectra: (a-d) parallel and (e-h) perpendicular components simulated with different transverse relaxation times of (a, e) ${\tau _R} = \infty $, (b, f) 100 fs, (c, g) 10 fs, and (d, h) 1 fs.

Download Full Size | PDF

Figure 8 shows the HH tempol waveforms for the (a) intraband and (b) interband components, simulated with a fixed input electric field of 2 V/nm and varied transverse relaxation times of 1, 10, and 100 fs. The transverse relaxation time represents almost no impact and minor impact on temporal waveforms of (a) the intraband current and (b) interband polarization, respectively.

 

Fig. 8. The temporal waveforms of HH components simulated with a fixed input electric field of 2 V/nm and varied transverse relaxation times of 1, 10, and 100 fs. The transverse relaxation time represents almost no impact and minor impact on temporal waveforms of (a) the intraband current and (b) interband polarization, respectively.

Download Full Size | PDF

Appendix B. Input-intensity dependence

The experimental input-intensity dependence of the HH yield is shown in Fig. 9 for (a) parallel and (b) perpendicular components of each harmonic order. The HH intensity for the n-th order follows the n-th power law at an input intensity lower than 0.1 TW/cm² but deviates from it at higher input intensity. The deviation from the n-th order dependence indicates transformation from the perturbative to nonperturbative regimes. In the experiments shown in Figs. 2, 3(a) and (b), we adopt an input intensity of 0.48 TW/cm² of the nonperturbative regime.

 

Fig. 9. Input-intensity dependence of the HH yield for each harmonic order: (a) parallel components measured at θ = 30° and (b) perpendicular components measured at θ = 0° and 15° for even and odd orders, respectively.

Download Full Size | PDF

Appendix C. Symmetry consideration of the orientation angle dependence

We briefly describe how one can explain the dependence of HHG from the ɛ-type GaSe on the input field orientation.

The ɛ-type GaSe crystal has a threefold rotational symmetry (C₃) and three mirror planes (σ_v). The mirror planes, or bc-planes, include Ga-Se chemical bonds. We assume that the input field is monochromatic and linearly polarized and propagates along the c-axis of the GaSe crystal. The input field is identified by the angle $\theta $ as ${\boldsymbol E}({t;\theta } )= E(t )({\cos \theta \hat{x} + \sin \theta \hat{y}} )$, where $E(t )$ is the amplitude of the electric field at time t, $\hat{x}$ and $\hat{y}$ are unit vectors along the x and y axes, respectively, and the x(y)-axis is chosen to be parallel to the crystallographic a(b)-axis (see Fig. 1(c)). In general, the polarization response is not instantaneous, and therefore, the induced polarization at time t depends on the input electric field at the whole past time $t^{\prime} \le t$, reflecting causality. Therefore, the polarization is regarded as a functional of the time-dependent field ${\boldsymbol E}({t;\theta } )$, which we describe as ${\boldsymbol P}({t;\theta } )\equiv {\boldsymbol P}[{{\boldsymbol E}({t^{\prime};\theta } )} ](t )$ for simplicity.

We consider that an input field ${\boldsymbol E}({t;\theta } )$ induces a polarization ${\boldsymbol P}({t;\theta } )$ and then consider a case where the input field is rotated by an angle 2π/3: ${\boldsymbol E}({t;\theta + 2\pi /3} )= {R_{{C_3}}}{\boldsymbol E}({t;\theta } )$. Here, ${R_{{C_3}}}$ denotes a transformation operation for the ${C_3}$ symmetry or $2\pi /3$ rotation. Because the crystal is invariant with the ${C_3}$ symmetry operation, so is the polarization response. Therefore, the induced polarization also will be rotated by 2π/3: ${\boldsymbol P}({t;\theta + 2\pi /3} )= {R_{{C_3}}}{\boldsymbol P}({t;\theta } )$, as shown in Fig. 10(a). This invariance in the polarization response is described as

(1a)$${P_\parallel }({t;\theta } )= {P_\parallel }({t;\theta + 2\pi /3} ),$$

(1b)$${P_ \bot }({t;\theta } )= {P_ \bot }({t;\theta + 2\pi /3} ),$$

where ${P_\parallel }({t;\theta } )$ and ${P_ \bot }({t;\theta } )$ denote polarization components parallel and perpendicular to the input field ${\boldsymbol E}({t;\theta } )$, respectively, as defined by ${\boldsymbol P}({t;\theta } )= {P_\parallel }({t;\theta } ){\hat{e}_\parallel } + {P_ \bot }({t;\theta } ){\hat{e}_ \bot }$ (${\hat{e}_\parallel }$ and ${\hat{e}_ \bot }$ are unit vectors parallel and perpendicular to ${\boldsymbol E}({t;\theta } )$, respectively, with ${\hat{e}_\parallel } \times {\hat{e}_ \bot }$ oriented along $\hat{x} \times \hat{y}$).

 

Fig. 10. Schematic explanation of how the induced polarization changes upon each of (a) C₃ rotation and (b) σ_v reflection of the input field. Gray triangles represent crystals with C₃ and σ_v symmetries. Red and green arrows represent the directions of the input field and the induced polarization, respectively.

Download Full Size | PDF

By a similar consideration for the reflection operation ${\sigma _v}$ with respect to the yz-plane, we obtain ${\boldsymbol P}({t;\pi - \theta } )= {R_{{\sigma _v}}}{\boldsymbol P}({t;\theta } )$, as shown in Fig. 10(b), which is described as

(2a)$${P_\parallel }({t;\theta } )={+} {P_\parallel }({t;\pi - \theta } ),$$

(2b)$${P_ \bot }({t;\theta } )={-} {P_ \bot }({t;\pi - \theta } ).\; $$

Here, ${R_{{\sigma _v}}}$ is a transformation operation for the ${\sigma _v}$ symmetry with respect to the yz-plane.

Next, we take into account that the monochromatic input field oscillates in time with a period T, $E({t + T} )= E(t )$, and changes its sign after the half period, $E({t + T/2} )={-} E(t )$. We also assume the same periodicity for the induced polarization: ${\boldsymbol P}({t + T;\theta } )= {\boldsymbol P}({t;\theta } )$. Considering that the electric field changes its sign either by the rotation of $\pi $ or by the time shift of $T/2$, we obtain ${\boldsymbol E}({t + T/2;\theta } )= {\boldsymbol E}({t;\theta + \pi } )$ and therefore ${\boldsymbol P}({t + T/2;\theta } )= {\boldsymbol P}({t;\theta + \pi } ).\; $This invariance is described by

(3a)$${P_\parallel }({t + T/2;\theta } )={-} {P_\parallel }({t;\theta + \pi } ),$$

(3b)$${P_ \bot }({t + T/2;\theta } )={-} {P_ \bot }({t;\theta + \pi } ).$$

Here, we see that the signs of the parallel and perpendicular components are reversed although the polarization vector itself is identical. This is simply because the directions of the unit vectors ${\hat{e}_\parallel }$ and ${\hat{e}_ \bot }$ are opposite between ${\boldsymbol P}({t + T/2;\theta } )$ and ${\boldsymbol P}({t;\theta + \pi } )$. By applying Eqs. (2a, 2b) to Eqs. (3a, 3b), we obtain

(4a)$${P_\parallel }({t + T/2;\theta } )={-} {P_\parallel }({t; - \theta } ),$$

(4b)$${P_ \bot }({t + T/2;\theta } )={+} {P_ \bot }({t; - \theta } ).$$

The temporal symmetry for each of the even-order harmonics ${P^e}$ and odd-order harmonics ${P^o}$ is described as

(5a)$${P^e}({t + T/2;\theta } )= {P^e}({t;\theta } ),$$

(5b)$${P^o}({t + T/2;\theta } )={-} {P^o}({t;\theta } ).$$

By combining Eqs. (5a, 5b), (3a, 3b), and (1a, 1b), we obtain

$$P_{{\parallel} , \bot }^e({t;\theta } )={+} P_{{\parallel} , \bot }^e({t + T/2;\theta } )={-} P_{{\parallel} , \bot }^e({t;\theta + \pi } )={-} P_{{\parallel} , \bot }^e({t;\theta + \pi /3} ),$$

$$P_{{\parallel} , \bot }^o({t;\theta } )={-} P_{{\parallel} , \bot }^o({t + T/2;\theta } )={+} P_{{\parallel} , \bot }^o({t;\theta + \pi } )={+} P_{{\parallel} , \bot }^o({t;\theta + \pi /3} ).\; $$

These two equations are summarized as

(6a)$$P_\parallel ^e({t;\theta } )={-} P_\parallel ^e({t;\theta + \pi /3} ),$$

(6b)$$P_ \bot ^e({t;\theta } )={-} P_ \bot ^e({t;\theta + \pi /3} ),$$

(6c)$$P_\parallel ^o({t;\theta } )={+} P_\parallel ^o({t;\theta + \pi /3} ),$$

(6d)$$P_ \bot ^o({t;\theta } )={+} P_ \bot ^o({t;\theta + \pi /3} ).$$

Similarly, by applying Eqs. (5a, 5b) to Eqs. (4a, 4b), we obtain

(7a)$$P_\parallel ^e({t;\theta } )={-} P_\parallel ^e({t; - \theta } ),$$

(7b)$$P_ \bot ^e({t;\theta } )={+} P_ \bot ^e({t; - \theta } ),$$

(7c)$$P_\parallel ^o({t;\theta } )={+} P_\parallel ^o({t; - \theta } ),$$

(7d)$$P_ \bot ^o({t;\theta } )={-} P_ \bot ^o({t; - \theta } ).$$

Recalling that Eqs. (2a, 2b) are valid for both odd- and even-order harmonics, we obtain

(8a)$$P_\parallel ^e({t;\theta } )={+} P_\parallel ^e({t;\pi - \theta } ),$$

(8b)$$P_ \bot ^e({t,\theta } )={-} P_ \bot ^e({t;\pi - \theta } ),$$

(8c)$$P_\parallel ^o({t;\theta } )={+} P_\parallel ^o({t;\pi - \theta } ),$$

(8d)$$P_ \bot ^o({t,\theta } )={-} \; P_ \bot ^o({t;\pi - \theta } ).$$

Equations (6a–6d) indicate how each polarization component changes upon the rotation operation of the input field:

• • The magnitudes of $P_\parallel ^e$ and $P_ \bot ^e$ change with a $\pi /3$ period, and the signs of them are reversed with every $\pi /3$.
• • $P_\parallel ^o$ and $P_ \bot ^o$ change with a $\pi /3$ period $.$

Equations (7a–7d) indicate how each polarization component changes upon the reflection operation of the input field with respect to the ac-plane:

• • $P_ \bot ^e$ and $P_\parallel ^o$ change in a symmetric manner.
• • $P_\parallel ^e$ and $P_ \bot ^o$ change in an anti-symmetric manner and therefore take zero at $\theta = 0$.

Equations (8a–8d) indicate how each polarization component changes upon the reflection operation of the input field with respect to the bc-plane:

• • $P_\parallel ^e$ and $P_\parallel ^o$ change in a symmetric manner.
• • $P_ \bot ^e$ and $P_ \bot ^o$ change in an anti-symmetric manner and therefore take zero at $\theta = \pi /2$.

These properties lead to the following features of the squared values.

• I. All of ${|{P_\parallel^e} |^2},\; {|{P_ \bot^e} |^2},\; {|{P_\parallel^o} |^2}$, and ${|{P_ \bot^o} |^2}$ represent a $\pi /3$ periodicity and symmetry with respect to $\theta = n\pi /3$ and $\pi /6 + n\pi /3$.
• II. ${|{P_\parallel^e} |^2}$ represents nodes at $\theta = n\pi /3$.
• III. ${|{P_ \bot^e} |^2}$ represents nodes at $\theta = \pi /6 + n\pi /3$.
• IV. ${|{P_ \bot^o} |^2}$ represents nodes at $\theta = n\pi /3$ and $\pi /6 + n\pi /3$.

Note that the measured HH intensity for each of the parallel and perpendicular components is proportional to ${|{P_\parallel^{o,e}} |^2}$ and ${|{P_ \bot^{o,e}} |^2}$, respectively. Indeed, the experimental/numerical results of the orientation angle dependence presented in Fig. 3 agree with features I-IV. Since we assume only crystal symmetry and a linearly polarized monochromatic input field, the abovementioned consequence does not reflect detailed information on the band structure. However, the consequence is generally valid regardless of whether the system is in the perturbative or nonperturbative regime or whether the dominant mechanism is the intraband current or interband polarization.

Funding

Japan Science and Technology Agency (Center of Innovation Program JPMJCE1313, CREST JPMJCR15N1, Research and Education Consortium); Japan Society for the Promotion of Science (JP18K14145, JP18K19030, JP19H02623, JP20H02651, JP20H05670, JP20K20556); Ministry of Education, Culture, Sports, Science and Technology (Exploratory Challenge on Post-K Computer, Q-LEAP JPMXS0118067246, Q-LEAP JPMXS0118068681).

Acknowledgments

The computation in this work was done using the facilities of the Supercomputer Center (the Institute for Solid State Physics, the University of Tokyo) and also using the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp180067). K. I. and S. A. are grateful to T. Shimura at the University of Tokyo for his constructive advice.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, “Observation of high-order harmonic generation in a bulk crystal,” Nat. Phys. 7(2), 138–141 (2011). [CrossRef]  

2. S. Ghimire and D. A. Reis, “High-harmonic generation from solids,” Nat. Phys. 15(1), 10–16 (2019). [CrossRef]  

3. S. Gholam-Mirzaei, J. Beetar, and M. Chini, “High harmonic generation in ZnO with a high-power mid-IR OPA,” Appl. Phys. Lett. 110(6), 061101 (2017). [CrossRef]  

4. T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, “Extreme ultraviolet high-harmonic spectroscopy of solids,” Nature 521(7553), 498–502 (2015). [CrossRef]  

5. G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, D. D. Klug, and P. B. Corkum, “All-optical reconstruction of crystal band structure,” Phys. Rev. Lett. 115(19), 193603 (2015). [CrossRef]  

6. M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Klostermann, A. Guggenmos, and E. Goulielmakis, “Multi-petahertz electronic metrology,” Nature 538(7625), 359–363 (2016). [CrossRef]  

7. Y. S. You, D. A. Reis, and S. Ghimire, “Anisotropic high-harmonic generation in bulk crystals,” Nat. Phys. 13(4), 345–349 (2017). [CrossRef]  

8. G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, “Linking high harmonics from gases and solids,” Nature 522(7557), 462–464 (2015). [CrossRef]  

9. S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, “High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure,” Nat. Commun. 7(1), 13105 (2016). [CrossRef]  

10. M. Sivis, M. Taucer, G. Vampa, K. Johnston, A. Staudte, A. Y. Naumov, D. M. Villeneuve, C. Ropers, and P. B. Corkum, “Tailored semiconductors for high-harmonic optoelectronics,” Science 357(6348), 303–306 (2017). [CrossRef]  

11. G. Vampa, B. G. Ghamsari, S. S. Mousavi, T. J. Hammond, A. Olivieri, E. Lisicka-Skrek, A. Yu Naumov, D. M. Villeneuve, A. Staudte, P. Berini, and P. B. Corkum, “Plasmon-enhanced high-harmonic generation from silicon,” Nat. Phys. 13(7), 659–662 (2017). [CrossRef]  

12. K. Imasaka, T. Kaji, T. Shimura, and S. Ashihara, “Antenna-enhanced high harmonic generation in a wide-bandgap semiconductor ZnO,” Opt. Express 26(16), 21364–21374 (2018). [CrossRef]  

13. G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, “Theoretical analysis of high-harmonic generation in solids,” Phys. Rev. Lett. 113(7), 073901 (2014). [CrossRef]  

14. T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto, M. Kuwata-Gonokami, and K. L. Ishikawa, “Trajectory analysis of high-order-harmonic generation from periodic crystals,” Phys. Rev. A 95(4), 043416 (2017). [CrossRef]  

15. M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B. Gaarde, “High-harmonic generation from Bloch electrons in solids,” Phys. Rev. A 91(4), 043839 (2015). [CrossRef]  

16. A. A. Lanin, E. A. Stepanov, A. B. Fedotov, and A. M. Zheltikov, “Mapping the electron band structure by intraband high-harmonic generation in solids,” Optica 4(5), 516–519 (2017). [CrossRef]  

17. Y. S. You, M. Wu, Y. Yin, A. Chew, X. Ren, S. Gholam-Mirzaei, D. A. Browne, M. Chini, Z. Chang, K. J. Schafer, M. B. Gaarde, and S. Ghimire, “Laser waveform control of extreme ultraviolet high harmonics from solids,” Opt. Lett. 42(9), 1816–1819 (2017). [CrossRef]  

18. M. Garg, H. Y. Kim, and E. Goulielmakis, “Ultimate waveform reproducibility of extreme-ultraviolet pulses by high-harmonic generation in quartz,” Nat. Photonics 12(5), 291–296 (2018). [CrossRef]  

19. O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, “Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations,” Nat. Photonics 8(2), 119–123 (2014). [CrossRef]  

20. F. Langer, M. Hohenleutner, U. Huttner, S. W. Koch, M. Kira, and R. Huber, “Symmetry-controlled temporal structure of high-harmonic carrier fields from a bulk crystal,” Nat. Photonics 11(4), 227–231 (2017). [CrossRef]  

21. K. Kaneshima, Y. Shinohara, K. Takeuchi, N. Ishii, K. Imasaka, T. Kaji, S. Ashihara, K. L. Ishikawa, and J. Itatani, “Polarization-resolved study of high harmonics from bulk semiconductors,” Phys. Rev. Lett. 120(24), 243903 (2018). [CrossRef]  

22. Y. Morimoto, Y. Shinohara, M. Tani, B.-H. Chen, K. L. Ishikawa, and P. Baum, “Asymmetric single-cycle control of valence electron motion in polar chemical bonds,” Optica 8(3), 382–387 (2021). [CrossRef]  

23. N. C. Fernelius, “Properties of gallium selenide single crystal,” Prog. Cryst. Growth Charact. Mater. 28(4), 275–353 (1994). [CrossRef]  

24. https://elk.sourceforge.io/.

25. J. P. Perdew and Y. Wang, “Accurate and simple analytic representation of the electron-gas correlation energy,” Phys. Rev. B 45(23), 13244–13249 (1992). [CrossRef]  

26. S. Jiang, S. Gholam-Mirzaei, E. Crites, J. E. Beetar, M. Singh, R. Lu, M. Chini, and C. D. Lin, “Crystal symmetry and polarization of high-order harmonics in ZnO,” J. Phys. B: At. Mol. Opt. Phys. 52(22), 225601 (2019). [CrossRef]  

27. S. Ghimire, A. D. DiChiara, E. Sistrunk, G. Ndabashimiye, U. B. Szafruga, A. Mohammad, P. Agostini, L. F. DiMauro, and D. A. Reis, “Generation and propagation of high-order harmonics in crystals,” Phys. Rev. A 85(4), 043836 (2012). [CrossRef]  

28. G. Vampa, Y. S. You, H. Liu, S. Ghimire, and D. A. Reis, “Observation of backward high-harmonic emission from solids,” Opt. Express 26(9), 12210–12218 (2018). [CrossRef]  

29. R. W. Boyd, Nonlinear Optics,3rd edition (Academic Press, 2008).

30. P. Xia, C. Kim, F. Lu, T. Kanai, H. Akiyama, J. Itatani, and N. Ishii, “Nonlinear propagation effects in high harmonic generation in reflection and transmission from gallium arsenide,” Opt. Express 26(22), 29393–29400 (2018). [CrossRef]  

31. O. Neufeld, D. Podolsky, and O. Cohen, “Floquet group theory and its application to selection rules in harmonic generation,” Nat. Commun. 10(1), 405 (2019). [CrossRef]  

32. M. Lindberg and S. W. Koch, “Effective Bloch equations for semiconductors,” Phys. Rev. B 38(5), 3342–3350 (1988). [CrossRef]  

33. Y. Sanari, H. Hirori, T. Aharen, H. Tahara, Y. Shinohara, K. L. Ishikawa, T. Otobe, P. Xia, N. Ishii, J. Itatani, S. A. Sato, and Y. Kanemitsu, “Role of virtual band population for high harmonic generation in solids,” Phys. Rev. B 102(4), 041125 (2020). [CrossRef]