High-harmonic generation (HHG) is one of the most intriguing topics among the phenomena arising from the interaction of matter with intense laser fields. HHG originating from the ultrafast electron dynamics of atoms or molecules can be used to study the field-driven physics between strong laser pulse and matter. It has irreplaceable applications in many research fields, such as coherent extreme ultraviolet (EUV) to X-ray^[1–3], chemical reaction tracing^[4], molecular orbital imaging^[5–7], and attosecond science^[8–12]. In recent years, the nonperturbative HHG has also been observed in solid materials^[13,14]. It has been proved to be a useful tool to reconstruct the electronic band structures and reveal the ultrafast electron dynamics in condensed matter all-optically^[13–19].

HHG from solids is fundamentally different from that in the gaseous medium, owing to a variety of electronic structures in the lattice of bulk crystals^{[13–17,20–29]}. Both the intraband current and interband polarization of the optically excited carriers can contribute to the emission of high-harmonics, despite the two mechanisms often intertwining^{[20–23,30,31]}. In 2015, Vampa et al. found that the generalized recollision between the field-driven electron-hole pairs is the primary mechanism of HHG in ZnO^[24], which brings new ideas for the study of solid-state HHG. In their study, two-color laser fields composed of the fundamental and second-harmonic generation (SHG) were adopted as a pump. The weak SHG can unbalance the electrons’ trajectories in adjacent half-cycles of driving field, and thus lead to the emission of even harmonics. Latterly, the modulation of high-harmonic emission was found to be sensitive to the symmetry properties of solid materials, which provides a method to analyze the mechanism of HHG^[25–27]. Recently, two-color laser fields driving HHG have been adopted to disentangle the electron-hole dynamics from the surface charge region and the crystal bulk^[27]. This has been used to establish the all-optical spectroscopy of a strongly driven crystal and reveal a laser-induced modification of the band structure^[32].

Three-dimensional (3D) topological insulators (TIs) feature unconventional two-dimensional topological surface states (TSSs). In TSSs, the carriers are massless Dirac fermions, which provide a dissipationless transport channel on the TIs’ surface^[33,34]. Very recently, the observations of HHG from 3D TIs have been reported by several groups, and experimental evidence of HHG originating from strong field dynamics of TSSs has been demonstrated^[35–37]. Despite the exciting progress, strong laser fields can inevitably excite transitions between the bulk states and surface states and, in consequence, affect the nonlinear response. Thus, how to disentangle the bulk and surface contributions of HHG is still an open question.

In this work, we investigated the HHG from ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ crystal, a prototypical 3D TI, in response to two-color laser fields. The modulation of even harmonics sensitively depends on the retarded phase, $\varphi$, between the fundamental and SHG fields. By adjusting the SHG intensity, two periodic modulations were observed. We found that the harmonic yield originating from the surface states modulates once per cycle of the SHG, while the harmonic yield from the bulk states modulates twice per cycle. The observations are well explained by the generalized recollision model based on the semiconductor Bloch equations^[24,38,39].

The experimental setup is shown in Fig. 1. The fundamental pulse was produced at the wavelength of 3820 nm from a home-built system that is composed of two identical three-stage optical parametric amplifiers (OPAs) and a difference frequency generation (DFG) module^[40]. Subsequently, the fundamental pulse passed through an ${\mathrm{AgGaS}}_2$ (AGS) crystal with the thickness of 400 µm to produce SHG at the wavelength of 1910 nm. Both the fundamental and SHG pulses were synchronized and collinearly focused onto a ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ crystal. A zinc sulfide (ZnS) plate was used to compensate for the group delay between the fundamental and SHG pulses. The dual-wavelength wave plate (DWP), a half-wave plate for the fundamental and a full-wave plate for the SHG, was adopted to adjust the polarization of the two-color laser fields. By tuning the angle of the DWP and a polarizer (WP25M-UB, Thorlabs), the SHG intensity was maintained much more weakly than the fundamental intensity ($<2\times {10}^{-2}$) as a perturbation factor. Phase delay of the two pulses was tuned by a pair of ${\mathrm{BaF}}_2$ wedges with an angle of 3°. In the measurements, the phase delay was adjusted by a step of $0.1\pi$. The ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ crystal was prepared by using the self-flux method described in Ref. [35]. The emitted high-harmonic spectra were measured by a grating spectrometer equipped with a CCD camera (Andor Shamrock 303i, Newton 920).

