Introduction

Nowadays, due to their unique transport properties, topological insulators (TIs) attract great attention1,2. Due to the symmetry-protected Dirac fermions with nontrivial topology on their surface, they are highly prospective as functional materials in future electronic and spintronic devices3. Moreover, various optoelectronic and thermoelectric applications can benefit from the combination of protected charge transport in surface states, and large Seebeck coefficients. For many of these applications it is crucial to understand the characteristic timescales of the relaxation dynamics of excited carriers, and in particular determine if these dynamics are different for topological surface states compared to the bulk.

The carrier dynamics in TIs have been addressed using various pump–probe techniques4,5,6,7,8,9,10,11,12. However, despite these experimental efforts, the exact mechanism and timescale of carrier relaxation of the Dirac fermions in surface states of TIs are still under debate. The reason for this is that it is challenging to experimentally disentangle excitation and relaxation channels of surface and bulk states. This is caused by the relatively small bandgap of TIs of the bismuth and antimony chalcogenide family, which is a few hundred meV13, meaning that optical excitation with near-infrared and visible light leads to interband transitions between the bulk bands. Thus, in order to optically address the surface states, energetically located within the bandgap, photons in the mid-infrared or terahertz regimes are required. Furthermore, the Fermi energy is typically located either in the valence or the conduction band for binary TI compounds, e.g. Bi2Se3 and Bi2Te3, such that bulk states are already populated. Whereas separating surface dynamics from bulk dynamics has recently been achieved at the cryogenic temperature of 5 K11, this is not the case at more technologically relevant temperatures. The main reason for this is the occurrence of phonon-assisted surface-to-bulk scattering above the Debye temperature (~180 K for Bi2Se3)7.

Here, we overcome these challenges by combining optical excitation with low photon energies and a TI sample with the Fermi energy located inside the bandgap. As a result, we are able to isolate the response of Dirac electrons in the surface states, without the contribution of bulk states. In particular, we use THz pulses with photon energies below 4 meV, and verify whether the observed carrier dynamics originate from surface states (SSs), bulk conduction states (BCSs), or bulk valence states (BVSs), using three different TI samples with the Fermi energy in the valence band (Bi2Te3), in the conduction band (Bi2Se3), and inside the bandgap (Bi1.4Sb0.6Te1.51Se1.49, i.e. BSTS). THz excitation of bulk-conducting Bi2Se3 and Bi2Te3 corresponds to intraband excitation of carriers in the bulk and at the surface simultaneously, such that surface and bulk relaxation dynamics cannot be disentangled in a straightforward manner. However, the BCSs of the BSTS are empty while the BVSs are completely filled, so that intraband transitions in the bulk of BSTS are Pauli blocked. As a result, in the BSTS sample THz excitation is dominated by the intraband excitation of surface states only. Using two THz-based nonlinear techniques, namely THz-pump optical-probe (TPOP) and THz high-harmonic generation (THz HHG), we find that relaxation of Dirac fermions at the surface of TIs occurs on a timescale of a few hundred femtoseconds, which is significantly faster than the relaxation of bulk carriers.

Results and discussion

Samples

Our three different samples, rhombohedral Bi2Se3, Bi2Te3, and Bi1.4Sb0.6Te1.51Se1.49 (BSTS) films were grown on (0001) Al2O3 substrates. Their film thicknesses are 30 nm, 490 nm, and 375 nm, respectively. The growth method is described in detail in Ref. 14. It was shown in15 that for topological insulators of the Bi2−xSbxTe3−ySey quaternary compounds, there is an optimal curve (Ren’s curve) in the composition diagram y(x). Near this curve, the properties of electronic edge states are most pronounced due to the fact that acceptors and donors mutually cancel each other. Compositions below Ren’s curve, as a rule of thumb, have predominant hole conductivity, while a composition above Ren’s curve exhibits electronic conductivity. Using THz transmission measurements, we estimate the Fermi level of Bi2Se3 to be located 26 meV above the bottom of the conduction band, and for Bi2Te3 we find it 30 meV below the top of the valence band. In case of BSTS, we find that the Fermi level is located inside the bandgap, being 110 meV below the conduction band, which is in agreement with earlier works13,15,16,17 (see Supplementary Fig. 13).

