Nonlinear optical observation of coherent acoustic Dirac plasmons in thin-film topological insulators

Low-energy collective electronic excitations exhibiting sound-like linear dispersion have been intensively studied both experimentally and theoretically for a long time. However, coherent acoustic plasmon modes appearing in time-domain measurements are rarely observed due to Landau damping by the single-particle continua. Here we report on the observation of coherent acoustic Dirac plasmon (CADP) modes excited in indirectly (electrostatically) opposite-surface coupled films of the topological insulator Bi2Se3. Using transient second-harmonic generation, a technique capable of independently monitoring the in-plane and out-of-plane electron dynamics in the films, the GHz-range oscillations were observed without corresponding oscillations in the transient reflectivity. These oscillations were assigned to the transverse magnetic and transverse electric guided CADP modes induced by the evanescent guided Lamb acoustic waves and remained Landau undamped due to fermion tunnelling between the opposite-surface Dirac states.

H eavily doped and highly photo-excited three-dimensional (3D) semiconductors with high carrier densities (B10 17 -10 19 cm À 3 ) can sustain quantized longitudinal collective oscillations (plasmons) resulting from Coulomb interactions. These oscillations occur at the plasma frequency o p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ne 2 =me 0 p at vanishing wave vectors (q-0), where n, e, m and e 0 are density, charge, mass of electrons and the permittivity of free space, respectively. Consequently, o p lies in the far-infrared/infrared range since the carrier density is low compared to 3D metals (B10 23 cm À 3 ) where o p lies in the ultraviolet range. These plasmon modes and their interactions with lattice phonons in semiconductors have been extensively studied using optical spectroscopy [1][2][3][4][5][6][7][8] and were usually referred to as conventional 3D plasmons with a parabolic dispersion oðqÞ ¼ o p þ 3u 2 F q 2 10o p , where u F is the Fermi velocity. In contrast, the conventional 3D surface plasmon (polariton) is a self-sustaining oscillation of the semi-infinite electron gas whose surface plasmon frequency is , where e sur is the dielectric constant of the surrounding medium 9 . The resulting frequencies of conventional 3D surface plasmons in metals and conventional localized 3D surface plasmons in metal nanostructures may appear in the visible region, allowing for numerous photonic applications 10 .
3D metals can also support monolayer-thick surface electronic states that form a two-dimensional (2D) electron-density layer 11 . 2D plasmons with energy vanishing as q-0 have been observed in charge inversion layers of p-type Si(100) 12 , quasi-2D systems of monoatomic metal layers on semiconductors [13][14][15] , and extrinsic free-standing graphene 16 . Initially it was proposed that 2D electron layer with density n 2D could support only high-energy optical plasmons with o OP ðqÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n 2D e 2 q=2me 0 p dispersion [12][13][14][15][16][17] , because low-energy acoustic plasmons (APs) with sound-like o AP (q) ¼ u p q dispersion (where v p is the sound phase velocity) are expected to be Landau damped by the underlying 3D electrons 18 . Alternatively, the excitation of APs has been theoretically predicted for two-component 3D plasmas with a large difference in masses of the light and heavy particles 19,20 and for two-component 2D spatially separated plasmas 18 . In both cases, APs are expected to behave as high-frequency modes with plasmon energy (:o AP ) dispersing away from the single-particle continua. APs are also expected to be excited on 3D metal surfaces due to the nonlocality of the 3D dynamical response 21 . Because APs may mediate the Cooper pairing of particles, the AP concept has been used in theoretical studies of superconductivity in high-Tc cuprates 22,23 , MgB 2 structures 24 and layered structures 25,26 .
The topological insulator (TI) Bi 2 Se 3 , having an insulating gap in the bulk (B0.3 eV) and intrinsic metallic-type Dirac surface states (SS) 34,35 , resembles the artificially created quasi-2D systems of monoatomic metal layers on semiconductors [13][14][15] , thus offering similar conditions for 2D plasmon excitation. In principle, low-energy (GHz-range) coherent APs for the two component Dirac fermions (that is, coherent acoustic Dirac plasmons (CADPs)) may exist due to a very anisotropic mass tensor of Bi-related materials (such as Bi and Bi 2 Te 3 ) 36 . However, this tensor for Bi 2 Se 3 is almost isotropic 37 , causing the CADPs to be sharply Landau damped by the high-density single-particle continua resulting from natural n-doping. This statement is also applicable to thin Bi 2 Se 3 films thicker than d ¼ 15 nm, despite being thinner, two THz-range (1 THzB4.14 meV) plasmon modes ('optical' and 'acoustic') can be excited, similarly to double layer structures, such as bilayer graphene [38][39][40][41] . The latter plasmon modes can be detected using the nonlinear optical techniques of four-wave mixing, difference frequency generation and stimulated Raman scattering [41][42][43] . These techniques deal with the anisotropy between in-plane and out-of-plane tensor components of nonlinear susceptibilities and hence avoid the photon-plasmon momentum mismatch that in the linear optical techniques is usually compensated using the waveguide evanescent Otto/Kretschmann prism couplers or grating structures 44,45 . Because the phase-matched plasmon frequency in the latter case is exclusively determined by the micro-ribbon width and the grating period 46,47 , the linear optical techniques seem to be less flexible and less sensitive to the actual (intrinsic) plasmon dynamics in the films, especially when d matches the range corresponding to direct (wave-functions overlap) and indirect (electrostatic) opposite-surface coupling (do6 nm and do15 nm, respectively) 41,[48][49][50][51] . It is worth noting that the term 'acoustic' used for the double layer structures is not always related to the excitation of real acoustic waves and is used only to distinguish between the antisymmetric ('acoustic') and symmetric ('optical') modes. These considerations indicate that unique conditions are required to observe CADPs in thin Bi 2 Se 3 films.
Here we report on the observation of transverse magnetic (TM) and transverse electric (TE) guided CADP modes in thin films of the TI Bi 2 Se 3 . A combined ultrafast linear-nonlinear optical pump-probe technique was used to simultaneously measure TR and transient second-harmonic generation (TSHG) in the reflection geometry exploiting various incident/outgoing probe light polarizations 52,53 . In addition, TSHG rotational anisotropy (TSHGRA) was measured as a function of pump-to-probe delay time. We show that TSHG and TSHGRA techniques allow for independent monitoring of temporal electron dynamics confined to the topmost Se atomic layer along the in-plane and out-of-plane directions of the films, thus providing evidence for the excitation of TM-and TE-guided CADP modes, respectively. We also demonstrate that the excitation of Landau undamped CADPs occurs at a unique condition which has been theoretically predicted for the coupled bilayer structures 54 and can be associated with fermion tunnelling between the opposite-surface Dirac SS. For a driving source of TM-and TE-guided CADP modes, we consider the strong resonant coupling of Dirac fermions to evanescent (non-propagating) guided Lamb waves excited within the same relaxation process 55 . This kind of acousto-plasmonic control is known to occur in complex nanostructures 56 , while the resonance-type acousticplasmon-to-acoustic-phonon coupling has been considered theoretically only for Bi 36 .

