Ultrafast photocurrents at the surface of the three-dimensional topological insulator Bi2Se3

Three-dimensional topological insulators are fascinating materials with insulating bulk yet metallic surfaces that host highly mobile charge carriers with locked spin and momentum. Remarkably, surface currents with tunable direction and magnitude can be launched with tailored light beams. To better understand the underlying mechanisms, the current dynamics need to be resolved on the timescale of elementary scattering events (∼10 fs). Here, we excite and measure photocurrents in the model topological insulator Bi2Se3 with a time resolution of 20 fs by sampling the concomitantly emitted broadband terahertz (THz) electromagnetic field from 0.3 to 40 THz. Strikingly, the surface current response is dominated by an ultrafast charge transfer along the Se–Bi bonds. In contrast, photon-helicity-dependent photocurrents are found to be orders of magnitude smaller than expected from generation scenarios based on asymmetric depopulation of the Dirac cone. Our findings are of direct relevance for broadband optoelectronic devices based on topological-insulator surface currents.

M any efforts in current solid-state research aim at pushing the speed of electronic devices from the gigahertz to the terahertz (1 THz ¼ 10 12 Hz) range 1 and at extending their functionalities by the spin of the electron 2 . In these respects, three-dimensional topological insulators (TIs) are a highly promising material class. Although having an insulating bulk, their surface is metallic due to a band inversion that is topologically protected against external perturbations. Bi 2 Se 3 is a model TI 3 as its surface features a single pair of linear Dirac-type electronic energy bands 4 with spin-velocity locking and forbidden 180°backscattering 5 . These properties are ideal prerequisites to induce large spin polarizations by means of surface currents.
Part of this considerable potential was demonstrated by recent works, which reported the exciting possibility of launching TI surface currents by simply illuminating the sample with light [6][7][8][9][10][11][12][13] . The direction of the photocurrent could be controlled through the polarization state of the incident light beam. The assignment to a surface process was supported by picosecond time-of-flight measurements 13 showing that the photoinduced carriers were propagating at a speed comparable to the band velocity of the Dirac states. There is, however, still an intense debate about mechanisms leading to TI surface photocurrents. Scenarios based on asymmetric depopulation of the Dirac cone 6 transitions into other higher-lying cones 13 , and asymmetric scattering of electrons 8 have been proposed. To directly resolve the generation of TI surface photocurrents, we need to boost the time resolution of the experiment from so far as B250 fs and longer [9][10][11][12][13] to the scale of elementary scattering events, which can be shorter than 10 fs.
Here, we use ultrabroadband terahertz (THz) emission spectroscopy [14][15][16] from 0.3 to 40 THz to probe the ultrafast evolution of photocurrents in the Ca-doped model TI Bi 2 Se 3 with unprecedented time resolution. On the basis of an analysis of their temporal structure and symmetry, we identify distinct current sources. First, a slow drift current of photoinduced bulk charge carriers along the TI surface field is found. Second, we observe that currents depending on the pump helicity are orders of magnitude smaller than expected from the photocurrent generation scenario based on asymmetric depopulation of the Dirac cone 6 . This remarkable result suggests a strong mutual cancellation of the contributions of the various optical transitions, much reduced matrix elements for surface-to-bulk transitions and/or relatively small pump-induced changes in the electron band velocity. Finally, for the first time, we observe a new type of photocurrent, a surface shift current, which originates from an instantaneous displacement of electron density along the Se-Bi bond. This charge transfer is localized in a surface region of B2 nm thickness, the natural confinement scale of topological surface states. Its relaxation time of 22 fs provides the timescale on which the optically excited surface carriers relax to an isotropic distribution. In terms of applications, the instantaneous electric field generated by the shift current could be used to drive highly spin-polarized THz electric currents at the TI surface along an easily tunable direction.