Terahertz-pump optical-probe measurements

We first measure the relaxation dynamics of these samples using transient reflectivity spectroscopy. We use laser pulses with 800 nm central wavelength and 100 fs pulse duration for probing, while single cycle THz pulses serve as a pump in THz-pump optical-probe (TPOP) experiments. As a reference, we perform optical-pump optical-probe (OPOP) measurements, where 800 nm pulses also serve as pump, instead of THz light. The experimental setups used for TPOP and OPOP measurements are described in the Supplementary information (Supplementary Fig. 1). Ultrafast pump–probe reflectivity measurements in the near- and mid-infrared ranges have been used earlier to study the dynamics of charge carriers4,8,9,10,18 and phonons19,20,21 in topological insulators. The optical reflectivity in TIs at near-IR frequencies is determined by interband transitions between the valence and conduction bands, and depends on the number of unoccupied states in the conduction band10. The sign of the optical reflectivity can be both negative and positive, where the former (latter) is mostly associated with carrier temperature (carrier density), see4,8,9,10. The negative reflectivity change is most likely the effect of Pauli blocking of interband excitations at the probe energy due to carrier heating, i.e. broadening of the Fermi–Dirac distributions of the different transitions that play a role at our probe energy. The probe-energy and probe-fluence dependent reflectivity of BSTS was addressed in detail in10, where a continuum consisting of the responses of transitions between different valence and conduction bands was shown to play a role. Similar effects likely play a role for the other TIs.

In previous experiments4,8,9,10, however, the mid- or near-IR pump pulses had photon energies above or close to the energy of the bandgap, such that the dynamics of bulk and surface carriers both contribute to the observed signals. The description of excited carrier dynamics in TIs is further complicated by the formation of metastable states at the conduction band bottom5,6. The slow population decay in these states can compete with the faster carrier cooling. In a few pump–probe studies, THz light was used as the pump, in order to study carriers11, and phonons21,22. However, in both cases bulk-metallic TIs were used. Therefore, the carrier relaxation dynamics observed in all of those studies always contained contributions from both the surface and bulk electronic systems.

In Fig. 1 we show normalized optical-probe reflectivity changes for the three topological insulators investigated in our work after excitation by broadband THz radiation with 400 kV cm−1 peak field strength. We point out that the sign of the observed reflectivity variations for all three cases is negative, in agreement with8,10. The penetration depth of the 1.5 eV probe pulse is few tens of nm10, thus the optical reflectivity signal does not carry information on the dynamics at the sample-to-substrate interface. The observed dynamics are significantly different for each of the three samples, and can be decomposed into several components: ultrafast build-up, initial decay and a Raman-excited coherent phonon response. The fast rise time of the THz response corresponds to the ultrafast thermalization of excited carriers through electron-electron scattering. As measured with optical techniques and time- and angle-resolved photoemission spectroscopy (tr-ARPES), the electron-electron scattering time constant has been determined to be lower than 200 fs for Bi2Se37. The relaxation dynamics of Bi2Se3 (Fig. 1a, TPOP) can be decomposed into two different exponential decays. The faster decay rate varies from 1.45 to 1.9 ps under excitation with 280 and 400 kV cm−1 fields (see Supplementary Fig. 8) and can be attributed to carrier cooling via electron–phonon scattering5,7,8,11. The slow decay component shows a time constant of about 190 ps (see Supplementary Fig. 7). The presence of a slow relaxation upon sub-bandgap excitation in Bi2Se3 was observed previously using tr-ARPES23 and was attributed to thermal excitation of the carrier population. As mentioned before, due to the coexistence of partially filled BCS and gapless Dirac SS24, it is not possible to separate relaxation processes occurring in bulk and surface states based on these measurements.

Fig. 1: Transient reflectivity spectroscopy.
figure 1