Results
The origin of oscillations in TR and TSHG signals. The TR and TSHG traces measured simultaneously for the 10 nm thick Bi 2 Se 3 film at various pump/probe powers and in the P in À S pump À P out and S in À P pump À S out polarization geometries are shown in Figs 1 and 2, respectively, where P and S denote the polarization of the incident probe and pump laser beams ('in' and 'pump') and the outgoing reflected fundamental/SHG probe beam ('out') in the plane of incidence (xz) and in the plane of the film (xy), respectively. The positive sign of the TSHG response refers to the quadratic form of the SHG intensity that is governed by both the purely surface and depletion-electric-field-induced bulk contributions because of the centrosymmetric nature of Bi 2 Se 3 52,53,57,58 . The SHG intensity can hence be expressed as where A is a proportionality coefficient, w S ð Þ ijk and w B ð Þ ijkl are the components of the surface second-order and bulk third-order susceptibility tensors, E dep dc is the dc(direct current)-type depletion electric field of frequency much lower than that of the driving electric field of the incident light E(o), and is the intensity of the incident probe beam. Consequently, the TSHG response monitors the variations of the E dep dc strength due to the dynamical spatial redistribution of photoexcited carriers in close proximity to the surface. Specifically, for the actual free carrier density in the film, E dep dc extends inward towards the film by B6.5 nm, thus being comparable to the incident and SHG light penetration depth (B10 and B7 nm, respectively) 41,[48][49][50][58][59][60][61][62] . This field can also be modulated by the periodic electric field associated with CADPs. Owing to indirect opposite-surface coupling in the film [48][49][50] , this modulation will always contribute to the nonlinear susceptibilities independently of whether the TM-or TE-guided mode is excited. However, by choosing the appropriate light polarization geometry, the in-plane and out-of-plane tensor components contributing to the TSHG response can be separately monitored. In contrast, an isotropic TR response is negative due to an absorption bleaching process associated with Pauli blocking, thus exposing exclusively the carrier population dynamics mainly through the out-of-plane refractive index modulation as a consequence of a normal incidence pumping geometry applied and the same spot size of the focused pump and probe beams (see the methods section) 48,55,62,63 .
Despite the different nature of TR and TSHG, both signals show a rise-and multiple decay-time behaviour that can be characterized by the corresponding rise-time (t R ) and decay-time constants (t D1 , t D2 and t D3 ) and by peak (TR-peak and TSHG-peak), background (TR-BG and TSHG-BG) and unrecovered background (TR-UBG and TSHG-UBG) intensities (Figs 1 and 2). We assigned the time constants of the TR signals to electron-electron thermalization occurring in both the bulk states and Dirac SS t TR À Á , which leads to a metastable population of the conduction band edge that continuously feeds a non-equilibrium population of Dirac SS alternatively to the carrier recombination t TR D2 À Á , and quasi-equilibrium carrier population in Dirac SS t TR D3 À Á 48,55,62 . The corresponding short time-scale rise-and decay-time constants for the TSHG responses have other physical interpretations, as discussed further below. The TR-peak, TR-BG, TR-UBG, TSHG-peak, TSHG-BG and TSHG-UBG intensities all increase with increasing probe beam power, but at different rates in accordance with their linear (for TR processes) and quadratic (for TSHG processes) power dependences 53,62 . The TR traces are similar to those reported previously 48,55,62,63 , while the TSHG traces show an oscillatory behaviour that is not present in the TR traces.
Specifically, the Fourier transform of the TSHG oscillations observed in the P in À S pump À P out polarization geometry yields a frequency B42 GHz (B0.174 meV) ( Supplementary Fig. 1). Using this frequency, one can fit the oscillatory part to a damped cosine function to yield a B50 ps damping coefficient (Fig. 1). The oscillations in the TSHG traces measured in the S in À P pump À S out polarization geometry are still present, although they occur at a lower frequency of B15 GHz (B0.062 meV) ( Supplementary Fig. 1) and with a larger damping coefficient of B80 ps (Fig. 2). In general, both the frequency and damping coefficient of oscillations observed in the P in À S pump À P out polarization geometry can be associated with the coherent acoustic phonon modes. However, these modes of similar frequencies were observed in the TR traces of the films with dZ40 nm, whereas the frequencies were significantly increased with decreasing d in the range from 40 to 15 nm because the film bulk acoustic wave resonator (FBAWR) modes begin contributing to the dynamics (Fig. 3a). Moreover, the FBAWR modes disappear completely for films with do15 nm (the range where oscillations in the TSHG traces appear) due to indirect opposite-surface coupling that leads to a gradual elimination of the out-of-plane refractive index modulation as the FBAWR regime crosses over to that of the evanescent guided Lamb wave excitation (Fig. 3a,b) 55 .
This change of the interaction regime with decreasing film thickness, which is usually associated with the contribution of the imaginary roots of the Rayleigh-Lamb equation, localize the reactive power which oscillates along the in-plane and out-of-plane directions due to symmetric and antisymmetric displacements of atoms with respect to the mid-plane of the film, respectively 64 . Because of the expected strong coupling of the resulting evanescent fields to Dirac fermions 36,56 , the corresponding TM-and TE-guided CADP modes can be excited with well-defined acoustic wavevectors and velocities (B2 Â 10 3 ms À 1 ). It is worth noting that the acoustic velocities Step size 1 ps Step size 250 fs are significantly suppressed in this case compared to the Fermi velocity in the TI Bi 2 Se 3 single crystals (B5 Â 10 5 ms À 1 ) due to the strong indirect opposite-surface coupling regime in the film 39,55 . This kind of coupling between the evanescent fields and Dirac fermions is known as surface plasmon resonance, despite the evanescent fields being created by the evanescent guided Lamb wave instead of the incident light wave as in the more usual Otto/Kretschmann prism couplers or grating structures where the phase-matched plasmon frequency is determined by the incident light wavelength and the grating period, respectively 44,45 . It should also be emphasized here that the low-frequency oscillations observed in the S in À P pump À S out polarization geometry are at least B2.5 times lower in frequency than any coherent acoustic phonon modes observed for thin Bi 2 Se 3 films 55 . Thus, oscillations in TSHG traces measured in the P in À S pump À P out and S in À P pump À S out polarization geometries for the 10 nm thick film can be assigned to the modulation of E dep dc by a periodic electric field associated with TE-and TM-guided CADPs oscillating along the out-of-plane and in-plane directions of the films, respectively. These oscillations should be distinguished from those resulting from the compressive/tensile strains induced by coherent phonons which were assigned to the local breaking of the crystal lattice symmetry 65 or those being due to an electric field modulation resulting from the interference between the incident laser beam and that reflected from propagating acoustic waves 66 . The latter processes are known to appear with identical oscillation frequencies in both TR and TSHG responses 65    TSHGRA patterns and depletion electric field screening. The TSHGRA patterns measured in the P in À S pump À P out polarization geometry reveal a spectacular change of threefold-tosixfold rotational symmetry with delay time (Fig. 1). The initial threefold rotational symmetry at the TSHG-peak intensity is similar to that observed for the stationary SHG responses from the single-crystalline and thin-film Bi 2 Se 3 samples measured in the P in À P out polarization geometry 53,57,58 . Because of the multicomponent nature of the SHG response involving hyperpolarizabilities of both the planar hexagon-arranged topmost Se layer 41,49 and Bi-Se bonds arranged into the rhombohedral unit cell along the threefold trigonal axis 53 , the dominating threefold rotational symmetry points out that the depletion-electric-field-induced bulk contribution dominates over the purely surface contribution. This statement is consistent with the corresponding TSHG intensity which can be expressed similarly to the stationary SHG intensity as 53 where B ¼ AI 2 (o) is a proportionality constant for the given experimental conditions, j is the crystal surface rotation angle, and the weight coefficients c i represent the partial contributions of the isotropic, threefold and sixfold rotational symmetry components, respectively. It should be noted that a more precise comparison of the TSHGRA patterns for the stationary SHG response and those taken at the TSHG-peak intensity suggests that the threefold rotation symmetry component is even enhanced in the latter case compared to the sixfold one. This enhancement will be discussed further below in more detail. Because the isotropic component is negligible 53 , the observed temporal dynamics of the TSHGRA patterns point to a suppression of the out-of-plane contribution during a characteristic time of B15 ps when the metastable population of the bulk states becomes high enough to completely screen E dep dc , that is, the exact sixfold rotational symmetry of the TSHG response from the topmost hexagonally arranged Se-Se bonds can be observed.
The rotational symmetry of the TSHGRA patterns measured for the 10 nm thick film in the S in À P pump À S out polarization geometry remains sixfold for any delay times applied and is identical to that observed for stationary SHG measured in the S in -S out polarization geometry (Fig. 2) 53,57,58 . The timeindependent TSHGRA patterns are consistent with the monocomponent nature of the TSHG response originating exclusively from the in-plane (xy) hyperpolarizability of the continuous hexagonal network of Se-Se bonds 53,57,58 , Because of the same origin of the sixfold rotational symmetry component of TSHGRA patterns measured in both light polarization geometries applied, their relative 30°rotation (Figs 1 and 2) results from the cofunction identity cos 2 (3j) ¼ sin 2 [3(j þ 30°)] (equations (2) and (3)) 53 . We note also here that (i) the single atomic layer nature of the TSHG response eliminates from consideration all the light re-absorption effects on the damping of the guided CADP modes, (ii) similar variations in the SHG rotational symmetry, although with increasing laser power, has been observed for the Si(001)-SiO 2 system due to the carrier-induced screening of the dc interfacial electric field 67 , and (iii) the TSHG response from Bi 2 Se 3 single crystals measured in the same polarization geometry was too weak to recognize any temporal changes in TSHGRA patterns 68 . The aforementioned dynamical screening of the depletion electric field and the resulting single atomic layer nature of the TSHG response unambiguously prove that TSHG oscillations are sourced by the Dirac fermion collective excitations.
A film thickness dependent condition for CADPs excitation. One of the intriguing findings is a very narrow d range where CADPs can be efficiently excited and the temporal threefold-tosixfold change of rotational symmetry occurs. Specifically, both effects are maximized for the 10 nm thick film, while being significantly diminished for thicker and thinner films (although still observable for the 8 and 12 nm thick films) ( Supplementary  Figs 1-3). This behaviour agrees with the d dependences of the TSHG-peak, TSHG-BG and TSHG-UBG intensities (Supplementary Fig. 3) and is consistent with the resonance-type enhancement previously reported for the stationary SHG intensity, which has been associated with the nonlinear excitation of plasmons in indirectly opposite-surface coupled Dirac SS 53 . However, these modes of double layer structures ('optical' and 'acoustic'), which are also known for the electrostatically coupled graphene bilayers [38][39][40] , have much higher frequencies B2.01 and B7.56 THz (B8.0 and B31 meV) 41,49 and should not to be confused with the CADPs driven by strong coupling between Dirac fermions and the evanescent guided Lamb waves discussed here. The coexistence of two kinds of plasmons in the TI Bi 2 Se 3 films is similar to that occurring on metal surfaces 11,21 and graphene 16 . The resonant enhancement of the stationary SHG and TSHG signals are hence due to high-energy plasmon excitation. In contrast, low-energy CADPs manifest themselves within the carrier relaxation dynamics when some unique thickness dependent condition is reached at which Landau damping of CADPs in Bi 2 Se 3 films becomes ineffective.
We associate this condition with the fermion tunnelling between the opposite-surface Dirac SS, which accompanies the indirect coupling regime for films with do15 nm and is naturally limited by the direct coupling regime for films with d o6 nm. The fermion tunnelling is expected to modify the linear AP dispersion o AP (q) ¼ u p q to where the plasmon energy gap D ¼ o AP (q ¼ 0) and the coefficients C 1 and C 2 nontrivially depend on the electron density, the fermion tunnelling amplitude, and the film thickness 54 . It should be especially emphasized here that despite this form of AP dispersion that has been suggested for 'acoustic' (antisymmetric) high-energy plasmon mode in the coupled bilayer structures, it can be directly applied to CADPs induced by the evanescent guided Lamb wave because the same dynamical dielectric function characterizes both phenomena. Using the frequencies of the FBAWR modes experimentally observed in these films 55 , and taking into account that the wavevector for the fundamental FBAWR mode is given by q ¼ 2p/l ¼ p/d, where l ¼ 2d is the acoustic wavelength (Fig. 3b), the FBAWR mode dispersion can be constructed as shown in Fig. 3a. The resulting dispersion is linear o FBAWR (q) ¼ u p q as q-0, whereas it breaks down for larger q when the evanescent guided Lamb wave regime becomes dominant (Fig. 3b). Assuming that the evanescent guided Lamb waves completely control the wavevectors and frequencies of collective electronic excitations in Dirac SS and taking into account the fermion tunnelling between the opposite-surface Dirac SS, the AP dispersion may reveal a minimum 54 . Figure 3a shows a qualitative modelling of the AP dispersion with D ¼ 0.92 meV and coefficients C 1 ¼ 5.1 and C 2 ¼ 7.684, where the patterned area represents the Landau damping region. Note that the values for C 1 and C 2 are not necessarily unique and are used here for qualitative purposes only, whereas the value of D is similar to the theoretically predicted value of B1 meV 54 . The energies of the observed CADP modes are therefore expected to occupy the area slightly above the minimum. Although this crude qualitative modelling requires further significant theoretical efforts to obtain a more complete explanation for the existence of CADPs, the model does provide a realistic explanation of why CADPs remain Landau undamped in the very narrow film thickness range centred at dB10 nm ( Supplementary Fig. 1). It is important to note that because the evanescent Lamb wave field is extended over the entire film and therefore affects equally all the free carriers independently wherever they are located (Fig. 3b), all the plasmons excited by the evanescent Lamb wave field are expected to be Landau damped on the very short time-scale of a few tens of fs, except for Dirac plasmons that remain undamped due to fermion tunnelling between the opposite-surface Dirac SS (Fig. 3a). When fermion tunnelling becomes inefficient for films thinner that B8 nm and thicker than B12 nm, Dirac plasmons are also damped. Because the TSHG technique that we used is sensitive exclusively to Dirac SS, and because the fermion tunnelling condition uniquely occurs in a certain thickness range of Bi 2 Se 3 films, the observation of CADPs becomes possible. Although the resonance-type acoustic-plasmon-to-acoustic-phonon coupling has been considered theoretically only for Bi so far 36 , we assume here that similar strong coupling may occur in Bi 2 Se 3 as well. At the very least, our experimental results clearly point out that this strong coupling does exist in a similar way as that observed in complex nanostructures 56 . We note also that because of the unusual behaviour of the 10 nm thick film, we grew a second sample using the same growth conditions, which also had similar resonant characteristics.
Ultrafast carrier SS-bulk-SS vertical transport. The simultaneous measurements of the TR and TSHG signals allow for a comprehensive study of the ultrafast carrier dynamics in thin films of the TI Bi 2 Se 3 , which shed light on the aforementioned screening of E dep dc and conditions under which CADPs can be excited. Figure 3c shows the normalized TR and TSHG traces for the 10 nm thick film. As we discussed above, the TR signal increases with t TR R $ 0:3 ps due to the electronelectron thermalization and subsequently decreases with t TR D1 $ 2:15 ps due to the electron-LO-phonon scattering 48,55,62,63 . In contrast, the TSHG signal increases more slowly with t SHG R1 $ 1:1 ps and decays faster with t SHG D1 $ 1:6 ps. Because the optical excitation used (B10 20 cm À 3 ) 62 exceeds the natural n-doping level in the films (B10 19 cm À 3 ) 41 , the initial E dep dc , which originates from the upward band bending at the surface due to the higher density of free electrons residing in Dirac SS than those in the bulk states 59 , can be screened/enhanced depending on a balance of photoexcited electrons between the bulk states and Dirac SS. Specifically, due to the direct optical coupling of incident 1.51 eV photons to Dirac SS 69 , electrons below the Fermi level of occupied Dirac SS (1SS) can be photoexcited into the unoccupied Dirac SS (2SS) with an efficiency comparable to or even higher than that of the common valence-to-conduction band transitions in the bulk (Fig. 3d). The longer t SHG R1 compared to t TR R confirms that the screening/enhancement dynamics of E dep dc develop only after the electron-electron thermalization establishes a Fermi-Dirac distribution of photoexcited electrons. Due to an overlap in energy between Dirac 2SS and the high energy bulk bands, the thermalized electrons initially populate the upper cone of Dirac 2SS with high efficiency 48 , thus enhancing the initial E dep dc , which can be presented by solving the corresponding Poisson equation, where z SS and z dep are the widths of the Dirac SS actual range (B3 nm-a half of the critical film thickness for direct opposite-surface coupling) 51 and the surface depletion layer (B6.5 nm) 48,59 , respectively; n SS and n B are the densities of electrons residing in Dirac SS and the bulk states, respectively, and e is the effective low-frequency dielectric constant. The enhancement of E dep dc appears as an initial rise of the TSHG signal. The subsequent decay of the TSHG response can hence be associated with E dep dc screening due to the spatial redistribution of photoexcited electrons towards the bulk, which occurs primarily through electron-LO-phonon scattering. The resulting increase of the bulk electron density finally changes E dep dc sign (equation (5)), the process which is accompanied by the change of the initial upward band bending to the downward band bending that permits 2D electron gas (2DEG) to coexist with an electron population in Dirac 1SS (ref. 61). These dynamics appear as a drop of the TSHG signal to the minimal intensity at the dip which is controlled by the total density of photoexcited electrons dynamically residing in the bulk. This behaviour is consistent with the temporal threefold-tosixfold change of rotational symmetry in the TSHGRA patterns (weight coefficient c 2 in equation (2)) that closely follows the decay trend of the TSHG response (Fig. 3c). The minimal TSHG intensity, which can be characterized by the ratio of the SHG intensity at B7 ps (a dip) to that at B1.7 ps (TSHG-peak), progressively increases with decreasing d from 12 to 6 nm for the P in À S pump À P out polarization geometry, whereas it remains almost constant for the S in À P pump À S out polarization geometry ( Supplementary Fig. 3). This behaviour agrees with an increase of the electron-hole recombination rate in Dirac SS with decreasing d, because the density of electrons dynamically residing in the bulk is inversely proportional to the rate of carrier recombination in Dirac SS 48 .
After the TSHG intensity dips, it increases again with a second rise-time of t SHG R2 $ 5:8 ps. Taking into account the scattering processes between 2DEG and Dirac 1SS (ref. 61), the second rise of the TSHG signal and the corresponding second peak indicate the balancing of electron density between the 2DEG and Dirac 1SS by electron transfer from the former to the latter (Fig. 3d) 48,55,62 . This dynamical increase of electron density in Dirac 1SS progressively weakens E dep dc to 0, the process which is self-consistently controlled by the 2DEG-to-Dirac-1SS electron scattering rate. Correspondingly, the weight coefficient c 2 is stabilized at the minimal level, indicating that the purely surface second-order susceptibility tensor dominates the TSHG response. Subsequently, the TR and TSHG signals decrease with longer decay-times of t SHG D2 % t TR D2 $ 250 À 280 ps. Because both signals reveal similar relaxation dynamics at this stage, we associate this decay with carrier recombination in Dirac SS 48 .
The discussed ultrafast carrier SS-bulk-SS vertical transport and the resulting quasi-equilibrium between electron populations in 2DEG and Dirac 1SS also establishes the unique and stringent condition for the excitation of CADPs by the evanescent guided Lamb waves developed within the same relaxation process and for the fermion tunnelling between the opposite-surface Dirac SS through the region where the massless fermions acquire a finite mass and hence prevents CADPs from Landau damping similarly to the coupled bilayer structures 54  are almost comparable to the inverse of the repetition rate of the laser used (12.5 ns), thus allowing a quasi-steady accumulation of photoexcited electrons in Dirac 1SS as long as the sample is illuminated 48,53,55,62 .

