Mapping propagation of collective modes in Bi2Se3 and Bi2Te2.2Se0.8 topological insulators by near-field terahertz nanoscopy

Near-field microscopy discloses a peculiar potential to explore novel quantum state of matter at the nanoscale, providing an intriguing playground to investigate, locally, carrier dynamics or propagation of photoexcited modes as plasmons, phonons, plasmon-polaritons or phonon-polaritons. Here, we exploit a combination of hyperspectral time domain spectroscopy nano-imaging and detectorless scattering near-field optical microscopy, at multiple terahertz frequencies, to explore the rich physics of layered topological insulators as Bi2Se3 and Bi2Te2.2Se0.8, hyperbolic materials with topologically protected surface states. By mapping the near-field scattering signal from a set of thin flakes of Bi2Se3 and Bi2Te2.2Se0.8 of various thicknesses, we shed light on the nature of the collective modes dominating their optical response in the 2-3 THz range. We capture snapshots of the activation of transverse and longitudinal optical phonons and reveal the propagation of sub-diffractional hyperbolic phonon-polariton modes influenced by the Dirac plasmons arising from the topological surface states and of bulk plasmons, prospecting new research directions in plasmonics, tailored nanophotonics, spintronics and quantum technologies.

R hombohedral layered compounds of bismuth with group VI elements as Bi 2 Se 3 , Bi 2 Te 3 and their alloys have attracted renewed interest as prime candidates for the development of photonic devices exploiting their topological nature, 1,2 thanks to their thickness-dependent bandgap and to the electronic dispersion of their topological surface states (TSS) 3 . These topological insulators (TIs) 1 are semiconductors in the bulk, exhibiting metallic conduction at the surface, activated by the TSS, which form a Dirac cone around the Γ point of a hexagonal Brillouin zone 2 .
TIs provide a versatile platform for investigating quantum phenomena and for applications in electronics and photonics: with almost the same high-absorbance of graphene 4 , they can exploit a tunable surface bandgap 5 ; their optical and electronic properties can be engineered by changing the material stoichiometry 6 ; and the two-dimensional electron gas (2DEG) arising from the TSS can support collective excitation (Dirac plasmon) at terahertz (THz) frequencies 7,8 . This opens intriguing perspectives for a variety of application fields, such as quantum computing 9 , spintronics 10 , and photodetection 6 .
Some of these applications may benefit from the exploitation of electronic collective excitations, like plasmons, 11,12 plasmonpolaritons 13 -hybrid light-matter modes involving the collective oscillations of charge carriers, phonons 14 or phonon-polaritons 15,16 , whose wavelength can be electrically controlled. However, capturing these modes on the TSS of a TI or even in the bulk is a very demanding task: it requires sophisticated optical techniques capable of imaging at their exact frequency or mapping their long-range propagation; very importantly, the correct identification of TSS requires the knowledge of their complete energy/momentum dispersion, that can be only unveiled in TI samples showing sufficiently high structural quality (both stoichiometric and crystallographic) 17,18 .
Amplitude-and phase-resolved scattering near-field optical microscopy (s-SNOM) 14,15,[19][20][21][22][23][24][25] can allow access to the spatial variation of complex-valued dielectric response of layered 2D materials, heterostructures and low dimensional systems, therefore providing the sub-diffractional spatial resolution required to investigate light-matter interaction at the nanoscale. Moreover, s-SNOM is the ideal tool to investigate polaritons since its sharp tip acts both as a launcher and as a detector of propagating collective modes 15 . A large variety of interferometric approaches has been developed in the last years, allowing amplitude and phase-resolved s-SNOM imaging in the visible and infrared spectral ranges 19,20 , including fs pulsed laser sources and electro-optic sampling detection in the THz range, 21 microwave circuitry in the sub-THz range 22,23 or interferometric techniques 24,25 strongly limited by the poor dynamic range of the cryogenically cooled bolometric detectors needed to measure the typically small s-SNOM THz signals.
At THz frequencies, the dielectric response of TIs based on bismuth chalcogenides is characterized by longitudinal (LO) and transverse (TO) optical phonons that identify frequency regions of hyperbolicity, within which the dielectric permittivities along the inplane and out-of-plane directions have an opposite sign and the isofrequency surfaces of the extraordinary rays are hyperboloids in the momentum space 26,27 . Hyperbolic phonon-polaritons modes have been intensively investigated in hexagonal boron nitride (hBN) 15,28,29 , where they occur at mid-infrared frequencies, due to their unique capability to support the propagation of light with momentum far exceeding the free-space value without evanescent decay 15,28 . Remarkably, in TIs, these polaritons, that have great potential for THz waveguiding, can interact with the electrons associated to the TSS resulting in doping tunable dispersion 30 .
In this letter, we exploit an innovative combination of hyperspectral time-domain spectroscopy (TDS) nano-imaging 31 and detectorless scattering near-field optical microscopy (s-SNOM) 14 , at three distinctive pumping energies in the far-infrared, to explore the rich physics of thin flakes of Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 , having different thicknesses. We select Bi 2 Se 3 because of its simple electronic structure, showing a topologically nontrivial direct energy gap of 0.3 eV 32,33 , about ten times higher than thermal excitation at RT 34 and higher than that of Bi 2 Te 3 , having an energy gap of about 0.15 eV 35 . We then investigate Bi 2 Te (3−3x) Se 3x because, on one side, the energy gap is expected to increase rapidly with x, reaching values >0.2 eV for 3x = 0.8 35 , and on the other side, the Fermi level of high Te-content alloys, like Bi 2 Te 2.2 Se 0.8 , crosses only the TSS Dirac cone, meaning that its TSS electrons behave as an ideal 2D electron gas 36 .
Theoretical proposals 2 on TSS require the chemical potential to lie at, or to approach, the surface Dirac point. However, the surface Fermi level of a TI depends on the detailed electrostatics of the surface and it is not necessarily at the Dirac point. In a naturally grown Bi 2 Se 3 , the bulk Fermi energy does not lie in the gap 33 . The material displays n-type behavior and usually exhibits significant contributions from the bulk carriers to the total conductivity 33,34,[37][38][39] . One reason for the shift of the chemical potential into the conduction band is donor doping, typically induced by Se vacancies 2 and antisite defects 40 . Another effect, which can interfere with the ideal TI response, occurs when the bulk bands cross the Fermi level near the material surface, generating a 2DEG of massive carriers close to the surface if the bands bend downwards 41 . This has been observed in similarly grown single-crystal Bi 2 Se 3 when exposed to the atmosphere after cleaving 42 , as a consequence of extrinsic defects or molecular adsorption at the surface 36 . The 2DEG is believed to be confined within 20 nm from the surface in Bi 2 Se 3 36 and coexists with the TSS. Engineering the Bi 2 Te (3−3x) Se 3x stoichiometry can therefore allow tuning the TI physics, on purpose.
Our combination of two different near-field techniques allows identifying optically activated TO and LO phonons, and reveals massive bulk plasmons, arising from the 2DEG originated from the band bending at the sample surfaces in Bi 2 Se 3 , and hybrid plasmon-phonon-polariton modes propagating from the edge of the Bi 2 Te 2.2 Se 0.8 flakes. Remarkably, this represents the first experimental demonstration of hybrid plasmon-phonon-polariton modes reported to date in topological insulators and paves the way for novel subdiffractional and tunable THz polaritonic waveguides.