THz-pump optical-probe (TPOP) reflectivity measurements in Bi2Se3 (a), Bi2Te3 (b) and BSTS (c) (blue line). The intensity transient of the driving THz field is shown in all panels (gray shaded line). The red line corresponds to the dynamics when using optical excitation with 1.5 eV photon energy (OPOP). The band diagram with the Fermi level of each sample, obtained from THz transmission experiments, is schematically shown in the insets.

Full size image

Measurements performed on the Bi2Te3 sample are presented in Fig. 1b. The Raman-excited optical phonon25, seen as periodic oscillation, is more pronounced than in Bi2Se321 and slightly red-shifted due to the heavier Te atoms (Supplementary Fig. 10). After a sharp increase in reflectivity variation, the non-oscillatory part decreases within 2 ps via electron–phonon scattering and then monotonously increases again reaching a plateau at a pump–probe delay of 15 ps. At longer delays, the reflectivity variation slowly decreases with a decay timescale of about 600 ps (shown in Supplementary Fig. 7). Such dynamics are very similar to those observed for 3 eV pump excitation25 that was attributed to the interplay between free carrier absorption and the band filling effect.

In comparison with Bi2Te3 and Bi2Se3, the carrier dynamics in BSTS (Fig. 1c) show a much faster decay after the ultrafast build-up caused by THz excitation. The slow response observed in the other samples immediately after pump excitation is almost negligible. The relaxation time in BSTS is short—a few hundred fs—which is very different compared to Bi2Se3, Bi2Te3 and also with respect to the dynamics in BSTS observed in studies employing above-bandgap excitation10,26. In all these cases, a pronounced initial relaxation component with a timescale of a few picoseconds is present. We ascribe the observed fast dynamics to the isolated response of Dirac fermions in the surface states of our TI system. We suggest that these fast surface-state relaxation dynamics we unveiled are the result of the efficient coupling of carriers to phonons5,27 and ultrafast surface-to-bulk scattering7, and could hence be a general property of this material class. These fast dynamics are likely related to the carrier temperature that can exceed 1000 K for the fluence we used (~100 μJ cm−2), according to Ref. 7. We note that the observed decay is significantly faster than the decay of surface-state carriers reported at 5 K, which was ~1.5 ps11. The reason for faster decay at room temperature than at 5 K is most likely related to phonons—either because of cooling directly to phonons, or because of phonon-assisted surface-bulk scattering that becomes faster at room temperature.

Optical-pump optical-probe measurements

To verify independently that the ultrafast relaxation in BSTS originates from intraband dynamics within the topological SS, we have performed degenerate optical pump–probe measurements using pump pulses with a photon energy of 1.5 eV (Fig. 1c). The experimental setup is described in the Supplementary Information. This scheme leads to the simultaneous excitation of bulk and surface states in all samples. In the case of Bi2Te3 and Bi2Se3 (Fig. 1a, b), we see almost identical relaxation dynamics for optical and THz excitation, confirming that the dominant contribution to the dynamics originates from bulk states in bulk-metallic TIs. For BSTS we, in contrast to the below bandgap THz excitation, observe also a much slower dynamic on picosecond timescales in agreement with earlier observations for above bandgap excitation28. The significantly faster decay after photoexcitation with below-bandgap THz photons, in comparison with excitation with above-bandgap photons, clearly suggests that the hot surface-state Dirac fermions have a much shorter lifetime. Apart from an initial decay, the OPOP measurements in Bi2Te3, Bi2Se3 and BSTS yield qualitatively very different results. This can be attributed to different ways in which the temperature and population relaxation can compete with each other, depending on the material’s band structure and wavelength of the probe pulse10. For example, in the case of BSTS, depending on the probe wavelength, the photo-induced reflectivity variations change the sign for probe energies below 0.95 eV, and not for higher energies, such as the 1.5 eV that we used10. Note that for all investigated TIs both TPOP and OPOP measurements show negative photo-induced reflectivity at 1.5 eV probe energy.

Terahertz harmonic generation