Figure 2 shows a typical HHG spectrum pumped by the two-color laser fields; as a comparison, the HHG spectrum for fundamental pulse only is also shown. In measurements, the SHG intensity is $3.2\times {10}^{-4}$ times of the fundamental intensity, which is not enough to independently produce HHG. Driven by the fundamental pulse alone, even harmonics should stem from inversion asymmetric states confined in the surface^[35,36]. As the weak SHG is added, notable changes in even harmonics are observed with the increasing yield of the fourth and eighth harmonic (H4, H8), while the sixth harmonic (H6) is opposite. Compared with the even harmonics, the odd harmonics’ intensity remains unchanged. This indicates that even harmonics are sensitive to the weak SHG field, which may originate from the surface states that lack inversion symmetry, while the odd harmonics may come from bulk states with inversion and time-reversal symmetry ^[24,27]. The variation of even harmonics can be interpreted qualitatively by the interference effects of HHG produced from asymmetric two-color laser fields^[41].

By changing the phase delay of the two-color fields, periodic modulations of HHG were observed experimentally, as shown in Figs. 3(a)–3(c). For weak SHG intensity of $0.025{\mathrm{GW}/\mathrm{cm}}^2$, the yield of H4 and H6 oscillates once per cycle of the SHG, that is, the modulation period equals ${\mathit{T}}_{2\omega }$. However, the modulation of H5 is absent in the same measurement. As the SHG intensity is increased to $1.2{\mathrm{GW}/\mathrm{cm}}^2$ ($1.6\times {10}^{-2}$ times that of the fundamental), the modulation of H5 appears. The intensity of H5 modulates twice per cycle of the SHG, that is, the modulation period is ${\mathit{T}}_{2\omega }/2$. The modulation depth of H6 is enhanced; however, the modulation period does not change. At the same time, the intensity modulation of H4 shows a distinctive feature: two extrema emerge in one SHG cycle. By Fourier transform, two kinds of modulation periods can be extracted: one is ${\mathit{T}}_{2\omega }$; the other is ${\mathit{T}}_{2\omega }/2$. Therefore, we speculate that for very weak SHG, even harmonics with ${\mathit{T}}_{2\omega }$ modulation originate from asymmetric surface states. With the increase of the SHG intensity, the contribution from the bulk states of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ begins to play a role. Therefore, the multiperiodic modulation of even harmonics can be attributed to the interplay between the surface states and bulk states.

In order to understand why there are two kinds of modulations for the even harmonics, we discussed the HHG mechanisms by using semiconductor Bloch equations (SBEs)^[28]. It is known that the modulations of HHG are closely related to the symmetry of solids^[25-27]. According to Ref. [28], the complex transition dipole moment (TDM) ${\mathit{D}}_{cv}(\mathit{k})$ that describes the interband coherence contributes to even harmonics from crystals that lack inversion symmetry^[28]. For the 3D TI, the symmetry is different for the electronic states inside the bulk and confined at the surface^[34,42]. For the bulk states with inversion symmetry, TDM follows the restrictions of $D_{||,cv}(-k)=D_{||,cv}^{*}(k)$ and $D_{||,cv}(-k)=\pm D_{||,cv}(k)$, which should be purely real or imaginary. For the surface states, however, the TDM that lacks inversion symmetry should take a complex value, $D_{||,cv}(k)=|D_{||,cv}(k)|e^{i\theta (k)}$; here $\theta (k)$ is the nontrivial transition dipole phase (TDP)^[28].

To find out how TDM affects the observed HHG modulation in our experiments, we focus on the interband polarization between conduction and valence bands in SBEs. The interband current generated by a single quantum orbital satisfies^[28]$${\mathit{J}}_{\text{inter}}(t)\propto D_{||,cv}^{*}(k)D_{||,cv}(k^{\prime })F(t^{\prime })e^{-i\mathrm{\Omega }t-iS(k,t,t^{\prime })-\frac{t-t^{\prime }}{T_2}},$$(1)where $k^{\prime }=k+A(t)-A(t^{\prime })$ is the crystal momentum of electrons excited from valence band to conduction band, $k$ is the crystal momentum of electrons recombining back to the valence band, $\mathit{A}(t)$ is the vector potential of the laser field $F(t)$, $S(\mathit{k},t,t^{\prime })$ is the classical action, $\mathrm{\Omega }=N\omega$ is the frequency of the $N$th harmonic with $\omega$ the fundamental frequency, and ${\mathit{T}}_2$ is the dephasing time. In generalized recollision, the electron-hole pair is created at a time of $t^{\prime }$ and accelerated in the laser field and then recollides with each other at a time of $t$, thus, bringing the emission of high harmonics.