Results and discussion
Samples. Flakes of thickness d in the range 15-200 nm and lateral size in the range 0.2-4 μm are obtained by mechanical exfoliation on a 350-μm-thick silicon wafer, with a 300 nm thin insulating SiO 2 toplayer. The micro Raman characterization of prototypical flakes of Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 is reported in Supplementary Fig. 1a, b, respectively. Bi 2 Se 3 and Bi 2 Te 3 crystallize in tetradymites with rhombohedral crystal structure endowed with D 5 3d space group symmetry and five atoms per unit cell. The structure consists of quintuple layers (QL) stacked via weak van der Waals forces, which permit an easy mechanical exfoliation, orthogonally to the trigonal c-axis. One QL corresponds to a 1.42 nm thick crystal of Bi 2 Se 3 and to 1 nm for Bi 2 Te 2.2 Se 0.8 6 . Within each QL, pictogen (Bi) and chalcogen atoms (Se, Te) are held together by ionic covalent bonds. The different bond strength and character result in a strong in-plane/outof-plane anisotropy in the thermal, electronic and dielectric properties of these materials.
Time-domain near-field spectroscopy: the role of hyperbolic phonons. Hyperspectral THz nano-imaging 31 is performed by coupling the s-SNOM microscope to a THz TDS system based on photoconductive (PC) antennas (Menlo TERA15 FC), as schematically depicted in Fig. 1a.
The broadband THz pulse produced by the transmitter antenna, with a spectrum spanning the 4-140 cm −1 range, is focused on the atomic force microscope (AFM) tip of the s-SNOM and the forwardscattered field is detected as a function of the time delay from the gate pulse by a second PC antenna, the receiver, working in reverse mode. This technique combines broadband spectral coverage with subdiffractional spatial resolution of~260 nm (see Supplementary Fig. 4) achieved by using an AFM tip with apex radius r = 40 nm, 80 μm cantilever length (Rocky Mountain Nanotechnology, 25PtIr300B) working in tapping mode, with an oscillation amplitude A~200 nm and a tapping frequency Ω~20 kHz.
The detected THz scattered field S n = s n e −iϕ , with amplitude s n and phase ϕ, is demodulated at the tapping frequency nth harmonics Ω n = nΩ with n = 2-4 to isolate the near-field contribution from the background due to far-field illumination of the tip shaft and sample 19 . Images of the near-field scattered amplitude s n , at a fixed time delay, for thick flakes (d > 75 nm) of Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 as a function of the demodulation order n, are reported in Supplementary  Fig. 2a, b. To find a good compromise between signal intensity and suppression of the far-field background, we investigate the second harmonic near-field signal. This choice is supported by the analysis of the approach curves ( Supplementary Fig. 3) at the different demodulation orders.
Nano-THz spectra at fixed positions on thick (d > 75 nm) Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 flakes and on a reference point on an Au marker are reported in Fig. 1b. Each spectrum is obtained by averaging 60 different spectra acquired with 100 ps time-scan, lasting 45 s, performed in a nitrogen atmosphere to prevent absorption by water vapor 43 . Due to the frequency-dependent scattering efficiency of the tip 44 , a signal-to-noise ratio (SNR) > 1.5 is achieved only in the range 15-90 cm −1 . In this range, the two crystals exhibit distinct contrast η Au , evaluated by the ratio of the second-order demodulated spectra s 2 scattered by the flake with that of the Au reference, as routinely done in nano-FTIR spectroscopy 45 . While both compounds show a broad peak in the probed range, the maximum η Au is reached at~80 cm −1 for Bi 2 Se 3 and at~60 cm −1 for Bi 2 Te 2.2 Se 0.8 . This scattering enhancement can be attributed to infrared-active bulk optical phonons [46][47][48] . The peak in the near-field scattering intensity of Bi 2 Se 3 occurs within the frequency band ω to ⊥ < ω < ω to // identified by the dominant long-wavelength TO optical phonon modes 30 : the E u TO (ω to ⊥ = 64 cm −1 ) and the A 2u TO (ω to // = 135 cm −1 ) phonons, involving atomic vibrations in the plane orthogonal (⊥) and parallel (//) to the trigonal c-axis, respectively. Within this frequency range, the real parts of the permittivity along the directions orthogonal ε ⊥ and parallel ε // to the c-axis have an opposite sign (Re{ε ⊥ } < 0, Re{ε // }>0) and Bi 2 Se 3 is expected to behave as a hyperbolic material of type II. We can use the vibrational modes of Bi 2 Te 3 to model the Bi 2 Te 2.2 Se 0.8 flake, due to the low dispersion 46 of the phonon modes of the Bi 2 Te (3−3x) Se 3x alloy for x < 1/3. The lower frequency of the Bi 2 Te 2.2 Se 0.8 peak can be traced back into the red-shift of E u and A 2u TO phonons with Three-dimensional model calculations of the near-field contrast η based on the well-established finite-dipole approximation 49,50 , considering multiple LO and TO phonon modes, both in // and ⊥ directions, are reported in the Supplementary Fig. 5a. The model well reproduces the experimental spectra, suggesting that Bi 2 Te 3 can reasonably mimic the near-field response of the Bi 2 Te 2.2 Se 0.8 . Moreover, since the exposed surfaces are parallel to the basal plane, the ability of detecting ⊥ modes implies that the component of the THz near-field orthogonal to the tip-tapping is non-negligible, as previously observed in anisotropic crystals as SiC 49 .
In order to disentangle the dielectric response from the bulk and the surface, we then combine the finite-dipole approach to describe the tip-sample interaction, with a multilayer model for up to four layers and we implement an inversion algorithm to extract the local dielectric function of the surface layer alone 51 (see Supplementary Note 3 and Fig. 6). We note that even if we include a finite bulk carrier density, both Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 maintain their hyperbolic behavior as testified by the opposite signs of the real parts of the bulk dielectric functions in the // and ⊥ directions (see Supplementary Fig. 6d, e). The comparison between the two models (the multilayer model in Fig. 1e and the bulk model in Supplementary Fig. 5a) clearly unveils that while in Bi 2 Se 3 the amplitude contrast is well reproduced by the bulk model and the TSS do not appear to contribute to the near-field response (black traces in Fig. 1e and Supplementary Fig. 5a), to describe the low frequency (ω < 50 cm −1 ) amplitude contrast measured in Bi 2 Te 2.2 Se 0.8 , we have to include a finite contribution from TSS to the near-field amplitude (see red traces in Fig. 1e and Supplementary Fig. 5a).
Detectorless near-field nanoscopy: role of the topological surface states. With the aim of isolating the nature of the modes activated through THz photoexcitation and their interaction with TSS, we then perform self-detection nanoscopy on a set of flakes with a progressively reduced thickness, employing single-frequency mode THz-quantum cascade lasers (QCLs), 14   and high spectral purity, with intrinsic linewidths as low as 100 Hz 52 ; this allows for selective interaction with THz resonant modes. In our self-detection experiment 14 , the QCL acts simultaneously as a powerful THz source and as a phase-sensitive detector, employing a detection scheme based on self-mixing (SM) interferometry 53 , therefore overcoming the need for bulky cryogenic detectors. This effect is based on the reinjection of a small fraction (10 −4 −10 −2 ) of the emitted field that coherently interferes within the laser cavity. The fraction of the laser radiation backscattered from the tip-sample system can be coupled back into the QCL cavity along the same optical path. The coherent superposition of reinjected THz field with the QCL intracavity field produces a perturbation of the laser voltage ΔV that depends on both the amplitude and the phase of the THz field scattered by the tip.
s-SNOM experiments allow to excite collective modes of finite momentum by measuring the local response, which is given by the superposition of different modes with a distribution of momenta extended up to q $ 1/r, rather than by a single-mode. Exploiting the higher power and sensitivity of the self-detection experiment with respect to the TDS system, we use a tip with a smaller radius, r = 10 nm, corresponding to a spatial resolution of about 30 nm (see Supplementary Fig. 8) and to a larger probed momentum range.
In order to sample the SM fringes and retrieve the phase shift experienced by the field in the scattering process, the length L of the external cavity, formed by the tip and output laser facet, is varied up to 350 μm (>2λ = 4πc/ω 0 ) with a delay line built on a linear translation stage with sub-μm spatial resolution (Fig. 2a).
A variable attenuator with transmission t is inserted in the optical path to reduce the optical feedback to the weak limit in which ΔV has a simple sinusoidal dependence on L, expressed, within the framework of the Lang-Kobayashi model 54 , as ΔV $ scosð 2ω 0 L c À ϕÞ, with s signal amplitude and ϕ phase of the scattered field. Lock-in detection of the nth harmonics of the tapping frequency Ω is then used to isolate the near-field component from the scattering background. Using a lock-in detection technique we have a noise level of 0.1 μV and a SNR > 40 on the n = 5 harmonics on thick flakes.
The scattering signal at the highest harmonic order (n = 5), acquired keeping fixed the cavity length L, is enhanced at the flakes compared to the substrate, see Fig. 2d, e. For thicker flakes with d > 40 nm, the amplitude of s 5 appears higher for Bi 2 Te 2.2 Se 0.8 (B1-B2 in Fig. 2c) than for Bi 2 Se 3 (A1-A2 in Fig. 2b) in agreement with the η Au values extracted via hyperspectral THz nano-imaging (Fig. 1c). However, the two materials show a different dependence on thickness: while negligible amplitude variation is observed upon reducing the thickness of Bi 2 Se 3 from 86 nm (flake A1) to 16 nm (A4) (see Fig. 2d), the s 5 amplitude drops in Bi 2 Te 2.2 Se 0.8 such that the 19 nm thick flake (B4) is barely distinguishable from the substrate (Fig. 2e).
We then extract line profiles from the maps in Fig. 2b, e while moving along a line orthogonal to the flakes edge and keeping fixed the cavity length L, to show the variation of the substrate/flake contrast with thickness in the two material systems (Fig. 2f, g).
The thickness-dependent signal change, revealed in the line profiles (Fig. 2f, g), could be caused by both a phase shift and by an actual variation of the amplitude of the SM fringes. To discriminate between the two effects, we resolve the phase of the scattered field (predicted in Fig. 1d): we acquire the SM fringes with 0.2 μm steps while scanning the sample along a line orthogonal to the substrate-flake edge. Figure 3a shows the retrieved amplitude signal at the third demodulation order s 3 , as a function of L and of the position on the sample labeled as X.
Considering the drop in the scattering efficiency of the tip with increasing THz frequencies, we choose the third-order demodulated signal s 3 to compare the data collected at the three different frequencies ω 0 . By analyzing the SM fringes at each X with sinusoidal fit functions (as in Fig. 3b), we extract the amplitude s 3 and the phase ϕ 3 along the scan trajectory (see Supplementary Fig. 9 for the phase profiles).
Since the absolute value of the SM signal depends on the threshold voltage of the THz-QCL, on its output power, and on the frequency-dependent scattering efficiency of the tip, the average scattering signal amplitude from the flake s 3 (flake), is divided by the signal from the substrate, s 3 (substrate), whose response is expected to show negligible variation with frequency. The refractive index variation of SiO 2 in the 60-90 cm −1 range is indeed <0.5% 55 . The thickness dependencies of the ratio η 3, which is accounting for the contrast between flake and substrate, and of the phase change Δϕ 3 , at the three pumping frequencies, are reported in Fig. 3c, e, and in Fig. 3d, f, for Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 , respectively. The near-field response appears very different in the two materials.
In Bi 2 Se 3 , the flat topography regions show roughly constant contrast with the SiO 2 substrate changing the flake thickness d, at all the probing frequencies, while only the signal at the flake edges varies with thickness. Looking at s 3 (ω 0 ), by keeping fixed the layer thickness, the contrast η 3 increases by 30-40% while increasing ω 0 from ω 0 = 66.7 cm −1 to ω 0 = 76.7 cm −1 in Bi 2 Se 3 (Fig. 3c). This is in disagreement with the calculated 2D near-field response of Bi 2 Se 3 that predicts 30 a peak in the s 3 amplitude at the TO phonon wavenumber ω to ⊥ = 64 cm −1 and a strong decrease of about 40% in the range 65-90 cm −1 . On the other hand, the observed increase of s 3 is compatible with the overall increase in the same range shown by the calculated contrast in Fig. 1e, and hence suggests a bulk phononlike nature for the dielectric response of our Bi 2 Se 3 samples, in agreement with what emerged from previous reports 56 . On the contrary, Bi 2 Te 2.2 Se 0.8 exhibits a visible decrease in the scattering intensity with flake thickness d at ω 0 = 66.7 cm −1 and ω 0 = 76.7 cm −1 , which becomes less pronounced at ω 0 = 89.7 cm −1 , such that, in all cases, the flakes become barely distinguishable (η 31 ) from the substrate when d < 20 nm (Fig. 3d). The thickness and frequency dependences of s 3 from Bi 2 Te 2.2 Se 0.8 are consistent with a recent model 30 of the near-field response formulated for Bi 2 Se 3 which predicts deeply sub-diffractional, highly directional hyperbolic phonon-polaritons interacting with the electrons of the TSS. In particular, the experimental frequency dependence of Bi 2 Te 2.2 Se 0.8 is retrieved if one includes a red-shift of the spectral features predicted by the model for Bi 2 Se 3 towards the optical phonon frequencies of Bi 2 Te 3 to describe Bi 2 Te 2.2 Se 0.8 , whose phonon frequencies reasonably approaches the ones of Bi 2 Te 3 46 . A detailed simulation of the near-field response s 3 of Bi 2 Te 3, as a function of the pumping frequency and of the flake thickness, is reported in the Supplementary Information (Supplementary Note 4 and Supplementary  Fig. 15a). Consistently, in the Bi 2 Te 2.2 Se 0.8 flakes (Fig. 3d), s 3 (ω 0 ) is almost frequency independent in the 66.7 cm −1 < ω 0 < 76.7 cm −1 range and then increases when ω 0 = 89.7 cm −1 , as we approach the peak of s 3 due to the TO phonon at ω to // = 95 cm −1 : A strong and sharp peak in near-field contrast is predicted at ω to ⊥ 30 , that in a Bi 2 Te 2.2 Se 0.8 falls outside of the investigated spectral range (ω to ⊥ = 50 cm −1 for the Bi 2 Te 3 ) 46,47 . However, differently from the peak at ω to // , this peak is related to the farfield factor rather than to the near-field tip-sample interaction.
Propagating collective excitations. Figure 4 shows the near-field maps (Fig. 4a, b) of the edges of the Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 flakes, acquired at ω 0 = 66.7 cm −1 for a fixed cavity length L, and the related line profiles extracted by averaging over the displayed window along the vertical axis (Fig. 4c, d).
The near-field scattering maps (Fig. 4a, b) show a visible signal enhancement at the flake edges with the brightest edge corresponding to the side from which the THz beam enters in the s-SNOM.
A strong signal enhancement at the flakes edges has been previously identified in hBN as a signature of the activation of localized phonon-polariton modes 57 . The peaks localized at the flake edges have distinct thickness dependence in the two materials: in Bi 2 Se 3 the peak broadens and increases its intensity with the thickness (Fig. 4c), while it appears always sharper, and with a less pronounced dependence of intensity on thickness in the Bi 2 Te 2.2 Se 0.8 (Fig. 4d). We ascribe those peaks to edge resonances of phonon-polariton modes 57 , whose interplay with the main photoexcited modes quickly damps them close to the edge. Indeed, while the peaks retrieved in Fig. 4c, d are visibly edge modes, we observe signal intensity modulations in regions on flat topography (Fig. 4c, d), characterized by a much weaker amplitude, which we attribute to the interference of propagating modes launched by the tip and reflected at the flake edges. In our case, these oscillations in Bi 2 Se 3 have comparable amplitude changing the thickness, while in Bi 2 Te 2.2 Se 0.8 they display a less pronounced intensity for thinner flakes, becoming not distinguishable in the thinner flake (d = 19 nm).
To unveil the nature of the aforementioned propagating photoexcited modes, we analyze these signal oscillations with the function:  Fig. 4c, d), 1= ffiffiffiffiffiffiffiffiffiffiffiffi ffi x À x 0 p is the geometrical decay and the exponential describes oscillations with periodicity π q p , where the wavevector q p is complex-valued q p = q 1 + iq 2 . The double period λ p = 2π q 1 represents the wavelength of the photoexcited modes, whose frequency is set by the energy conservation law to match that of the impinging photon 59 . Figure 5a shows the obtained λ p as a function of the flake thickness for Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 , respectively, at the three probed frequencies. In Bi 2 Se 3 , λ p is almost thickness independent (Fig. 5a) at all probing frequencies. Vice versa in Bi 2 Te 2.2 Se 0.8 , there is a visible dependence of the oscillation period on d (Fig. 5b): λ p decreases while reducing d.
The thickness-independent trend, extrapolated in Bi 2 Se 3, is compatible with photoexcited massive plasmons formed by the bulk charge carriers at the flakes surface, which give rise to a 2DEG, as a consequence of the band bending at the Bi 2 Se 3 interface with adjacent materials 36 .
Conversely, the extrapolated thickness dependence of λ p in Bi 2 Te 2.2 Se 0.8 is expected either in presence of massive plasmons related to bulk charge carriers or for Dirac plasmons related to TSS 56,60 . It is worth mentioning that, in the limit of small wavevectors (qd < 1), Dirac plasmons dispersion can be approximated to the thickness-independent first-order dispersion 8 , while at larger wavevectors both the dependence on the thickness and on the bulk dielectric function emerge 60,61 .
We then investigate the experimental mode dispersion (Fig. 5c, d), for both material systems, and compare it with the theoretical trends predicted under the hypothesis discussed above. A quadratic dependence of the real part of the wavevector q 1 = 2π/λ p on the frequency ω p =ω 0 is observed for both materials. In Bi 2 Se 3 the dispersion curve, calculated for massive plasmons formed by the 2DEG charge carriers 36 , perfectly matches our experimental data (Fig. 5c), confirming our claim. The dispersion of massive plasmons 30 is evaluated using the relation ω 2 p ¼ e 2 2ε 0 ε r n M m eff q 1 , where ε r~6 .3 is the average permittivity of the materials surrounding the flake (silica/ silicon substrate and air), e is the electron charge, ε 0 is the vacuum permittivity, n M is the carrier density associated with the 2DEG and m eff = 0.15 m 0 is the effective mass 34 of Bi 2 Se 3 , with m 0 electron mass. Good agreement with the data is obtained considering a surface carrier density n M = 3 × 10 10 cm −2 , which indicates charge-carrier depletion at the surface, in agreement with previous studies 62, 63 on Bi 2 Se 3 .
The theoretical dispersion curves predicted in Bi 2 Te 2.2 Se 0.8 for massive bulk plasmons and for TSS Dirac plasmons are shown in Supplementary Fig. 12a, b, respectively. In the case of massive bulk plasmons, we use the same dispersion relation as for Bi 2 Se 3 massive plasmon but with the effective mass 64 of Bi 2 Te 3 m eff = 0.044 m 0 , and a thickness-dependent density for the bulk massive electrons n M (d). The thickness dependence is here encoded into the carrier density term n M = n s + C· n bv d 0.5 , which depends on the sheet carrier density n s , on the bulk carrier density n bv , and for which we consider a square-root dependence on the thickness d as in ref. 65 . A bulk carrier density n bv = 3 × 10 18 cm −3 in Bi 2 Te 2.2 Se 0.8 is taken from the value experimentally retrieved with direct Hall measurement in Bi 2 Te 2.4 Se 0. 6 66 , the coefficient C = 0.8 × 10 −5 m 0.5 is taken from ref. 65 and estimation of n s = 2 × 10 12 cm −2 is obtained by multiplying n bv from ref. 66 for the thickness of the flake used for the n bv determination (6 QL~8 nm). On the contrary, the dispersion of Dirac plasmons associated with TSS 56,60 is: where n D is the Dirac carrier density (assumed equal to n s to calculate the dispersion), h is the Planck constant, and ε TI (ω) is the Bi 2 Te 3 permittivity 56 . The frequency dependence of ε TI is crucial to account for the thickness dependence of the Dirac plasma frequency 67 , while here we use an average of the dielectric constants along the basal plane ε ⊥ and along the c-axis ε // , computed as: Our model for the dielectric permittivity of Bi 2 Se 3 and Bi 2 Te 3 is detailed in Supplementary Note 4. The data of Supplementary Fig. 12a, b indeed clearly show that for Bi 2 Te 3 none of these trends matches with our experimental data (Fig. 5d), and that, for both assumptions, we would expect much lower q values. This reflects in propagating plasmon wavelengths λ p (Supplementary Fig. 12c, d) significantly longer (>100 μm), impossible to be captured via near-field nanoscopy.
We cannot, therefore, ascribe the oscillations observed in Bi 2 Te 2.2 Se 0.8 to propagating bulk plasmons or Dirac plasmons. To shed further light on the nature of these modes we then perform numerical simulations of the collective mode energy dispersion in our Bi 2 Te 2.2 Se 0.8 flakes, following the procedure described in ref. 30 . We look for the dispersion of phonon-polariton modes in the presence of charge carriers populating the TSS, described by the surface sheet conductivity of Dirac fermions 68 , assuming that the chemical potential is located in the bulk bandgap, resulting in a negligible bulk conductance. The charge carriers in the TSS are expected to form Dirac plasmons that hybridize with the phononpolariton modes to form hyperbolic plasmon-phonon-polaritons.
The near-field signal can be considered as a weighted average of the surface reflectivity for p-(or TM-) polarized light r P (q, ω 0 ) over q. The dispersion of collective modes is identified by the maxima, at real arguments q and ω, of the imaginary part of the r P (q, ω), Imfr P ðq; ωÞg, which we calculate including the strong anisotropy of the frequency-dependent dielectric permittivity of Bi 2 Te 2.2 Se 0.8 (see Methods).
The false color maps of the function Im{r P (q, ω)} in Fig. 6 provide a convenient visualization of the collective mode spectra.
The bright lines in Fig. 6a-c are the dispersion curves of the collective modes calculated for Bi 2 Te 2.2 Se 0.8 flakes of different thickness, assuming a constant chemical potential μ = 200 meV, which is chosen to fall within the electronic bandgap of the material (~280 meV 35 ). The Im{r p } map calculated for chemical potentials μ = 50 meV, μ = 150 meV, and μ = 250 meV and for flake thickness 118 nm is reported in the Supplementary Fig. 14.
In bismuth chalcogenide compounds, the electrochemical potential is indeed expected to approach the bottom of the conduction band as a consequence of the doping due to vacancies or antisite defects 69 . The smearing of the dispersion lines gives an idea of how damped the modes are.
The comparison with our experimental data (white crosses in Fig. 6a-c), show that our photoexcited modes all falls in the upper half of the hyperbolic band ω to ⊥ < ω < ω to // and their q values correspond to those expected for hybrid modes resulting from hybridization of hyperbolic phonon-polaritons with TSS Dirac plasmons into combined hyperbolic plasmon-phonon-polaritons modes (HP 3 modes), and differ significantly from those expected when only Dirac plasmons come into play (see Supplementary  Fig. 12b). The effect of the hybridization is encoded in the frequency shift of the dispersion curves when there is finite doping (see Supplementary Fig. 14). The damping of these modes is a consequence of the strong electron interaction with the charge carriers and increases with doping. From Fig. 6 we can also observe that the hyperbolic phonon-polariton modes fade away while decreasing the flake thickness. It is worth mentioning that the carrier density plays a fundamental role in dictating the observable activated excitations in our geometry, i.e., with wavelengths within the physical size of the exfoliated flakes. Indeed, in the case of Bi 2 Se 3 , the expected lower surface carrier density, combined with the larger effective mass m eff , entails the activation of short wavelength λ p (submicron) massive plasmons, whereas these are not observable in Bi 2 Te 2.2 Se 0.8 because the higher carrier density and lower m eff push λ p in the 0.5-1 mm range. Therefore, the ability to controllably tune the carrier density, e.g., by electrostatic gating, could open the possibility to selectively ignite the different phenomena.  Supplementary Fig. 12 and therein-discussed in detail. Error bars indicate the 95% confidence bands from the fit of interference patterns. In the plots, the Fermi velocity is v F = 0.5 × 10 8 cm/s as for Bi 2 Te 3 from ref. 48 , the electron scattering rate is γ e = 1 THz as in ref. 4 , the vacuum and substrate permittivity are ε 0 = 1, and ε sub = 6.3. The white symbols correspond to the experimental data of Fig. 5d.
Analogous activation of mixed plasmon and hyperbolic phonon-polaritons has been predicted for hBN thin slabs 30 , with the difference that the hyperbolic bands of hBN fall in the midinfrared range, at much higher frequencies (around 1500 cm −1 ). Technological applications based on mid-infrared hyperbolic phonon-polaritons demonstrated for hBN, such as hyperlensing 70 , waveguiding 30 , and light focusing 71 , can be translated to the THz frequency range using Bi 2 Te 2.2 Se 0.8 flakes.
The hybridization of phonon-polaritons with Dirac plasmons introduces additional control knobs, namely flake thickness and doping, to tune the HP 3 mode dispersion and accordingly the THz response of Bi 2 Te 2.2 Se 0.8 flakes. Moreover, since HP 3 modes are highly directional 30 , a change in the doping can be exploited to vary the propagation direction of the HP 3 .
In conclusion, we have provided the first direct experimental evidence of hyperbolic plasmon-phonon-polaritons modes at THz frequencies, and of the activation of massive bulk plasmons associated to band bending induced 2DEG, providing novel insights on the rich physics of TI materials. The interaction of the probed phonon-polaritons modes with the electrons from the TSS suggests possible ultrafast optical control 16 by above bandgap photoexcitation to ideally switch on and off the hybridized modes, opening intriguing technological perspectives in plasmonics 72 , ultrafast photonics, spintronics 10 , nanophotonics 73,74 , and quantum optics 75 .