Results
Ultrafast photocurrent amperemeter. A schematic of our experiment is depicted in Fig. 1a. A femtosecond laser pulse is incident on the specimen and launches a transient charge current density j(z,t). This photocurrent, in turn, emits an electromagnetic pulse with transient electric field E(t), in particular, covering frequencies up to the THz range, as expected from the inverse duration of the femtosecond stimulus. The measurement of E(t) over a large bandwidth (0.3 to 40 THz) permits extraction of the sheet current density with ultrafast time resolution. As detailed in the 'Methods' section, this approach allows us to separately determine the current component J x directed along the x axis and the component J yz , which is a linear combination of the roughly equally weighted Cartesian components J y and J z (see Fig. 1a and the 'Methods' section). By virtue of a generalized Ohm's law, the currents J x and J yz are, respectively, connected to the s-polarized electric-field component E x and the perpendicular, p-polarized component E yz directly behind the sample (Fig. 1a). The THz near-fields E x and E yz are obtained by measuring the THz far-field using electro-optic sampling, resulting in the electro-optic signals S x and S yz , respectively (see the 'Methods' section). THz waveforms are acquired at the shot-noise limit of the setup and for various settings of the pump polarization and sample azimuth f (Fig. 1b).
We use this approach to study a freshly cleaved, n-type, Ca-doped Bi 2 Se 3 single crystal in ambient air (see the 'Methods' section). While dipolar photocurrents in the inversion-symmetric crystal bulk (space group D 5 3d ) cancel, optical excitation can in principle launch a current at the surface (space group C 3v ) 17 . The surface region can be thought to consist of the air-crystal interface with locally relaxed lattice structure, which overlaps with the Dirac surface states (thickness of B2 nm) 18,19 , followed by a space-charge region with bent bulk bands (thickness of tens of nanometres) 20,21 . As the sample thickness (4300 mm) is much larger than the pump penetration depth (24 nm; ref. 20), the front but not the back surface of the sample is probed. Further details on sample properties and characterization can be found in the 'Methods' section.
In what follows, we will show that our broadband current measurements allow us to discriminate different types of photocurrents and their generation in the various surface regions. This goal is achieved by first identifying two dominating components in the THz emission signal using symmetry analysis. Finally, on the basis of their temporal structure and symmetry, the two underlying photocurrent components are assigned to microscopic generation scenarios.
Raw data. Typical THz electro-optic signal waveforms S(t) obtained from our Bi 2 Se 3 sample are shown in Fig. 1b. The THz waveforms depend sensitively on the setting of the THz polarization (x versus yz), the pump polarization and the sample azimuthal angle f. The signal amplitude grows linearly with increasing pump power, without any indication of saturation (inset of Fig. 1c). This behaviour implies that the number of excited carriers is proportional to the incident photon number.
As detailed in the following, we make the striking observation that the x-and yz-polarized components of the emitted THz field (and, thus, J x and J yz ) behave very differently in terms of their magnitude (Fig. 1b), temporal shape (Fig. 1b), their behaviour after sample cleavage (Fig. 2) and their dependence on the sample azimuth f (Fig. 3b). First, as seen in Fig. 1b, S yz exhibits much larger amplitude than S x but evolves significantly more slowly. Accordingly, the amplitude spectrum of S x exhibits a larger bandwidth than S yz (Fig. 1c). This trend becomes even clearer when we apply an inversion procedure to these data to extract the transient THz fields E x and E yz directly behind the sample (see the 'Methods' section). The resulting spectral amplitudes are displayed in Fig. 1d as a function of angular frequency o and show that |E x (o)| is much broader than |E yz (o)|, indicating much faster temporal dynamics. The spectrum exhibits features such as the dips of |E x (o)| at o/2p ¼ 2 THz and 4 THz whose origin becomes clear further below.
Second, to investigate the impact of surface modification on S x and S yz , we freshly cleave the sample and subsequently acquire THz signals continuously over 42 h with the sample exposed to air. While the shape of the THz waveforms does not undergo measurable modifications, their global amplitude increases by a factor of E2 in the course of time (Fig. 2). Note this rise proceeds within 30 min for S x but significantly slower (within 100 min) for S yz . We will later use this different evolution speed of S x and S yz to draw conclusions concerning the degree of surface localization of the currents J x and J yz . In contrast to S x and S yz , measurable changes of the sample reflectance at a wavelength of 790 nm are not observed, thereby ruling optical degradation of our sample out. In addition, we did not observe temporal changes in the signal symmetry, which is discussed in the next section.
Signal symmetries. In addition to their different amplitude and temporal structure, S x and S yz also depend very differently on the sample azimuth f. To quantify this behaviour, we measure waveforms S x (t, f) and S yz (t, f) for an extended set of f-values, example traces of which are shown in Fig. 3a. To reliably extract an average signal amplitude for each f, we project the time-domain signal on a suitable reference waveform (see the 'Methods' section). The resulting signal amplitude as a function of f is displayed in Fig. 3b. While both S x and S yz exhibit a 120°-periodic component of comparable magnitude, S yz has a much larger and dominant component independent of f.
The three-fold rotational symmetry of the THz signals is fully consistent with the symmetry groups of sample surface and bulk as shown by a detailed analysis of the second-order conductivity tensor (see the 'Methods' section) 3 . Importantly, it allows us to significantly reduce the large amount of experimental data contained in S(t, f): for a given THz polarization (x or yz) and pump polarization, each two-dimensional set S(t, f) can be written as a linear combination of just three basis functions (see the 'Methods' section), Therefore, three basis signals A(t), B(t) and C(t) fully characterize the entire data set S(t, f). They are, respectively, obtained by projecting S(t, f) onto the mutually orthogonal functions 1, sin(3f) and cos(3f) (see the 'Methods' section). Extracted waves We begin with considering the impact of the pump helicity on the photocurrent. The bottommost curve in Fig. 3c represents the f-independent component A x (t) of the difference of the signals taken with right-handed ( g ) and left-handed ( h ) circularly polarized pump light. The amplitude of this waveform is comparable to the noise floor. In other words, a helicitydependent yet simultaneously f-independent THz signal is small and below our detection threshold. This notion is consistent with time-domain raw data (blue versus green trace in Fig. 1b) and the absence of an offset in the f-dependence (blue curve of Fig. 3b). An analogous behaviour is observed for THz signals S yz (see Supplementary Fig. 2). We note that such small magnitude of the pump-helicity-dependent and f-independent photocurrent does not contradict the previously reported observation of timeintegrated currents 6 as will be addressed in the 'Discussion' section.
Ultrafast photocurrents. Figure 3c leads to another important conclusion of our symmetry analysis: regardless of the pump polarization, all signals S x and S yz are, respectively, dominated by just one fast and one slow waveform. We use these signals to extract the underlying source currents (see the 'Methods' section), which are displayed in Fig. 4a. After an initial onset, both J x and J yz change sign, indicating a backflow of charge. Note, however, J x proceeds much faster than J yz : the rise time from 10 to 90% of the current maximum is 16 and 40 fs for J x and J yz , respectively. Subsequently, J x decays with a time constant of 22 fs, while J yz decays within 700 fs.
To determine the origin of J x and J yz based on their ultrafast dynamics, we briefly review known photocurrent generation mechanisms 6,8,[22][23][24][25][26][27] . In general, optical excitation transfers electrons from initial states i j i into final states f j i (Fig. 4b), followed by relaxation processes such as scattering into other states, phonon emission and recombination 28 . Photocurrents can arise in both regimes, that is, during the optical transition and during the subsequent relaxation. As our pump photon energy (1.57 eV) is much larger than the Bi 2 Se 3 band gap, numerous vertical interband transitions are allowed 29 (Fig. 4f) and expected to outnumber the contribution of phonon-or impurity-assisted nonvertical transitions 30 . In the subsequent relaxation regime, currents can arise from, for instance, scattering by a noncentrosymmetric potential 8 , asymmetric recombination 31 and carrier acceleration in an intrinsic surface field (drift current) 32,33 . In all the cases, inversion symmetry needs to be broken to obtain a macroscopic dipolar net current.
Drift current. As seen in Fig. 4a, the slow current J yz (t) has a rise time (40 fs) significantly slower than the duration of the excitation pulse (E20 fs). Therefore, it cannot arise from the initial optical transition. In fact, previous works on Bi 2 Se 3 assigned the slow J yz component to a carrier drift in the surface field, consistent with the strong dependence of J yz on the doping level of Bi 2 Se 3 (refs 9-11). This notion is further supported by additional observations made in our experiment: first, the initial electron flow is directed toward the sample surface, along the direction of the space-charge field of our effectively n-doped sample. Second, the 40 fs rise time (Fig. 4a) is on the order of the bulk Drude scattering time of our sample (B18 fs, see the 'Methods' section) that limits charge acceleration in the surface field.
Following its initial rise, J yz is found to change sign. Plasma oscillation 34 of the charge carriers cannot account for this feature because the short Drude scattering time (B18 fs, see the 'Methods' section) would strongly attenuate such dynamics in less than 100 fs, in contrast to our observation. We consequently assign the sign change of J yz to the backflow of charge that accompanies the overall relaxation of the photoexcited system back to the equilibrium state. From photoemission studies 28 , the relaxation of optically excited bulk carriers is known to occur on a timescale of 1 ps, in agreement with the relaxation time of J yz (see Fig. 4a).
Shift current. Having assigned the slow, dominant part of current J yz , we now focus on the very fast, sub-100 fs dynamics of the photocurrent J x . Concerning immediate photocurrent generation by an optical transition i j i ! f j i, Sipe and colleagues 22,35 used perturbation theory and identified three distinct mechanisms: injection currents, shift currents and optical rectification 22,36 . Injection currents J inj arise because initial and final state of the perturbed electron have different band velocity. An example is the asymmetric band depopulation scenario 6 shown in Fig. 4b: a circularly polarized pump excites electrons from the Dirac cone into higher-lying states with different band slope (group velocity). Therefore, for short enough excitation, J inj should rise instantaneously to a magnitude that scales with the average Time t (ps) Time t (ps) velocity change Dv and the density N of the excited electrons. In this simplified model, the resulting current is where the initial sheet charge density s inj ¼ eNDz inj is proportional to the thickness Dz inj of the emitting sheet, and Y(t) is the unit step function. Note that relaxation processes such as electron scattering are not covered by the theory of ref. 22. We have introduced them phenomenologically by an exponential decay with time constant t inj . Backflow of electrons is diffusive 37 and ignored on the short timescales considered here. Finally, the convolution (denoted by Ã) with the pump intensity envelope I p (t) (normalized to unity) accounts for the shape of the pump pulse, resulting in a current with the typical temporal shape shown in Fig. 4c. Shift currents 38 , on the other hand, arise when the electron density distribution of the excited state f j i is spatially shifted with respect to i j i (Fig. 4d). For short excitation, this process leads to a step-like charge displacement Dx sh Y(t) whose temporal derivative is proportional to the shift current J sh . With arguments analogous to the injection case, we obtain with s sh ¼ eNDz sh . This model implies J sh initially follows the profile of I p (t) and becomes bipolar if the relaxation time t sh is comparable to or longer than the pump duration (Fig. 4e). Finally, optical rectification can be understood as a transient, nonresonantly driven charge displacement that follows the intensity envelope of the pump pulse. It arises from all transitions between the initial and final states whose energy difference is different from the incident photon energy of 1.57 eV (ref. 39). To evaluate the relative importance of optical rectification in our experiment, we compared the emitted THz amplitude of our TI sample with that of a pump-transparent ZnTe(110) crystal, which is known for relatively strong optical rectification at the pump photon energy used here 40 . We find that, normalized to the thickness of the emitting crystal, the THz signal S x from the TI sample is about two orders of magnitude larger than that from ZnTe. Therefore, optical rectification in Bi 2 Se 3 is expected to make a negligible contribution, and we are left with considering ultrafast injection and shift currents.
Note the characteristic shape of J inj and J sh is very distinct: unipolar (Fig. 4c) versus bipolar asymmetric (Fig. 4e). Having understood how the temporal shape of a current is intrinsically linked to its origin, we now look for such fingerprints in our data (Fig. 4a). Indeed, we find that the measured photocurrent J x (Fig. 4a) has bipolar asymmetric temporal shape: the signature of a shift current. In addition, fitting equation (4) to J x yields excellent agreement (Fig. 4a) for a pump duration of 23 fs, t sh ¼ 22 fs and Dx sh Dz sh E36 Å 2 . In this procedure, we use the excitation density (N ¼ 6.9 Â 10 24 m À 3 ) as inferred from the absorbed pump fluence (4 mJ cm À 2 ), the pump photon energy (1.57 eV) and the pump penetration depth (24 nm at 1/e intensity) 20 . The amplitude spectrum of the calculated J sh peaks at 6 THz (not shown), consistent with the peak position of the measured E x amplitude spectrum (Fig. 1d).
One could argue that the current J x is not a pure shift current (Fig. 4e) but still contains a component arising from a rapidly decaying injection current (Fig. 4c). We can exclude this possibility based on a symmetry analysis and photocurrent theory work 22 . For example, for linear pump polarization at 45°w ith respect to the plane of incidence, the resulting current J x is solely related to the second-order conductivity tensor element s xxy (see the 'Methods' section and Supplementary Note 1). As dictated by symmetry, this element equals s xyx and becomes, therefore, real-valued when the photocurrent frequency approaches o ¼ 0. From the microscopic analysis of Nastos and Sipe 22 , it follows that this tensor element and the related photocurrent J x are not due to an injection current, provided that photocurrent frequencies o higher than the current relaxation rate are considered (see the 'Methods' section). Surface localization. Our photocurrent measurements and analysis directly reveal an ultrafast shift current and a drift current in the time domain. It is so far, however, unclear to which extent these currents are localized at the surface. Since the photocurrent observed here is a quadratic nonlinear-optical effect (see inset of Fig. 1c), it only flows in regions where inversion symmetry is locally broken 41 . In our sample, this breaking arises from two perturbations of the inversion-symmetric bulk: (i) the surface, which extends over the first 1 to 2 quintuple layers (B1 to 2 nm) as indicated by the thickness of surface states 18 and the depth of surface lattice relaxation, 19 and (ii) the space-charge field E SC , which points along the surface normal and extends B15 nm into the depth of the sample, as implied by an estimate analogous to ref. 20.
Since the slow component of J yz is a drift current of photoexcited carriers in the space-charge field [9][10][11] , it is also localized in a depth of B15 nm. In addition, its amplitude changes during sample aging must arise from gradual modifications of E SC (Fig. 2). If J x were also dominated by E SC , its amplitude should evolve analogously to J yz following sample cleaving. In contrast, we observe that J x evolves about five times faster than J yz (see Fig. 2) such that J x cannot arise from perturbation (ii) but rather from the remaining perturbation (i). Therefore, the shift current is localized in the first one to two quintuple layers of the Bi 2 Se 3 surface. The resulting Dz sh of less than 2 nm and the above extracted estimate for Dx sh Dz sh imply that the shift distance Dx sh is on the order of 1 Å.
We note that the THz emission spectra exhibit sharp features at 2 and 4 THz (Fig. 1c,d) which coincide with the frequencies of longwavelength bulk phonon modes at 1.9, 2.1 and 4 THz (ref. 42). While the 1.9 THz mode is infrared-active and can thus absorb pump-generated THz radiation, the 4 THz mode is exclusively Raman active in the inversion-symmetric bulk. Therefore, the presence of the 4 THz feature in the emission spectrum suggests this mode is infrared-activated in the TI surface region where inversion symmetry is locally broken 43 . This effect further underlines the surface sensitivity of THz emission spectroscopy.

