Probing and controlling terahertz-driven structural dynamics with surface sensitivity

Intense, single-cycle terahertz (THz) pulses offer a promising approach for understanding and controlling the properties of a material on an ultrafast time scale. In particular, resonantly exciting phonons leads to a better understanding of how they couple to other degrees of freedom in the material (e.g., ferroelectricity, conductivity and magnetism) while enabling coherent control of lattice vibrations and the symmetry changes associated with them. However, an ultrafast method for observing the resulting structural changes at the atomic scale is essential for studying phonon dynamics. A simple approach for doing this is optical second harmonic generation (SHG), a technique with remarkable sensitivity to crystalline symmetry in the bulk of a material as well as at surfaces and interfaces. This makes SHG an ideal method for probing phonon dynamics in topological insulators (TI), materials with unique surface transport properties. Here, we resonantly excite a polar phonon mode in the canonical TI Bi$_2$Se$_3$ with intense THz pulses and probe the subsequent response with SHG. This enables us to separate the photoinduced lattice dynamics at the surface from transient inversion symmetry breaking in the bulk. Furthermore, we coherently control the phonon oscillations by varying the time delay between a pair of driving THz pulses. Our work thus demonstrates a versatile, table-top tool for probing and controlling ultrafast phonon dynamics in materials, particularly at surfaces and interfaces, such as that between a TI and a magnetic material, where exotic new states of matter are predicted to exist.