Methods
Growth. Single-crystalline ingots of Bi 2 Se 3 and Bi 2 Te 2.2 Se 0.8 were grown from melt by the vertical vacuum/inert atmosphere optimized Bridgman−Stockbarger method, following methods described in ref. 76 . The resulting samples show extremely high chemical stability and do not manifest any effect related to degradation under air exposure.
Total reflectivity of Bi 2 Te 2.2 Se 0.8 . The sample system is modeled as a stack of three materials: the air superstrate (labeled as 0) with permittivity ε 0 = 1, from which the incident THz beam comes, the Bi 2 Te 2.2 Se 0.8 flake of thickness d (labeled as 1) with sheet conductivity, σ SS ðq; ωÞand the undoped silicon/silica substrate (labeled as 2) with permittivity ε s = 6.3. Following ref. 30 , the P-polarization reflectivity r P of this layered system is given by: r P ¼ r 12 r 01 þ r 10 À 1 À Á À r 01 e ðÀ2ik == dÞ r 10 r 12 À e ðÀ2ik == dÞ ð1Þ where k // = ffiffiffiffiffi ε ? p ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω 2 c 2 À q 2 ε == q is the out-of-plane momentum inside the flake, and ε ⊥ and ε // are the in-plane and out-of-plane dielectric permittivities calculated from the parameters available in the literature for Bi 2 Te 3 as detailed in the supplementary information. The term r 12 is the reflectivity at the flake/substrate interface and r 01 , r 10 are the reflectivities at the air/flake interface as coming from the air and from the flake side, respectively. These latter are computed as follows in cgs units: k s and k 0 are the out-of-plane momenta at the substrate and in the air and are defined in analogy to that in the flake but considering that they are isotropic materials: . Finally, the sheet conductivity σ SS is evaluated as: with P polarizability of Dirac fermions 30,68 , in which the Fermi velocity of the material v F = 0.5 × 10 8 cm/s and the chemical potential μ enter. By using the same σ SS in r 01 and in r 12 we are assuming the same sheet conductivity at both the interfaces air/flake and flake/substrate. The sheet conductivity is nonzero only for finite doping and determines the appearance of the Dirac plasmon modes and their hybridization with hyperbolic phonon-polaritons to form HP 3 excitations.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding authors upon reasonable request.

Code availability
The codes and simulation files that support the plots and data analysis within this paper are available from the corresponding author upon reasonable request.