Discussion
Summarizing our results, we have shown that our ultrabroadband THz emission data are fully consistent with the notion that (i) the photocurrent J x arises from an instantaneous photoinduced shift of charge density by B1 Å in an B2 nm thick surface region of Bi 2 Se 3 . The displacement relaxes on a very fast timescale of 22 fs. The much slower current J yz is dominated by a drift current of optically excited carriers in the surface field. (ii) A pump-helicity-dependent and simultaneously azimuthindependent photocurrent is smaller than our detection threshold of 10 18 e m À 1 s À 1 . This assertion is also valid for other injectiontype transport scenarios such as photon-drag currents 44 . It is instructive to discuss these observations and compare them with previous works.
Finding (i) represents the first observation of a surface shift current, which was predicted by Cabellos et al. 45 very recently. We note that we also observe signatures of ultrafast surface shift currents in other TI samples, including thin films 46 of Bi 2 Se 3 and Bi 2 Te 3 (not shown). We emphasize that revealing the time-domain fingerprint of shift currents relies on the 20 fs time resolution of our experiment. Longer pump pulses can easily obscure this signature, even in materials with broken bulk inversion symmetry 36 . The surface shift current is probably the source of the linear photogalvanic surface currents that have been reported previously 6 , but not assigned.
Our results show that the displacement of bound charges occurs in a sheet with thickness Dz sh B2 nm, which is the thickness of the layer where the Dirac states are expected to dominate charge transport 18 . This notion is consistent with reports 19 showing that only the first quintuple layer exhibits inversion asymmetry on the order of 10%. The shift distance Dx sh B1 Å compares well with reported charge shifts on the order of the bond length (B3 Å) in noncentrosymmetric semiconductors 38 .
The three-fold azimuthal symmetry of J x (Fig. 3b) suggests the electron density is displaced along the 120°-ordered p-type Se-Bi bonds 47 . In fact, previous studies have shown that the electron density associated with the Dirac states is known to shift gradually from Se toward Bi atoms when energies below and above the the Dirac point are considered 3,47,48 . As Bi and Se atoms lie in different layers, the Bi-Se bond forms an angle of about 45°with respect to the sample surface normal. Therefore, the shift current also has a z-component with a strength comparable to J x , consistent with the sharp peak present in J yz at t ¼ 0 fs (see Fig. 4a and Supplementary Fig. 1d).
In principle, the ultrafast Se-Bi charge transfer can be driven by all kinds of optical transitions involving surface states. Examples of possible transitions induced by the 1.57 eV pump photons are shown in Fig. 4f: (1) surface-to-bulk, (2) bulk-tosurface and (3) surface-to-surface (intercone) transitions. Following excitation, the charge carriers undergo multiple scattering processes, finally resulting in an isotropic distribution. The time constant (22 fs) of this process found here is consistent with other measurements of anisotropy relaxation at the TI surface: ultrafast optical Kerr effect (time constant of B25 fs) 49 and equilibrium transport experiments (substantially smaller than 300 fs, see Supplementary Note 2) 50,51 . Note that due to spinmomentum locking at the TI surface, the decay of the anisotropy of the charge-carrier distribution is strongly connected with the decay of the transient surface spin polarization 52 . The timescales of these processes are much shorter than that of the energy relaxation of the surface charge carriers (B1 ps; ref. 28).
Result (ii), the absence of a pump-helicity-dependent photocurrent, is surprising and imposes significant constraints on the generation mechanism of this current. We first note this result is consistent with the photocurrent magnitudes found by previous electrode-based time-integrating 6 and picosecond-resolved 13 measurements. From these works, it follows that under excitation conditions similar to ours, the pump-helicitydependent photocurrent reaches a peak value that is slightly below the detection threshold 10 18 e m À 1 s À 1 of our setup (see Supplementary Note 3).
It is instructive to compare the upper photocurrent limit set by our experiment to a recently suggested microscopic scenario 6 in which the pump-helicity-dependent photocurrent arises from asymmetric depopulation of the Dirac cone by optical transitions into rapidly decaying bulk states (Fig. 4b). On the basis of this injection-type scenario, we use equation (3) to estimate the initial ballistic sheet-current density as Nev D Dz D , where v D ¼ 0.5 nm fs À 1 is the band velocity in the Dirac cone 28 , Dz D ¼ 2 nm the thickness of the Dirac states 18 and N is the bulk excitation density. The resulting magnitude of 10 22 e m À 1 s À 1 is four orders of magnitude larger than the maximum current measured in our experiment.
Possible reasons why the photocurrent magnitude predicted by this plausible scenario so drastically exceeds the actually measured pump-helicity-dependent photocurrent are as follows. First, matrix elements for bulk-surface optical transitions are much smaller than for bulk-bulk transitions. In other words, light absorption is much more pronounced in the bulk than at the surface, in contrast to the tacit assumption of homogeneous absorptivity made in the estimate above. Second, the band velocity of initial and final state are approximately equal (see equation (3)). Finally, there is a great deal of cancellation when summing over the contributions of all optical transitions. Indeed, as shown theoretically, a zero net current results when only optical transitions within the Dirac cone are considered 53 . Experiments using pump pulses with tunable pump photon energy 54,55 will likely provide more insights into the nature of the pump-helicity-dependent photocurrent, in particular, in terms of intracone and below-gap optical transitions 24,53,54 .
In conclusion, we have measured the dynamics of ultrafast photocurrents on the surface of the Ca-doped three-dimensional model TI Bi 2 Se 3 with a time resolution as short as 20 fs. We find that the peak amplitude of pump-helicity-dependent photocurrents is much smaller than predicted by a recent model based on asymmetric optical transitions between Dirac-cone and bulk states 6 . For the first time, we have observed a surface shift photocurrent, which arises from a charge displacement at the TI surface. The fast decay time (22 fs) of this current attests to the rapid loss of charge-carrier surface anisotropy. On a more applied note, the shift current is potentially interesting for ultrafast optical manipulation of the TI surface, thereby ultimately modifying its topological properties 56 . The local electric field generated by the shift current could be used to drive highly spin-polarized THz electric currents at the TI surface. The direction of this secondary current could easily be controlled by the pump-beam polarization, thereby opening up new possibilities to access the TI surface states. Finally, our results highlight broadband THz emission spectroscopy as a novel and highly sensitive probe of surfaces.