Considering that the SHG is very weak and only perturbs the dynamical phase of the electron wave packet, the classical action is given by^[38]$$S_2(t,\varphi )=S_1(t)-{\int }_{t^{\prime }}^t\mathrm{d}\tau v(\tau ,t^{\prime })·{\mathit{A}}_2(\tau ,\varphi ),$$(2)where ${\mathit{A}}_2(\tau ,\varphi )=A_2\sin (2\omega \tau +\varphi )$ is the vector potential of the SHG, $S_1$ the action related to the fundamental pulse only, $\mathit{v}(\mathit{k})={\nabla }_{\mathit{k}}{\epsilon }_g(\mathit{k})$ the group velocity, and ${\epsilon }_g(\mathit{k})$ the momentum dependent bandgap. The second term in Eq. (2) is an additional dipole phase induced by the SHG, $\sigma (t,\varphi )={\int }_{t^{\prime }}^t\mathrm{d}\tau v(\tau ,t^{\prime })·{\mathit{A}}_2(\tau ,\varphi )$, where $t^{\prime }$ and $t$ are the excitation and the recombination time of electron-hole pairs, respectively.

We calculated the interband current by summing up the transient current from adjacent half-cycles, which is similar to the derivation in Ref. [28], $$\begin{gathered} {\mathit{J}}_{||,\text{inter}}(t)\propto [e^{i\sigma (t,\varphi )}-e^{-iN\pi }·e^{-i\sigma (t,\varphi )}] \\ ·D_{||,cv}^{*}(k)D_{||,cv}(k^{\prime })F(t^{\prime })e^{-iN\omega t-iS(k,t,t^{\prime })-\frac{t-t^{\prime }}{T_2}}, \end{gathered}$$(3)where the additional dipole phase is the origin of even harmonics, and can be written in the form of $\sigma (t,\varphi )={\sigma }_a\cos \varphi +{\sigma }_b\sin \varphi =\mathrm{\Phi }(t)\cos [\varphi -\alpha (t)]$^[24,38], where $\mathrm{\Phi }(t)=\sqrt{{\sigma }_a^2+{\sigma }_b^2}$ and $\alpha (t)=a\tan ({\sigma }_b/{\sigma }_a)$. Clearly, the additional dipole phase oscillates with the relative phase delay, $\varphi$. The amplitude of the $N$th harmonic is proportional to ${\mathrm{\Delta }}_N=[e^{i\sigma (t,\varphi )}-e^{-iN\pi }·e^{-i\sigma (t,\varphi )}]D_{||,cv}^{*}(k_s)D_{||,cv}(k_s^{\prime })$. Then, the high-harmonic yield can be calculated by $S_{\mathrm{HHG}}(\mathrm{\Omega })\propto {|{\mathrm{\Delta }}_N|}^2$. For HHG from bulk states of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ with real-valued TDM, the spectral intensity of odd harmonics is $S_{\mathrm{odd}}(\mathrm{\Omega })\propto {\cos }^2\sigma$, and that of even harmonics is $S_{\text{even}}(\mathrm{\Omega })\propto {\sin }^2\sigma$. The additional dipole is far less than 1, due to the very weak SHG field. Thus, the spectral intensity of odd harmonics can be simplified as $$S_{\mathrm{odd}}\propto 1-{\mathrm{\Phi }}^2(t){\cos }^2[\varphi -\alpha (t)].$$(4)The spectral intensity of even harmonics is $$S_{\text{even}}\propto {\mathrm{\Phi }}^2(t){\cos }^2[\varphi -\alpha (t)].$$(5)Therefore, the modulation periods of both even and odd harmonics from bulk states are all $T_{2\omega }/2$.

For the surface states, the dipole phase plays an indispensable role. Because the SHG is very weak, its influence on TDM can be ignored^[24,28,38]. Taking the expression of complex TDM, we can get ${\mathrm{\Delta }}_N=\{e^{i[\mathrm{\Delta }\theta +\sigma (t,\varphi )]}-e^{-iN\pi }·e^{-i[\mathrm{\Delta }\theta +\sigma (t,\varphi )]}\}·|D_{||,cv}(k)|·|D_{||,cv}(k^{\prime })|$, where $\mathrm{\Delta }\theta =\theta (k^{\prime })-\theta (k)$. The spectral intensity of the odd harmonics is $S_{\mathrm{odd}}(\mathrm{\Omega })\propto {\cos }^2(\mathrm{\Delta }\theta +\sigma )$, and that of the even harmonics is $S_{\text{even}}(\mathrm{\Omega })\propto {\sin }^2(\mathrm{\Delta }\theta +\sigma )$. Considering $\sigma$ is very small, we can get $$S_{\mathrm{odd}}\propto {\cos }^2\mathrm{\Delta }\theta -\sin (2\mathrm{\Delta }\theta )·\sigma (t,\varphi )-\cos (2\mathrm{\Delta }\theta )·{\sigma }^2(t,\varphi ),$$(6)and $$S_{\text{even}}\propto {\sin }^2\mathrm{\Delta }\theta +\sin (2\mathrm{\Delta }\theta )·\sigma (t,\varphi )+\cos (2\mathrm{\Delta }\theta )·{\sigma }^2(t,\varphi ).$$(7)