In order to obtain additional insights into the relaxation dynamics of Dirac fermions in TI surface states, we study THz high-harmonic generation (HHG). Strong THz high-harmonic generation (HHG) has been observed in a number of Dirac materials: in graphene29,30,31, in the 3D Dirac semimetal Cd3As232,33, and in the prototypical TI Bi2Se334. Comparing the nonlinear THz response of samples with nontrivial (Bi2Se3) and trivial ((Bi0.9In0.1)2Se3) topology at the surface, it was shown that THz HHG originates from the surface states of the TI34.

The underlying mechanism for THz HHG by Dirac fermions can be described by a thermodynamic nonlinearity mechanism that relies crucially on the ultrafast modulation of the THz absorption. This mechanism is enabled by efficient heating and subsequent cooling of the electronic system on fs to ps timescales upon interaction with strong THz fields29,30,35. Therefore, the THz HHG behavior is intricately related to the characteristic timescales of carrier dynamics of Dirac fermions. In particular, we point out that when cooling takes several picoseconds, this leads to a reduction of the nonlinearity coefficient and strong saturation effects for increasing incident field strengths, as observed recently for graphene under strong THz excitation31 (see also Supplementary Fig. 12).

In Fig. 2, we show measurements of fundamental and third-harmonic (THG) fields measured by electro-optic sampling for Bi2Te3, together with the respective amplitude spectra. The details of the THz HHG experimental setup are shown in the Supplementary information. We observe a clear third-harmonic signal from Bi2Te3 at 1.5 THz (Fig. 2b) upon driving the sample with the 0.5 THz fundamental with a peak field strength of 100 kV cm−1. In Fig. 2c the resulting amplitude spectra of the fundamental beam and the THG in both TI and substrate are shown. In these experiments, the THG polarization is collinear with the polarization of the linear fundamental beam polarization. The third-harmonic field strength scales down linearly with fundamental beam ellipticity reaching zero in the case of circularly polarized pump (see Supplementary Figs. 5, 6). As we investigated the THz THG in transmission geometry, the observed harmonics can originate from both front (air–TI interface) and back surfaces (TI–substrate interface) of the sample.

Fig. 2: THz third-harmonic generation in Bi2Te3.
figure 2

Time-domain waveforms of the fundamental radiation electric field (a) and the third-harmonic field generated in Bi2Te3, and the Al2O3 substrate for reference (b). Amplitude spectra of the fundamental radiation and THG in Bi2Te3 and substrate (c).

Full size image

In Fig. 3a, we show the peak field strengths of the third-harmonic signal generated in Bi2Se3, Bi2Te3 and BSTS as a function of the field strength of the fundamental radiation. For comparison we show THz third-harmonic signal measurements from p-doped graphene with around 1013 cm−2 charge carrier concentration. For fundamental field strengths between 60 and 100 kV cm−1 the third-harmonic signal from graphene exhibits clear saturation behavior and begins to scale linearly with the fundamental field, rather than the expected cubic scaling for third-harmonic generation in the perturbative regime. Even in the lowest range of driving fields, between 15 and 25 kV cm−1, third-harmonic generation in graphene shows a power law scaling with an exponent of 2.3, indicating that the perturbative regime is only maintained at driving fields below 10 kV cm−1. These observations for graphene are in agreement with earlier publications29,31,33. The strong saturation at the highest field strengths has been ascribed to a combination of effects: the first reason for this is the occurrence of saturable absorption: if the carrier temperature increases, the graphene conductivity decreases and the absorption of THz light is reduced. Eventually, the carrier temperature can become too high such that no incident THz light is absorbed. In that case, further increasing of the field strength will not lead to any additional generation of a higher harmonic signal. The second reason for saturation is that the increase of temperature scales sub-linearly with the incident THz power36,37, because the heat capacity scales with temperature. This sub-linearity also gives rise to saturation. The third reason is that at elevated temperatures cooling slows down, because of a hot phonon bottleneck effect. With slower cooling, the heating-cooling dynamics oscillate with a smaller amplitude, thus giving rise to saturation in the HHG.