evant for studying topological insulators (TIs), materials that are bulk insulators, but have a conducting surface state with a Dirac cone band structure [11,12]. THz time-domain spectroscopy has already shown the ability to separately measure the linear response of the surface and bulk states in TIs (e.g., in Bi 2 Se 3 ) [22-24], revealing both a Drude response and an infrared (IR)-active phonon. These low energy excitations play an important role in the high mobility, spin-polarized surface currents characteristic of TIs, with numerous applications in optoelectronics and spintronics [11,13], but their nonlinear response is still relatively unexplored.
A powerful approach for addressing this is to probe the THz-driven response with optical second harmonic generation (SHG). Much like x-ray diffraction (XRD), the rotational anisotropy of the SHG signal, or its intensity as a function of the crystal orientation, yields information about the lattice, electronic and magnetic symmetries of a crystal [6,7,25]. However, SHG also provides surface and interface selectivity in crystals with bulk inversion symmetry, where the bulk contribution to the SHG signal, but not the surface or interface contributions, vanishes [8]. Previous work demonstrated that static optical SHG is indeed a probe of the surface in Bi 2 Se 3 [9]. Furthermore, temporally resolving the SHG signal after optical excitation made it possible to separately track carrier [26] and phonon dynamics [10] at the surface and in the bulk, revealing distinct responses.
Here, we demonstrate a unique approach for understanding and controlling THz light-matter interactions by combining resonant excitation of a polar phonon mode with optical SHG. In the TI Bi 2 Se 3 , the resulting dynamic symmetry changes lead to a large (up to ∼50%) modulation of the SHG signal that consists of two components: one that originates from the surface, oscillating at the phonon frequency (ν p = 1.95 THz), and another at 2ν p = 3.9 THz that comes from tran-sient phonon-induced symmetry breaking in the bulk. Finally, we coherently control the phonon oscillations by tuning the time delay between a pair of driving THz pulses, enabling us to either increase or decrease the oscillation amplitude. More generally, our combination of THz excitation with SHG probing represents a new approach for directly driving structural dynamics through low energy excitations and probing the resulting ultrafast response using nonlinear optics.
The Bi 2 Se 3 sample used in these experiments was a 25 nm thick film, deposited by molecular beam epitaxy on a (0001) sapphire substrate [27]. The Fermi level in this sample is slightly above the bulk conduction band minimum. Fig. 1a shows a schematic of our THz-pump, SHGprobe experimental setup (more details are provided in Methods). The THz beam (peak E-field ≈150 kV/cm) is p-polarized and focused tightly at the Bi 2 Se 3 sample with an incident angle of θ = 30 • . After a time delay τ (relative to the peak of the THz pulse), a co-propagating and colinearly polarized 1.55 eV (800 nm) probe pulse arrives, and the p-polarized component of the reflected 3.1 eV (400 nm) SHG signal produced at the sample is collected with a photomultiplier tube (PMT). We can then measure the photoinduced symmetry changes after THz excitation by measuring the SHG signal while rotating the sample about the azimuthal angle, φ. Finally, we note that all of the measurements shown here were performed at room temperature, but data taken at lower temperatures was consistent with these results (see the supplemental information).
We began by blocking the THz pump beam to characterize the static SHG pattern, I 2ω (φ), of our Bi 2 Se 3 sample (Fig. 3a, with more detail in the supplemental information), confirming that the pattern corresponds to C 3v symmetry, in good agreement with previous measurements on single crystals [9]. We then characterized the linear THz response of our Bi 2 Se 3 film (Fig. 1b). This also agrees well with previous work [22], showing the Drude response and the IR active in-plane E 1 u phonon at ν p =1.95 THz. More detail is given in Fig. 1c, which shows a quintuple layer of Bi 2 Se 3 and illustrates the lattice motions corresponding to the E 1 u phonon.   comparison. Finally, to better interpret these experiments, we also measured the dependence of ∆I 2ω on the magnitude of the incident THz E-field E T Hz (t), which is plotted in Fig. 3b. This shows that the ν p component depends linearly on E T Hz (t), while the 2ν p component depends The most compelling observation from these figures is that resonantly driving the E 1 u phonon with a THz pulse causes the SHG intensity to oscillate in time at frequencies of ν p and 2ν p , with their relative amplitudes depending on φ (Fig. 2d). Physically, this occurs because THz excitation of this phonon transiently breaks lattice symmetry, directly modifying the time-dependent SHG signal [6,7]. However, two questions remain: do the THz-induced changes happen in the bulk or at the surface, and what is the origin of the 2ν p oscillations? It is known from previous work that the static SHG signal in Bi 2 Se 3 originates from the surface [9,28]. Resonantly driving the E 1 u phonon then changes the surface crystal symmetry, which should transiently modify the surface SHG. However, the intense THz field also penetrates into the bulk, breaking inversion symmetry and potentially "turning on" SHG there when the lattice ions are displaced from equilibrium by the phonon excitation.
For more insight, we consider the nonlinear polarization at 2ω, P 2ω , where ω corresponds to 1.55 eV and 2ω to 3.1 eV, in the presence of the THz E-field: Just as in static E-field induced SHG, the THz-field-induced changes can be described as a fourwave-mixing process that depends on the third order nonlinear susceptibility χ (3) [29, 30]. Then, calculating the SHG intensity from I 2ω ∝ (P 2ω ) 2 , [6] and computing the difference (∆I 2ω ) with and without the THz E-field, we find, in agreement with ref.
[30], that Eq. 2 consists of two terms, one linear and one quadratic in the amplitude of E T Hz , much like the data shown in Fig. 3b. As with the linear susceptibility, χ (3) can be decomposed into a sum of resonant and non-resonant contributions [6,31], which at THz frequencies in Bi 2 Se 3 are the phonon and Drude responses given by χ Drude , respectively (plus a constant originating from higher frequency optical transitions). Since E T Hz is resonant with a phonon mode (which is unique to our experiment), the phonon contribution to χ (3) will dominate the nonlinear response.
χ (3) in Eq. 2 can then be written as a product of four linear susceptibilities at each of the relevant frequencies as follows [6]: A is the anharmonicity and χ (1) phonon is the phonon contribution to the linear THz susceptibility, excluding the Drude response (i.e., the low frequency tail in Fig. 1b). This enables us to conclude that the spectral shape of χ (3) E T Hz (from Eq. 2) will follow χ phonon in Bi 2 Se 3 can be described by a complex Lorentzian [22] whose Fourier transform is proportional to B(t)cos(2πν p t), where B(t) is a decaying exponential that depends on the linewidth of the Lorentzian. Substituting this into Eq. 2 reveals that the first term oscillates at ν p , while the second term depends on cos 2 (2πν p t) and thus has both a DC component as well as an oscillatory component at 2ν p . This is consistent with the spectra shown in Fig. 2b. These two terms therefore correspond to the ν p and 2ν p components in our pump-probe measurements, respectively (Figs. 2 and 3). For more information on this derivation, see the supplemental information.
A more detailed understanding of the physical origin of these two terms comes from considering that the linear contribution contains χ (2) , which is zero in the bulk. Therefore we conclude that this component comes only from the surface. However, the second term depends only on which is non-zero in the bulk and thus must originate from the entire thickness of the film. This is the transient bulk contribution to the SHG discussed above, occuring when the photoexcited phonon breaks inversion symmetry. Eq. 3 also gives two reasons for the large magnitude of our THz-pump, SHG-probe signal: one, because the linear susceptibilities at 1.55 and 3.1 eV are very large (the dielectric constant ǫ at 1.55 eV and 3 eV is ∼ 29 from the supplemental material in [9]), and second, because the resonance at 1.95 THz substantially increases χ (1) (ω T Hz ).
The symmetry of the THz-induced SHG changes shown in Fig. 3a is consistent with that expected from the phonon. This can be seen from the schematic depicting the E 1 u phonon mode in Fig. 1c. Fixing the THz E-field direction along the a-axis while rotating the crystal (as in the experiment), each C-site Se ion reaches either a maximal or minimal displacement with respect to the B-site ions every 60 • , repeating itself every 120 • . This yields an SHG pattern with three-fold symmetry, as shown by the magenta curve in Fig. 3a. In contrast, the 2ν p component is always positive, since there is no static SHG signal from the bulk for the induced dipole moment to subtract from. This means that a transient phonon displacement can only create bulk SHG, resulting in the six peaks shown by the blue curve in Fig. 3a.
We can coherently control the E 1 u lattice vibration by driving it with a pair of THz pulses, as shown in Fig. 4, similar to what was done previously with magnons [17]. Setting the delay between the THz pump pulses to be τ T Hz =250 fs, or half of the phonon oscillation period, makes the oscillations driven by each THz pulse out of phase with one another, reducing the amplitude of the modulations in the SHG signal (Fig. 4a). Choosing instead τ T Hz =500 fs causes both THz pulses to drive the phonon oscillations in phase with one another, leading to an increase in the oscillation amplitude (Fig. 4b). We note that the oscillations do not completely cancel out ( Fig.   4a) or double in amplitude (Fig. 4b) since the two THz pump pulses do not have exactly the same time-domain waveforms and the phonon coherence time is finite. However, arbitrary control of the lattice vibrations could be achieved using a THz pulse shaper to produce an arbitrary amplitude, phase and polarization [5,32].
Finally, we note that although these measurements were performed on a topological insulator, they reveal THz-driven phonon dynamics at the surface and bulk, and contain no apparent signature Note that this means that the exact time delay between the E T Hz (t) waveform measured on the GaAs wafer and the ∆I 2ω signals measured on Bi 2 Se 3 is not known. The susceptibility shown in Fig. 1(b) was measured using THz time-domain spectroscopy. The second THz pulse used in Fig. 4    the input 1.55 eV polarization ρ and the crystal angle φ. The incident angle for the 1.55 eV beam was 45 • . The measured images are shown in Fig. S1a-b for both S and P output polarizations at 3.1 eV. Cuts from these images (Fig. S1c-f) for input S and P polarizations at 1.55 eV are in excellent agreement with that measured on a crystal after 30 min of exposure to air [1] (see the discussion in Sec. 4 about the surface accumulation layer). The P in , P out I 2ω curve shown as the orange curve in Fig. 3a of the main text, measured with the 1 kHz system used for the terahertz (THz) experiments, also shows the same SHG symmetry. The only difference was in that case, a slightly lower contrast between the minima and maxima was observed. This is likely due to the larger ∼ 200 µm 1.55 eV beam spot size used in the THz experiments, which averages over more grains (see [2] for a discussion of the grain size in Bi 2 Se 3 films). These can be directly compared to the curves of Fig. 1 in ref. [1]. These measurements were made at 300 K.