In both Eqs. (6) and (7), the first terms are contributed from the nontrivial dipole phase, and the second and third terms describe the modulation due to the SHG field that oscillates as a function of phase delay $\varphi$. The modulation periods of HHG from the surface states are all ${\mathit{T}}_{2\omega }$. Similar modulation has been observed in recent experiments; however, no in-depth explanation has been given^[26,27]. Our theoretical analysis further explains the relationship between the modulation of even harmonics and the symmetry restrictions. We summarize these symmetry restrictions of even harmonics in Table 1.

The analytical expression of the intraband current is relatively complex; however, it has a form similar to that of the interband polarization^[28]. Through the above analytical analysis, we interpret the origin of the two kinds of HHG modulation.

Our analysis above focuses on the contributions of surface states and bulk states to HHG. It confirms our speculation that for very weak SHG, even harmonics with ${\mathit{T}}_{2\omega }$ modulation originate from asymmetric surface states, and with the increase of the SHG intensity, the contribution from the bulk states of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ begins to play a role. It is known that the contribution from surface states of TIs contains TSSs and two-dimensional electron gas (2DEG) in the process of HHG^[35]. Here we carried out the simulation by solving the SBEs with the Hamiltonian of TSSs and bulk states derived from the tight-binding model (TBM)^[43]. We calculated the HHG from TSSs and bulk states, respectively. The grid points of two-dimensional momentum grid are set to $400\times 400$ in the length gauge. The wavelength of the fundamental is set to 3820 nm, the intensity is set to $40\mathrm{GW}/{\mathrm{cm}}^2$, and the wavelength of the SHG is set to 1910 nm. The dephasing times are set at ${\mathit{T}}_1=200\mathrm{fs}$ and ${\mathit{T}}_2=5\mathrm{fs}$.

Figure 4(a) shows the high-harmonic spectrum from TSSs with and without the SHG in pump pulses, respectively. The intensity of H4 increases, and the intensity of H6 decreases in the presence of the SHG compared with the absence of the SHG, which is consistent with the experimental observations displayed in Fig. 2. Therefore, we speculate that the even harmonics in the experiment mainly come from TSSs. For comparison, in Fig. 4(b) we plot the high-harmonic spectrum from bulk states driven by the fundamental and two-color laser fields. Unsurprisingly, even harmonics vanish when the SHG is absent, whereas even harmonics can be produced by adding very weak SHG. These indicate that the weak SHG breaks the symmetry of the driven field and induces asymmetry in electron-hole trajectory; thus, the even harmonics are produced^[23].

We also calculated the high-harmonic yield that modulates relative to the phase delay $\varphi$, as shown in Figs. 5(a)–5(d). For HHG from the TSSs, the modulation periods of H4 and H5 are all equal to ${\mathit{T}}_{2\omega }$. The modulation depth of harmonic yield scales linearly with the increase of the SHG intensity [Figs. 5(a) and 5(b)]. Figures 5(c) and 5(d) show the modulation of HHG from the bulk states; however, the modulation periods of H4 are ${\mathit{T}}_{2\omega }/2$ and the modulation depth is proportional to the SHG intensity. Besides, the modulation of H5 begins to emerge as the SHG intensity reaches $40{\mathrm{MW}/\mathrm{cm}}^2$ [Fig. 5(d)]. As the SHG intensity further increases, the total yield of H5 decreases, which is qualitatively consistent with the experimental measurements in Fig. 3(b). Therefore, the HHG modulation in the numerical calculation verifies the analytical analysis above.

To conclude, the experimentally observed even harmonics with the modulation period of ${\mathit{T}}_{2\omega }$ can be attributed to the coupling between the nontrivial dipole phase of the surface states and the additional dynamical phase accumulated in the acceleration of the electron-hole pairs. The ${\mathit{T}}_{2\omega }/2$ modulation emerges in even harmonics originating from the bulk states as the SHG intensity further increases. Our findings are helpful in distinguishing the strong field-driven nonperturbative dynamics inside topologically protected surface states and bulk states. High-harmonic emission provides a novel ultrafast spectroscopy for probing the strong-field physics in a topological matter. In addition, it also gives new inspiration to manipulate the transportation of helical Dirac fermions by using a weak perturbing laser field. At the same time, we notice that the surface states of TIs contain rich physical properties such as TSSs and 2DEG, so how to distinguish their contribution to HHG is still a problem to be explored^[35,42].