Discussion
The observations presented in this article allow us to highlight several approaches which would be interesting to a broad audience of scholars exploring coherent low-energy collective electronic excitations in TIs. One of them is the application of nonlinear optical techniques allowing avoiding the photon-plasmon momentum mismatch for the excitation of surface plasmons. This approach seems to provide an advantage compared to the linear optical techniques, for which the waveguide evanescent Otto/Kretschmann prism couplers or grating structures deposited on the surface of TIs are required. In the latter case, which is widely used, the phase-matched plasmon frequency is limited only to the range determined by the micro-ribbon width and the grating period 46,47 . Moreover, the deposition of the micro-ribbon grating can hypothetically break the time-reversal symmetry of massless fermions in Dirac SS and hence strongly affect not only the intrinsic plasmon dynamics in the films, but also a unique metallic-type nature of Dirac SS, turning them into a trivial insulator state with massive fermions.
As we mentioned here, there are only a few reports discussing a time-domain detection of the slab thickness dependent oscillations which were associated with coherent APs 32,33 . However, these oscillations were detected using photocurrent autocorrelation and TR techniques, both being incapable, in principle, of distinguishing between the surface and bulk contributions and hence allowing for various interpretations. This statement is extremely important for TIs and Dirac systems in general, where the surface contribution is of major interest. Therefore, only applying surface and depletion electric field sensitive techniques, such as TSHG, the unambiguous time-domain detection of CADPs in TIs seems to be possible. Moreover, we have demonstrated here an advantage of simultaneous measurements of TR and TSHG responses, which due to their different origins provide us with more comprehensive information about carrier dynamics in TIs.
In this article we have also demonstrated a successful application of TSHGRA technique which allows for monitoring the crystal rotational symmetry with delay time. The observed spectacular change of threefold-to-sixfold rotational symmetry with delay time resulting from the dynamical screening of the surface depletion electric field in the TI Bi 2 Se 3 is unique and has never been reported previously. One of the attempts to use this technique was less successful due to extremely noisy SHG signals 68 .
Using TSHGRA technique, we have demonstrated that the dynamical elimination of the depletion-field-induced SHG response allows for observation of the exact sixfold rotational symmetry associated with the topmost hexagonally arranged Se-Se bonds in Bi 2 Se 3 films. The corresponding oscillations observed in this regime hence unambiguously characterize CADPs. The sensitivity of this technique exceeds that of the time-resolved and angle-resolved photoemission spectroscopy (TrARPES), which is known to monitor a few nanometres range in close proximity to the surface.
Another remarkable finding discussed in our article is the experimental proof of the validity of the theoretical prediction made by Das Sarma and Hwang 54 regarding plasmon excitation in coupled bilayer structures, which includes the modification of the AP dispersion when fermion tunnelling between the opposite-surface Dirac SS arises with decreasing film thickness.
The effect resonantly appears for the film thickness range limited by that associated with the direct coupling between the opposite-surface Dirac SS (6 nm) and that associated with the indirect coupling (15 nm). This resonance-type behaviour in Bi 2 Se 3 films centred at B10 nm is consistent with those appeared in the enhancement of high-frequency THz-range optical phonon modes 41 , nonlinear susceptibilities 53 , electron-phonon scattering strength 49 and electron-phonon relaxation dynamics 48,49,62 . All these observations point out that indirect opposite-surface coupling has a significant impact on carrier dynamics in this thickness range of the thin-film TI Bi 2 Se 3 .
Finally, we note that the CADP modes observed in Bi 2 Se 3 thin films using the TSHG technique could be expected across a wide range of 2D materials because the excitation of acoustic phonons and electron tunnelling are common stages of various electronic processes occurring in these materials.