THz-pump, SHG-probe data at 10 K
The THz-pump, SHG-probe measurements (∆I 2ω ) shown in the main text were also performed at different temperatures between 300 and 10 K. The results at 10 K are shown in Fig. S2. Since we could not rotate the crystal within the cryostat at lower temperatures, we instead measured the input polarization (ρ) and pump-probe delay (τ ) dependences of ∆I 2ω . It is more difficult to extract the pump-induced symmetry change from the ρ dependence, so in the main text we present only the room temperature data, where we were able to rotate the crystal. Nevertheless, the low temperature data shows the same features seen at room temperature. It also shows oscillations in ∆I 2ω following the phonon frequency ν p as well as at 2ν p , which are lower in frequency (ν p =1.95 THz at 300 K and 1.87 THz at 10 K) and spectrally narrower at 10 K compared to 300 K, consistent with static transmission measurements (see, e.g., ref. [3] for the temperature dependence of the phonon absorption). The same linear and quadratic dependence of the ν p and 2ν p terms, representing surface and bulk contributions, respectively, was also seen down to 10 K (see Fig. S2d).

Simulating the nonlinear THz response from the linear susceptibility
As discussed in the main text and following these references [4,5], the third order nonlinear susceptibility χ (3) is a sum of resonant and non-resonant contributions. In this context, at THz frequencies we refer to the resonant contribution as χ phonon and χ (1) Drude , are comparable in size [6], the resonant contribution to χ (3) is typically an order of magnitude or more larger near a resonance [4].
Therefore, the THz-driven change in the second harmonic signal ∆I 2ω , described in Eqs. 2-3 in the main text, is dominated by the χ phonon contribution to χ (3) . Then, following Eq. 3 from the main text, χ THz pulse spectral width, and also that the linear susceptibilities at ω and 2ω are slowly varying, the spectral shape of χ phonon E T Hz will follow that of χ (1) phonon .
To show that this is indeed the case, we extracted the Drude and phonon contributions from the total χ (1) , which was determined from a linear THz transmission measurement at 300 K, as shown in Fig. 1b where d is the film thickness and ǫ 0 is the vacuum permittivity. Next, to simulate ∆I 2ω , we take the product χ (e) (f) Figure S4 | Simulating the nonlinear response ∆I 2ω . a, The incident THz E-field E T Hz (t) and b, its spectrum |E T Hz (ν)|. c, The temporal shape, and d, spectral amplitude of χ (1) phonon E T Hz (t). e, The temporal shape and f, spectral amplitude of (χ (1) phonon E T Hz (t)) 2 .
χ (1) phonon was extracted from the measured linear susceptibility, as shown in Fig. S3. The grey dashed lines in d and f are at 1.95 and 3.9 THz, respectively.
4 Accounting for the χ (3) contribution from the surface accu-