Methods
Sample details. Single crystals of Ca-doped Bi 2 Se 3 were grown by the Bridgman-Stockbarger method by pulling a sealed quartz ampoule in a vertical temperature gradient. A fresh surface is obtained by cleaving using adhesive tape. The dimensions of the sample used for THz emission spectroscopy are approximately 2 Â 3 Â 0.3 mm 3 .
For characterization of the sample by angle-resolved photoelectron spectroscopy (ARPES), the sample was studied directly after cleaving under ultrahigh-vacuum conditions (pressure 2 Â 10 À 10 mbar) and, a second time, after exposing the sample to an ambient-pressure atmosphere of N 2 for a few seconds. ARPES measurements on the N 2 -exposed surface confirm the presence of Dirac surface states with the Fermi energy located roughly 100 meV above the bulk conduction band minimum (see Supplementary Fig. 3). From these data, a conduction-band electron mass of 0.115 bare electron masses is inferred 57 .
Hall measurements 58 yield a bulk hole density of 1.34 Â 10 17 cm À 3 and a mobility of 275 cm 2 V À 1 s À 1 . Using these values, the effective mass inferred from ARPES data and the Drude formula, we extract a velocity relaxation time of 18 fs in the bulk material.
Ultrafast THz emission setup. Laser pulses (duration of E20 fs, centre wavelength of 790 nm, energy 1 nJ) from a Ti:sapphire oscillator (repetition rate of 80 MHz) are focused onto the sample (beam diameter of 200 mm full-width at half intensity maximum) under 45°angle of incidence, resulting in an average intensity o0.3 kW cm À 2 , well below sample damage threshold. The specularly emitted THz pulse is focused onto an electro-optic crystal in which the THz electric field is detected by broadband electro-optic sampling 40 . We use a (110)-oriented GaP crystal (thickness of 250 mm) owing to its relatively flat and broadband response function 16 . The only exceptions are the measurements of the two-dimensional data set S(t, f) ( Fig. 3 and Supplementary Fig. 1), which are sped up by using a (110)-oriented ZnTe crystal (thickness of 300 mm), which exhibits an enhanced detector response at the expense of reduced bandwidth. To calibrate the direction of the measured photocurrent, we use a photoconductive switch as a reference emitter in which the direction of the initial photocurrent burst is along the direction of the external bias field.
Optical wave plates are used to set the polarization state of the pump pulse to linear (with arbitrary rotation angle) or circular. Here, it is essential to avoid any optically birefringent elements (including mirrors) between wave plate and sample. A THz wire-grid polarizer (field extinction ratio of 10 À 2 ) allows us to measure the x-and yz-components E x and E yz of the THz electric field separately, thereby disentangling current components J x and J yz , the latter being a linear combination of J y and J z (see Fig. 1a and equation (8)). To ensure the electro-optic THz detector has an identical response to E x and E yz , a wire-grid polarizer with 45°orientation is placed in front of it. Variation of the sample azimuth f is performed by attaching the sample to a computer-controlled rotation stage.
From THz signals to THz fields. To proceed from the measured electro-optic signal S(t) to the THz electric field E(t) directly above the sample surface, we note that there is a linear relationship between the two quantities. For example, in the frequency domain, the THz field component E x and the corresponding signal S x are connected by the transfer function h(o) through the simple multiplication A completely analogous relationship is valid for E yz and S yz . We measure the transfer function of our setup by using 50 mm thick GaP(110) as a reference emitter placed before the TI sample that has been substituted by an Ag mirror. Details of the shape of the transfer function are shown in Supplementary Fig. 4 and discussed in the Supplementary Note 4.
From THz fields to THz currents. To finally obtain the source current J(t) from the THz electric field E(t) measured directly above the sample surface, we make use of the following generalized Ohm's law 34 : Here, o/2p is the THz frequency, Z 0 E377 O is the vacuum impedance, n(o) is the refractive index of Bi 2 Se 3 taken from ref. 50, a ¼ 45°is the angle of incidence, and is a weighted sum of the currents flowing along the y and z directions. For Bi 2 Se 3 and frequencies above 5 THz, the weighting factor sin a= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 À sin 2 a p of J z is on the order of 0.3. The inverse Fourier transformation of the resulting current spectra yields the currents in the time domain (Fig. 4a).
We note that the overall magnitude of the extracted currents is subject to an estimated uncertainty on the order of 3, which arises from the cumulated effects of uncertainties of beam parameters (such as beam diameter and divergence), of the precise optical material properties (in particular, close to phonon frequencies) and from the deviations of the optical and THz beam from a perfect Gaussian profile.
Waveform mean amplitude. To characterize a complete THz waveform S(t) by a single mean amplitude S, an often-used procedure is to calculate the root-mean square R dt S 2 t ð Þ Â Ã 1=2 . The drawback of this method is that it is based on nonlinear operations, which annihilate phase information and make the identification of additive signal components difficult. Therefore, calculation of a mean amplitude S of a given signal should be linear with respect to S(t). Our solution is to project S(t) on a suitable and normalized reference waveform S ref (t) by means of the scalar product where N ¼ hS ref , S ref i 1/2 normalizes S ref /N to unity. As required, this operation is linear with respect to S(t) and, therefore, does not mix up additive signal components. In the case of our two-dimensional data set S(t, f), we use the most intense signal (with respect to f) as reference. This choice is arbitrary, but we checked that other reference waveforms yield qualitatively identical results within our signal-to-noise ratio. According to equations (2) and (9), the basis functions A(t), B(t) and C(t) of the data set S(t, f) are then obtained by multiplying S(t,f) with 1/2p, sin(3f)/p and cos(3f)/p, respectively, and subsequent integration from f ¼ 0 to 2p.
Azimuthal symmetry analysis. The inset of Fig. 1c shows that the THz signal grows linearly with the pump power, that is, with the square of the pump field F(t). Therefore, the resulting current density can phenomenologically be described by the general nonlocal relationship 59 Here, F j is the jth Cartesian component of the pump field, and s ijk is a third-rank tensor field describing the quadratic material response. Upon sample rotation described by the matrix R i 0 i , the nonlinear response function transforms according to When we focus on rotations about the sample normal by an angle f, the matrix elements R i 0 i are given by linear combinations of the functions 1, cosf and sinf or, equivalently, 1, exp(if) and exp( À if). Since the R i 0 i show up in third order in the response transformation described by equation (11), the transformed s 0 i 0 j 0 k 0 f ð Þ is a third-order mixed polynomial with respect to 1, exp(if) and exp( À if), that is, a linear combination of the terms exp(imf) where the integer m runs from À 3 to þ 3. From these seven terms, however, only exp( À 3if), exp(3if) and 1 remain since our sample is invariant under rotation by 2p/3 ¼ 120°about the surface normal. In other words, s 0 i 0 j 0 k 0 f ð Þ is a linear combination of the functions 1, sin(3f) and cos(3f).
As seen from equations (1) and (6)(7)(8), the THz field emitted by the current-density distribution j is linear with respect to j with a proportionality constant that depends on the material only through the linear-optical constants at THz frequencies. Since for our sample these constants are invariant under azimuthal rotations, the f-dependence of the THz signal is inherited by that of j and, thus, s 0 i 0 j 0 k 0 . To summarize, the THz signal depends on the sample azimuth f through a linear combination of the three terms 1, sin(3f) and cos(3f) (see equation (2)). This conclusion is consistent with our experimental observations (see Fig. 3b and Supplementary Fig. 1c).