Methods
Samples. Bi 2 Se 3 thin-films of 6,8,10,12,15,20,25,30,35 and 40 nm thick were grown on 0.5 mm Al 2 O 3 (0001) substrates by molecular beam epitaxy, with a 10 nm thick MgF 2 protecting capping layer, which was grown at room temperature without exposing the film to atmosphere. The quality of the films, their structure and thicknesses were determined from X-ray reflectivity, reflection high energy electron diffraction, and X-ray diffraction measurements 49 . The Hall conductivity measurements performed at 2 K showed that all films have n-doping level in the range of 0.5-3.5 Â 10 19 cm À 3 41 .
Experiments. The TR and TSHG measurements were performed using a multifunctional pump-probe setup that includes a Ti:Sapphire laser with an output average power of the fundamental laser beam of 2.5 W, pulse duration t L ¼ 100 fs, centre photon energy 1.51 eV and repetition rate 80 MHz 53 . The pump beam was at normal incidence and the probe beam was at an incident angle of B15°, focused through the same lens to a spot diameter of B100 mm. A pump beam average power was 640 mW. A probe beam average power was in the range of 100-580 mW. The TR and TSHG responses were measured simultaneously in the reflection geometry with various step sizes in delay times between the pump and probe pulses using a photodiode and photomultiplier tube, respectively. The signals were collected with a lock-in amplifier triggered at the pump beam modulation frequency of 800 Hz. The TSHGRA measurements were performed at various pump-to-probe delay times by rotating the sample mounted on a rotation stage about the surface normal with a step size of 2°. The linear polarization of the beams was either P (in the plane of incidence) or S (in the plane of the film). Four different light polarization geometries for the incident laser beam (in), the outgoing reflected fundamental and SHG beams (out), and the normal incidence pump beam were used to measure TSHGRA patterns: P in À S pump À P out , S in À S pump À P out , P in À P pump À S out , and S in À P pump À S out . Experiments were performed in air and at room temperature. No film damage was observed for the laser powers used in the measurements reported here.
Data availability. The data that support the findings discussed in this article are available from the corresponding author on request.