mulation layer
As mentioned in the footnote in the main text (Ref. 28) and shown in ref. [1], the static surface SHG in Bi 2 Se 3 , which originates from χ (2) , also contains a χ (3) contribution from a static electric field E s due to Se vacancies in the first ∼ 2 nm of the film thickness, referred to as the surface accumulation layer. This static electric field contributes to the SHG, since χ (2) = χ (2) surf ace + χ (3) E s . While χ (2) surf ace from surface symmetry breaking comes from only about a 1 nm thickness, the thickness of a Bi 2 Se 3 quintuple layer, the accumulation layer contribution comes from about the first 2 nm [1,2], depending on the penetration depth of E s . As observed in ref. [1], since surface oxidation increased the amplitude of the SHG by a factor of ∼ 2 and did not change its azimuthal angle dependence, the surface accumulation layer effectively just makes the region where the symmetry is broken at the Bi 2 Se 3 /air interface about twice as thick. Following the discussion in Sec. 3, we also note that χ (3) in this case is off-resonant, since the static field E s is not resonant to any phonons [8], and not modulated by the THz-excited E 1 u phonon. Considering this, if we explicitly include these two contributions to the SHG in Eqs. 2-3 of the main text, shown in Eqs S1-S2, it does not change the conclusions.
We note that phonons excitations normal to the surface could result in a modulation of χ (3) E s , effectively adding a resonant contribution to χ 3 . Such an effect would have occurred at a distinct frequency and was not observed in our measurements, since out-of-plane phonons were not resonantly excited by our THz pulse.