Fig. 3: Fluence dependence of the THz HHG.
figure 3

a The peak values of third-harmonic fields generated in topological insulators and graphene obtained for different driving fields of fundamental radiation at 500 GHz central frequency (dots). Solid lines are linear fits providing the scaling law: 1.06 ± 0.08 for graphene at highest field and 2.3 at lowest fields, 3.08 ± 0.03 for Bi2Te3, 2.99 ± 0.06 for BSTS, and 2.98 ± 0.07 for Bi2Se3. b Fields of fifth harmonic generated in BSTS, the slope of the linear fit is 5.2 ± 0.3; The inset shows amplitude spectra of fundamental and fifth harmonic.

Full size image

Compared to graphene, the TI samples show fundamentally different behavior: all topological insulators demonstrate a clear cubic dependence of the third-harmonic field across the whole accessible field regime meaning that THz-induced nonlinear processes are far away from the saturation regime. In order to verify this perturbative behavior further, we show in Fig. 3b the field strength of the fifth harmonic as a function of the fundamental field strength for the BSTS sample. Using a power-law fit, the fifth-harmonic (FHG) fields scale with an exponent of 5.2 ± 0.3 as a function of fundamental field, which is clear evidence that the FHG process is also in the perturbative regime. This purely perturbative behavior all the way up to incident field strengths of 140 kV cm−1 can be ascribed (partially) to the sub-picosecond cooling dynamics we observed in our TPOP experiments: the ultrafast surface-state relaxation of the TI samples prevents heat accumulation effects within the comparably long duration (several picoseconds) of the THz excitation pulse. An additional reason for the extended perturbative regime for the tested TIs could be the lower nonlinearity, which is likely due to the lower Fermi velocity compared with graphene38. This allows for driving the TI samples with higher fields before saturation is observed. The influence of both relaxation times and the Fermi velocity on the THz THG saturation is shown in the Supplementary information (Supplementary Fig. 15). According to the thermodynamic model, both relaxation and Fermi velocity determine the THz HHG saturation. Long relaxation and high Fermi velocity are required to have high THz nonlinearity at low driving fields, while at driving fields above 100 kV cm−1 ultrafast cooling is the most critical parameter to reach high THz HHG yield.

The observation of perturbative HHG in TIs is highly interesting, because it could mean that a higher conversion efficiency can be obtained for increased incident field strengths. We observe that the highest overall conversion efficiency for graphene is about 0.5% in the field, which is consistent with other studies29,31. However, saturation effects at high fields may impose an efficiency limit. For the TI samples, we observe maximum third-order conversion efficiencies of 0.13, 0.08, and 0.03% in the field for Bi2Te3, BSTS, and Bi2Se3, respectively. A maximum conversion efficiency for the fifth harmonic generation of around 0.014% at 140 kV cm−1 fundamental field strength is observed. These efficiencies are lower than the one obtained for graphene. However, the purely perturbative scaling indicates that higher conversion efficiencies than for graphene are achievable by further increasing the incident field, which is of great technological interest.

Perspective and summary

In conclusion, we experimentally isolated the dynamical response of Dirac fermions in surface states of TIs that are excited by THz light, and observe ultrafast relaxation dynamics on a timescale of a few hundred femtoseconds. This decay is significantly faster than the dynamics of photoexcited carriers in bulk states, and likely originates from efficient phonon-assisted scattering either via surface intraband or surface-bulk interband transitions. In agreement with such fast relaxation, we observe no saturation effects in THz harmonic generation up to the highest field of 140 kV cm−1, where the graphene harmonic response is already strongly saturated. Our findings are of high technological relevance since they indicate that the nonlinear conversion efficiencies in the investigated TIs can, in contrast to graphene, potentially be scaled to higher values e.g. by employing metamaterials to locally enhance the electric fields31. Furthermore, the observed fast relaxation of surface-state Dirac fermions could open up new opportunities for ultrafast spintronic devices based on topological insulators.