Tielrooij et al. Light: Science & Applications (2022)11:315 Official journal of the CIOMP 2047-7538 https://doi.org/10.1038/s41377-022-01008-y www.nature.com/lsa ARTICLE Open Access Milliwatt terahertz harmonic generation from topological insulator metamaterials Klaas-Jan Tielrooij 1,2 ✉, Alessandro Principi3, David Saleta Reig 1, Alexander Block 1, Sebin Varghese 1 , Steffen Schreyeck4, Karl Brunner4, Grzegorz Karczewski4,5, Igor Ilyakov6, Oleksiy Ponomaryov6, Thales V. A. G. de Oliveira 6, Min Chen 6, Jan-Christoph Deinert 6, Carmen Gomez Carbonell1, Sergio O. Valenzuela 1,7, Laurens W. Molenkamp4,8, Tobias Kiessling4, Georgy V. Astakhov 6 ✉ and Sergey Kovalev 6 ✉ Abstract Achieving efficient, high-power harmonic generation in the terahertz spectral domain has technological applications, for example, in sixth generation (6G) communication networks. Massless Dirac fermions possess extremely large terahertz nonlinear susceptibilities and harmonic conversion efficiencies. However, the observed maximum generated harmonic power is limited, because of saturation effects at increasing incident powers, as shown recently for graphene. Here, we demonstrate room-temperature terahertz harmonic generation in a Bi2Se3 topological insulator and topological-insulator-grating metamaterial structures with surface-selective terahertz field enhancement. We obtain a third-harmonic power approaching the milliwatt range for an incident power of 75 mW—an improvement by two orders of magnitude compared to a benchmarked graphene sample. We establish a framework in which this 1234567890():,; 1234567890():,; 1234567890():,; 1234567890():,; exceptional performance is the result of thermodynamic harmonic generation by the massless topological surface states, benefiting from ultrafast dissipation of electronic heat via surface-bulk Coulomb interactions. These results are an important step towards on-chip terahertz (opto)electronic applications. Introduction materials with massless Dirac fermions— charge carriers One of the most canonical manifestations of the non- with a linear energy-momentum dispersion relation—can linear interaction of electromagnetic waves with matter is lead to efficient THz harmonic generation. This has been harmonic generation, where the energy of generated shown with the prototypical massless Dirac-fermion photons is a multiple of the energy of incident photons. material graphene6, with topological insulators (TIs)7,8, Until a few years ago, the most efficient way of generating and with Dirac semi-metals9,10. Currently, the THz non- harmonics in the terahertz (THz) regime was using two- linear photonics “toolbox” includes the possibility to dimensional quantum well systems1–3. During the past electrically tune harmonic generation11, and the ability to few years several other material systems, including employ ultrathin metamaterials that generate an SrTiO34 and superconductors5 have shown strong THz increased harmonic power by orders of magnitude, as a nonlinearities. In particular, it has been demonstrated that result of field enhancement in a metallic grating located nearby the nonlinear Dirac-fermion system12. This has resulted in a record-high third-order susceptibility Correspondence: Klaas-Jan Tielrooij (klaas.tielrooij@icn2.cat) or χ(3) above 10–8 m2 V–2 in the THz regime for a grating- Georgy V. Astakhov (g.astakhov@hzdr.de) or Sergey Kovalev (s.kovalev@hzdr.de) graphene metamaterial12, which is more than three 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), BIST and CSIC, orders of magnitude larger than the highest value Campus UAB, Bellaterra, Barcelona 08193, Spain 2 obtained with quantum wells1–3. In order to fully exploit Department of Applied Physics, TU Eindhoven, Den Dolech 2, 5612 AZ Eindhoven, The Netherlands these ultrahigh THz nonlinearities towards practical Full list of author information is available at the end of the article © The Author(s) 2022 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Tielrooij et al. Light: Science & Applications (2022)11:315 Page 2 of 8 technologies, there are fundamental and technological We study the generation of THz third harmonics, where obstacles that need to be overcome. Important funda- generated photons have triple the energy of incident mental questions pertain to the mechanism that leads to photons, using narrowband THz pulses with a funda- THz harmonic generation in massless Dirac-fermion mental frequency f = 0.5 THz, and a duration of around systems, and to whether THz harmonic generation in 10 ps, generated at the superradiant THz facility TELBE TIs originates mainly from the surface states or from bulk located at Helmholtz Zentrum Dresden Rossendorf13 (see carriers without linear dispersion. For technological “Methods” for details). We focus this light onto a sample applications, it is necessary to significantly increase the and measure the transmitted waveform using electro- amount of generated harmonic power. optic sampling (EOS) with 100 fs laser sampling pulses and a 2 mm thick ZnTe nonlinear crystal. This provides Results information on the temporal evolution of the electric field. Here, we address these issues by comparing THz third- Fourier transformation of this time-domain signal then harmonic generation in a 102 nm thick Bi2Se3 TI sample provides the spectral components of the transmitted light, with the response from monolayer graphene serving as including signal at the third harmonic, 3f = 1.5 THz. We benchmark sample (Fig. 1). The TI sample is optimized point out that in these centrosymmetric materials even towards low bulk carrier density (see Methods). In Fig. S1 harmonics are forbidden. we also show results for 50 nm thick Bi2Te3. We chose In Fig. 1c we show the power of the third-harmonic these material systems because they combine the presence signal for increasing incident fundamental power for of massless Dirac fermions with ultrashort cooling times8, regions with and without metal grating. Both graphene which could be beneficial for harmonic generation. On and TI samples show enhancement due to the metal each sample, we fabricated a region with a metallic grating grating, similar to ref. 12. However, the behavior at high (Fig. 1a, b). The reason for adding a metallic grating is that incident power is remarkably different: Whereas the sig- it gives rise to field enhancement inside the gaps, while nal from graphene completely saturates at high pump not sacrificing too much active area. This leads to intensities, the signal from the TI keeps increasing, increased harmonic generation, as shown in ref. 12. The eventually surpassing the signal produced in graphene by metallic grating consists of metal strips with width wmetal, orders of magnitude. We quantify the maximum con- separated by a gap with width wgap. The graphene sample version efficiency for the TI sample in Fig. 1d, and we find has a fill factor of F ¼ wmetal þwgap  90% and the TI sample wmetal that the third-harmonic signal is ~25 times smaller than has F ≈ 85%. the fundamental. Taking into account the reduced EOS a b Bi2Se3 c Graphene 10–2 metamaterial w/ grating Bi2Se3 metamaterial 3f 3f w/o grating 10–3 Third-harmonic power (a.u.) f f 10–4 10–5 10–6 d 0 Graphene 10 EOS sensitivity Amplitude (a.u.) 10–7 Sample transmission 10–8 10–2 10–9 Electric field 0.0 0.5 1.0 1.5 10–5 10–4 10–3 10–2 10–1 Frequency (THz) Incident power (W) Fig. 1 THz harmonic generation in Bi2Se3 vs. graphene. a, b Schematic experimental layout of third-harmonic generation studies in graphene (a) and Bi2Se3 (b) metamaterials. c The THz third-harmonic power is enhanced by the metallic gratings for both graphene (black squares and lines) and Bi2Se3 (blue circles and lines). However, the power saturates strongly for graphene, and therefore eventually Bi2Se3 outcompetes graphene by around two orders of magnitude. The largest pump intensity corresponds to an incident peak field strength of 300 kV cm–1. d Transmitted spectral amplitude (normalized) after the Bi2Se3 metamaterial sample (blue line), with a 3f signal that is ~4% of the fundamental signal at f. Both the sample transmission (red line) and electrooptic sampling (EOS) sensitivity of the 2 mm thick ZnTe detection crystal (black line) are lower at 3f than at f. Taking into account the lower EOS sensitivity, the field-conversion efficiency is 8%, corresponding to 0.5 mW of generated harmonic power Tielrooij et al. Light: Science & Applications (2022)11:315 Page 3 of 8 a Graphene metamaterial Bi2Se3 metamaterial Electric field Electric field Fundamental Fundamental (a.u.) (a.u.) 0 20 40 0 20 40 Time (ps) Time (ps) b c Graphene Bi2Se3 9 kV cm–1 35 kV cm–1 Electric field (a.u.) Electric field (a.u.) 75 kV cm–1 75 kV cm–1 125 kV cm–1 125 kV cm–1 300 kV cm–1 300 kV cm–1 0 20 40 0 20 40 Time (ps) Time (ps) d e Fig. 2 Time-domain evolution of harmonic generation. a Time-domain dynamics (normalized) of a typical fundamental pulse, as measured through electro-optic sampling. b, c Time-domain dynamics (normalized) of generated third-harmonic signal for graphene (b) and Bi2Se3 (c), for four different incident peak field strengths. Both samples are measured in the region with a metal grating for field enhancement. Graphene displays complex dynamics because of the interplay between electronic and phononic heat accumulation and heat-induced (near-)complete THz transparency, which leads to saturation of harmonic generation. Due to the efficient surface-to-bulk dissipation of electronic heat in Bi2Se3, no heat accumulation occurs, and therefore the envelope of the time-domain signal basically follows that of the incident, fundamental waveform. d, e Schematic illustration of heat accumulation in graphene (d) and efficient surface-to-bulk dissipation in the TI (e) sensitivity at 3f, we infer a field conversion efficiency of significantly, the cooling to acoustic phonons becomes ~8%. With an input power of Pf = 75 mW, this means that slower, and optical phonons can reheat the electronic the TI sample is generating around P3f = 0.5 mW of third- system, resulting in slower cooling14. This leads to heat harmonic power. accumulation and to such high electron temperatures that In order to understand why saturation is much stronger THz absorption becomes nearly zero. Such saturable THz in graphene than in the TI, we examine the electric field absorption in graphene was observed in ref. 15 and can be signals in the time domain. Figure 2a shows the temporal ascribed to several phenomena, as reviewed in ref. 16. evolution of the electric field for a typical incident wave- During periods of negligible THz absorption, we expect form at 0.5 THz. For graphene using the weakest incident negligible THz harmonic generation. Only when the field strength of 9 kV cm–1, and filtering the signal at the temperature has dropped some picoseconds later can the fundamental frequency, the envelope of the transmitted system absorb incident THz light and generate THz pulse is similar to that of the incident waveform (Fig. 2b). harmonic signal. The delicate balance between THz- However, the dynamics change dramatically upon electronic heating and cooling through optical and increasing the incident field strength: the signal in the acoustic phonons during the interaction with the incident time domain becomes stretched, and even has some long- multicycle THz pulse thus explains the dips in the term oscillatory behavior (Fig. 2b). We ascribe this to heat envelope of the signal, and the strong decrease in overall accumulation effects: when the electronic system is heated harmonic generation efficiency. When we examine the Tielrooij et al. Light: Science & Applications (2022)11:315 Page 4 of 8 a 1.0 b 40 c 12 Surface Bulk Field enhancement factor M 0.8 Bulk Relaxation time cool (ps) 30 9 Absorption A (%) 0.6 20 6 100 µm 0.4 Surface 10 3 Surface 0.2 Bulk 0.0 0 0 500 1000 1500 2000 0 20 40 60 80 100 0 20 40 60 80 100 Electron temperature Te (K) Fill factor F (%) Fill factor F (%) Fig. 3 Surface-bulk interplay. a Calculation results for the cooling time of surface-state carriers by dissipation of electronic heat to bulk carriers, for two sets of input parameters. The inset is a schematic of the microscopic mechanism where surface-state carriers lose energy by Coulomb interaction with bulk carriers, without surface-to-bulk charge transfer. The obtained timescale of this “Coulomb cooling” is in agreement with a recent study of ultrafast carrier cooling of TI surface states8. b, c RCWA simulation results for total THz absorption (b) and field enhancement (c) for a 3D simulation using the properties of bulk carriers (blue, orange, green lines; corresponding to permittivities of 10, 180, 360, respectively) and a 2D simulation using the properties of the surface-state carriers (red lines). The grating effect is much larger for the surface-state carriers, where significant field enhancement occurs. For our sample with fill factor F = 85% (see microscope image in the inset of panel c, where the brighter stripes are the metal grating), we find fivefold field enhancement envelopes of the signal for Bi2Se3 (Fig. 2c), we do not bulk), while allowing for energy exchange. We perform observe any distortions, which means that, in contrast to the calculations with two sets of input parameters, where graphene, the TI sample does not suffer from any det- the Fermi velocity and bulk effective mass are 106 m s–1 rimental heat accumulation effects (Fig. 2d, e). As a and 0.14me, or 0.4 × 106 m s–1 and 0.21me, where me is the result, the TI exhibits much weaker saturation, leading electron mass. For the bandgap, we take 200 meV, and for to a significantly higher harmonic conversion efficiency the permittivity we take ε = 10. The Fermi energy is in the at higher incident powers and thus a larger harmonic bandgap with a carrier density of 1012 cm–2 for the surface output power. We also note that within the framework states. With these parameters, which are based on of thermodynamic harmonic generation, first discussed refs. 7,21, we find a cooling time for surface-state carriers for massless Dirac carriers in graphene6,17, fast electron around 300 fs for an electron temperature between 500 cooling is beneficial for an efficient harmonic genera- and 2000 K (Fig. 3a). This is in good agreement with a tion. Within this framework, continuous electron recent THz-pump optical-probe experiment8. Impor- heating-cooling dynamics occur during the interaction tantly, surface-to-bulk “Coulomb cooling” is faster than with an incident multi-cycle THz pulse. These dynamics electron-phonon coupling in bulk, which is typically lead to a modulating electron temperature, which cor- found to occur on a timescale of a few picoseconds7,22,23. responds to a modulated THz absorption, due to the We point out that this mechanism will likely be less electron-temperature-dependent THz conductivity16,18. efficient when the Fermi energy is located in the bulk of Faster cooling leads to a larger modulation of the elec- the TI. The picture of increased harmonic conversion tron temperature and therefore a larger harmonic gen- efficiency due to efficient surface-to-bulk electronic heat eration efficiency. dissipation is thus consistent with the framework of thermodynamic harmonic generation, specifically enabled Discussion by the massless Dirac surface states, which exhibit ultra- We ascribe the more efficient heat dissipation in TIs to fast cooling. highly efficient Coulomb scattering between surface-state We provide further evidence that the harmonic signal is charges and bulk charges—a process similar to the generated predominantly from the massless Dirac surface interlayer dissipation occurring in multilayer graphene19. states in TIs by simulating the effect of the metallic We calculate the cooling time that corresponds to this grating on the surface and on the bulk, using rigorous surface-bulk interaction mechanism using a Boltzmann coupled wave analysis (RCWA) simulations24 (see equation approach (see ref. 20 for details). The employed “Methods” for details). In the former case, we simulate the electron-electron collision integral conserves particle TI as a purely 2D conducting sheet with the electronic number in each of the separate systems (surface and properties of the surface states, and in the latter case, we Tielrooij et al. Light: Science & Applications (2022)11:315 Page 5 of 8 use a 3D conducting volume with the electronic proper- even larger third-harmonic output powers will be ties of the bulk. The results in Fig. 3b show that the achieved with or without the grating. Nevertheless, it is absorption in bulk is already quite high, and therefore clear that we have not yet seen the limit of what is pos- cannot be increased much by the metal grating. For fill sible, as the generated harmonic power was still increasing factors above 60%, the overall absorption even decreases, up to the largest incident power that was available. because of the reduced active area. For the surface-state In conclusion, we have demonstrated that topological carriers, the intrinsic absorption is lower, because the insulators generate efficient and high-power THz third- carrier density is lower, and therefore the absorption can harmonic radiation, with a field conversion efficiency of be enhanced substantially by the grating. Figure 3c shows ~8% and a power approaching the mW level. We interpret the field enhancement factor Mfield, which is significantly the results within the framework of harmonic generation higher for the surface-state carriers (Mfield ≈ 5) than for through the thermodynamic nonlinearity mechanism of the bulk carriers (Mfield ≈ 2.6). massless Dirac fermions, as introduced in refs. 6,15. We Since the metallic gratings lead to surface-selective attribute the saturation effects in graphene to heat accu- enhancement, this enables us to distinguish between mulation in both electronic and phononic systems. This surface-generated and bulk-generated harmonic signals heat accumulation does not occur in the TI, because of by comparing the experimentally observed harmonic ultrafast dissipation of electronic heat through “Coulomb enhancement with those calculated for harmonics from cooling” of surface states to the bulk. Therefore, the THG the bulk and from the surface. The power enhancement power generated by the TI is surpassing that from gra- factor for the third harmonic in the perturbative regime is phene by orders of magnitude for high incident powers. P given by M3f ¼ P3f3f; no; grating grating ¼ Mfield 6  ð1  FÞ (see ref. 12). Combining the TI with a metal grating results in a har- The (1 - F) term represents the effect of the reduced active monic power enhancement of two orders of magnitude area due to the metal grating, and the 6th power comes for low to intermediate input powers. This is consistent from the scaling between the generated third-harmonic with harmonic generation by surface states, rather than power and the incident field strength. For intermediate the bulk, as the harmonic enhancement by the grating is incident powers, around 1–10 mW, we estimate that this surface-selective. The effectiveness of the grating is of scaling is closer to the 4th power, as the third-harmonic technological importance, as it enables reaching similar power scales roughly quadratically with incident power output intensities for smaller input fields, and can be (Fig. 1c). Therefore, we expect the harmonic power optimized further. The results obtained here for the enhancement factor M3f to be Mfield 4  ð1  FÞ, which gives prototypical TIs Bi2Se3 and Bi2Te3 are likely representa- ~7 for the bulk and ~100 for the surface states. The tive of the broad class of TI materials with massless Dirac experimentally observed harmonic enhancement factor carriers. Our findings thus pave the way for on-chip, for Bi2Se3 (Fig. 1c) and Bi2Te3 (Fig. S1) at an incident technologically friendly high-power THz photonics. power around 1 mW is M3f = 100–300. This is therefore in much better agreement with harmonic generation Materials and methods from surface states than from bulk states, which is also Sample fabrication, Bi2Se3 supported by the conclusion of ref. 7. Moreover, it is The Bi2Se3 sample was grown by molecular beam epi- clear that with a larger fill factor, even larger third- taxy on a semi-insulating two-inch InP(111)B wafer. As harmonic enhancement, and possibly output intensities, substrate we choose the nearly perfect lattice matched will be possible. InP(111)B with a surface roughness root mean square of We note that some signatures of saturation do occur for about 2 nm to achieve a high-quality film with completely the TI-based samples, and in particular for the TI-grating suppressed twinning, which is otherwise a ubiquitous metamaterial (Fig. 1c). This also explains why the effect of defect in Bi2Se3 thin films25. The sample is grown at a the metallic grating is quite weak for the highest incident substrate temperature of 300 °C by co-deposition of 6 N Bi power: the grating works best in the perturbative regime, and 6 N Se by effusion cells at a growth rate of about 1 nm where the third-harmonic power scales with the sixth per second in Se-rich conditions25. After growth, the power of the field enhancement factor. There are a sample is cooled to 140 °C in Se flux to prevent Se out- number of reasons for the observed non-perturbative diffusion at elevated temperatures. The 102 nm thick behavior. First of all, the electron temperature of two- Bi2Se3 film was characterized by X-ray diffraction, trans- dimensional, massless Dirac fermions intrinsically scales mission electron microscopy, atomic force microscopy, sub-linearly with the incident power12. Furthermore, it is X-ray reflectometry, and magneto-optics, as shown in possible that the bulk of the TI becomes increasingly refs. 25,26. Hall measurements at a temperature of 4 K metallic for increasing incident powers, leading to showed an n-type carrier density of 9 × 1017 cm–3 and a absorption of THz harmonic signal that is generated at carrier mobility of 2000 cm2 V–1 s–1 (assuming a one- the front surface. Importantly, it remains to be seen if carrier model)25. This sample has been exposed to air Tielrooij et al. Light: Science & Applications (2022)11:315 Page 6 of 8 since 2012 and has not shown any indication of degra- sample. Using off-axis parabolic mirrors the THz pump dation, which is promising for technological applications. pulse was focused onto the sample surface and had about 400 µm full-width-at-half-maximum spot size. The Spiricon Sample fabrication, Bi2Te3 PyrocamIII was used to image the THz beam profile on the Untwinned Bi2Te3 films were grown on BaF2 (111) sub- sample. The THz radiation transmitted through the sam- strates (1 × 1 cm2) by molecular beam epitaxy under ultra- ples was refocused onto a 2 mm thick ZnTe crystal and high vacuum (UHV) conditions in a dual-chamber reactor standard THz time-domain spectroscopy was used to (Createc) with a base pressure of 2 × 10−10 Torr27. The measure the temporal profile of the radiation. In order not growth was carried out at 400 °C by co-evaporation of Bi and to saturate the electro-optic sampling crystal with the fun- Te using solid cracker cells with a rate of 0.3 QL min–1, while damental radiation, we either used 1.5 THz bandpass filters monitoring surface reconstruction in-situ by 15-KeV reflec- after the sample (for Figs. 1c, 2), or attenuated the funda- tion high-energy electron diffraction (RHEED). Calibration mental and third-harmonic radiation with a pair of wire was implemented using RHEED oscillations and X-ray grid polarizers, having the same attenuation for funda- reflectivity. Prior to growth, the Bi and Te fluxes were sta- mental and harmonic frequencies (Fig. 1d). For the electro- bilized using a quartz-crystal microbalance. The cell/cracker optical sampling, 100 fs laser pulses were used for probing. temperatures were 670/1100 °C for Bi and 261/600 °C for Te. The timing jitter between the TELBE source and the laser The Te/Bi flux ratio was kept around 2, slightly larger than probe pulses was on the order of a few tens of fs. To vary the stoichiometric ratio of 1.5, to compensate for Te deso- the fundamental beam power we used two wire grid rption during growth. The final film thickness was 50 nm. polarizers: the first one (rotatable) to control the incident After growth and to prevent the Bi2Te3 from oxidizing, the THz power and the second one to keep the THz beam sample was transferred to a second UHV chamber through a always vertically polarized. UHV channel, where a thin capping layer, consisting of 3 nm of Al, was deposited and let to oxidize in air. Finally, the Simulations of grating enhancement sample was taken to an atomic layer deposition reactor The absorption and field enhancement in the TI-grating (Cambridge NanoTech) where a final layer of 10 nm Al2O3 structure were calculated via rigorous coupled wave analysis was grown for subsequent fabrication steps. (RCWA) with a 1D periodic pattern under TM polarization via the periodically patterned multilayer package24,28. The Sample fabrication, gratings geometry consisted of a 1 mm thick InP substrate, followed Gold gratings with areas of 1.5–2.5 mm2 and fill factors by the TI, capped with a 10 nm Al2O3 layer. On top of the F up to 85% were fabricated onto Al2O3-capped TI Al2O3 there is the 50 nm thick gold grating, the only pat- samples by direct write laser lithography (Kloé Dilase terned part, with a fixed grating period of 33 µm, and the fill 250). For Bi2Se3 samples, positive photoresist AZ5214E factor F defined as the fraction of the period filled by (MicroChemicals GmbH) was spin-coated on the TI and the gold. Above and below the structure there are semi- baked at 100 °C for 1 min. The photoresist was then infinite layers of air. The complex permittivities ε ¼ ε0 þ iε00 exposed with optimal dose and developed to reveal the TI were taken as εAu = −1.13 × 105 + 3.48 × 107 i for gold29, surface under the exposed regions. For Bi2Te3 samples, εAl2 03 ¼ 9:48 þ 0:047i for Al2O330, and εInP = 12.5 for InP31. the gap regions were exposed and the pattern was inver- We performed calculations considering the TI layer with a ted by post-exposure baking at 120 °C for 2 min followed finite thickness (3D), as well as an interface only (2D), i.e., by ultraviolet exposure (Kloé KLUB 3) with a power bulk vs. surface, respectively. In the case of the 3D imple- density of 35 mW cm–2 for 10 s. Inductively coupled mentation, the real part of the TI permittivity was taken as plasma-reactive ion etching (Oxford PlasmaPro 100 ε0 ¼ 10; 180; 360 (blue, orange, and green lines in Fig. 3b, Cobra) was used to reduce the presence of residues with a c). We varied this because we found no consensus on this mild oxygen plasma prior to e-beam evaporation (AJA value in the literature. Figure 3 shows that this value does Orion) of Ti/Au (5/50 nm). Last, the lift-off was carried in not play an important role in the final result. The imaginary acetone (>2 h) aided with mild ultrasounds (10 s). part was calculated via ε00TI ¼ σ 3D =ðε0 ωÞ, and σ 3D ¼ eμTI nTI , with the mobility μTI = 2000 cm2 V–1 s–1, carrier density THz measurements nTI = 9 × 1017 cm–3 (from ref. 25), and frequency For THz harmonic measurements we used the super- ω = 2π × 0.5 THz. Here, ε0 is the vacuum permittivity and e radiant accelerator-based source TELBE. TELBE was tuned is the elementary charge. In the 2D implementation, the TI to generate narrowband THz radiation with 0.5 THz central was modeled as a single conducting interface, with a 2D frequency and about 10 ps pulse duration. To block any conductivity of σ2D = (1 kΩ)−1 (from ref. 25). In both cases, radiation outside the pump pulse bandwidth two bandpass the resulting absorption was calculated for normal incidence filters with 0.5 THz central frequency and 20% bandwidth and the enhancement factor was calculated with respect to were inserted into the beam path before the investigated the case of no gold grating. Tielrooij et al. Light: Science & Applications (2022)11:315 Page 7 of 8 Calculation of surface-to-bulk “Coulomb cooling” Finally, =mχ s ðq; ωÞ and =mχ b ðq; ωÞ are the imaginary The temperature dynamics of electrons in the surface parts of the density–density response functions of surface state of a topological insulator is studied by considering their and bulk states, respectively, which are calculated with ðsÞ ðsbÞ kinetic equation, i.e., ∂t fk;λ ¼ J k;λ . The collision integral standard methods of linear-response theory32. These ðsbÞ J k;λ describes the non-retarded and number-conserving describe the absorption and emission spectra of the two Coulomb interaction between surface and bulk states: systems, and therefore their product is proportional to the P PPR1       ðsbÞ J k;λ ¼ 2 hA2π2 L 2 ðssÞ ðsÞ 0 δ εk;λ  ε ðsÞ ð bÞ 0 þ ω δ ε 0  ε ðbÞ  ω δ k 0  ~0  q k 1 dωV 0 q;k ;k F 0 ~ k;λ;kþq;λ kþq;λ k ;η ~ 0 k ;η 0 k k ð1Þ z 0 0 0 k 0 ;k~ ;q λ η;η h       i ðsÞ ðbÞ ðsÞ ðbÞ ðsÞ ðbÞ ðsÞ ðbÞ ´ fk;λ fk 0 ;η 1  fkþq;λ0 1  fk~0 ;η0  1  fk;λ 1  fk 0 ;η fkþq;λ0 fk~0 ;η0 where the factor 2 upfront accounts for the spin degeneracy joint probability of a virtual (non-retarded) photon to be of bulk states, ω and q are the transferred energy and emitted by surface states and absorbed by bulk ones. The 0 momentum parallel to the surface, k 0k and k~k are the com- cooling time τ is finally obtained by recasting eq. (2) in 0 ponents of the three-dimensional momenta k and k~ parallel 0 the form Cs ∂t Ts ¼ ðTs  Tb Þ=τ. In the numerical calcula- to the surface, A is the surface area of the topological insu- tions, we set the surface Fermi velocity υF = 0.5 × 106 m s–1, lator and Lz its extension in the third dimension. In Eq. (1), the surface and bulk electron densities to ns = 1016 m–2 and ðssÞ Fk;λ;kþq;λ0 is the squared matrix element of the surface- nb = 0 (which implies that μb = 0), respectively, Tb = 300 K, electron density operator between incoming and outgoing m = 0.21me, where me is the bare electron mass, scattering states. The matrix element of the screened Cou- Δ = 100 meV, and the undoped-topological-insulator dielec- lomb interaction Vq;k 0 ;~k0 is obtained by (i) solving the elec- tric constant εb = 10. trostatic problem for surface electrons in proximity to bulk ones and (ii) by integrating the result over the incoming and Acknowledgements 0 outgoing bulk states of momenta k 0 and k~ . The surface and Parts of this research were carried out at ELBE at the Helmholtz Zentrum Dresden-Rossendorf e.V., a member of the Helmholtz Association. Simone bulk electrons are described with two distinct Fermi dis- ðsÞ ðbÞ Zanotto is acknowledged for developing the RCWA code dedicated to the tribution functions, fk;λ and fk 0 ;η , characterized by chemical treatment of the two-dimensional conducting interface. We acknowledge potentials µs and µb, and temperatures Ts and Tb, respec- discussions with G. Kh. Kitaeva and K. A. Kuznetsov. K.J.T. acknowledges tively. Here, k and k 0 are the respective particle momenta, funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 804349 (ERC StG CUHL), RYC fellowship λ = ±1 denote the two surface bands forming a Dirac cone, No. RYC-2017-22330, and IAE project PID2019-111673GB-I00. ICN2 was while η = ±1 denote the bulk conduction and valence bands. supported by the Severo Ochoa program from Spanish MINECO Grant No. SEV- ðsÞ The surface bands are approximated as εk;λ ¼ λhυF jk j, 2017-0706. A.P. acknowledges support from the European Commission under the EU Horizon 2020 MSCA-RISE-2019 programme (project 873028 where υF is the surface Fermi velocity, while the bulk ones are HYDROTRONICS) and from the Leverhulme Trust under the grant RPG-2019- ðbÞ εk 0 ;η ¼ η½Δ þ jk 0 j =ð2mÞ. Here, Δ is half the bulk band gap 2 363. S.S., K.B., T.K., and L.W.M. acknowledge support by the SFB1170 (DFG project ID 258499086). C.G. and S.O.V. acknowledge support from the and m the bulk effective mass. For simplicity, we assume bulk European Union’s Horizon 2020 FET-PROACTIVE project TOCHA under grant bands to be spherically symmetric and having identical agreement 824140. T.K., L.W.M., and G.V.A. acknowledge the financial support effective masses. Furthermore, we assume them to have no from the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter (EXC 2147, 390858490). spin structure. We obtain the cooling time by multiplying the ðsÞ kinetic equation by the energy of the surface state, εk;λ , and summing over k and λ. The equation then assumes the form Author details 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), BIST and CSIC, Z Z 1 h i Campus UAB, Bellaterra, Barcelona 08193, Spain. 2Department of Applied 4 d2 q Cs ∂t Ts ¼ Vq2 dωω nðsÞ ðωÞ  nðbÞ ðωÞ =mχ b ðq; ωÞ=mχ s ðq; ωÞ Physics, TU Eindhoven, Den Dolech 2, 5612 AZ Eindhoven, The Netherlands. πh ð2πÞ 2 0 3 School of Physics and Astronomy, University of Manchester, M13 9PL ð2Þ Manchester, UK. 4Physikalisches Institut (EP3), Universität Würzburg, Am Hubland, 97074 Würzburg, Germany. 5Institute of Physics, Polish Academy of where Cs is the heat capacity of surface states while Science, Al. Lotnikow 32/46, PL-02668 Warsaw, Poland. 6Helmholtz-Zentrum  1 Dresden-Rossendorf, Bautzner Landstr. 400, 01328 Dresden, Germany. 7ICREA, ns=b ðωÞ ¼ ehω=kB Ts=b  1 are Bose-Einstein distributions Pg. Lluís Companys 23, 08010 Barcelona, Spain. 8Institute for Topological at the surface and bulk temperatures, respectively. Insulators, Am Hubland, D-97074 Würzburg, Germany Tielrooij et al. Light: Science & Applications (2022)11:315 Page 8 of 8 Author contributions 10. Kovalev, S. et al. Non-perturbative terahertz high-harmonic generation in the S.K., G.V.A., and K.J.T. initiated the project; S.S., K.B., and G.K. designed and three-dimensional Dirac semimetal Cd3As2. Nat. Commun. 11, 2451 (2020). fabricated the Bi2Se3 films, T.K. G.V.A., and L.W.M. coordinated all efforts on 11. Kovalev, S. et al. Electrical tunability of terahertz nonlinearity in grapheme. Sci. the Bi2Se3 films; C.G. and S.O.V. prepared the Bi2Te3 films; D.S.R. and S.V. Adv. 7, eabf9809 (2021). performed the fabrication of metal gratings on the TI samples, supervised 12. Deinert, J. C. et al. Grating-graphene metamaterial as a platform for terahertz by K.J.T.; S.K., I.I., A.P., T.O., M.C., and J.C.D, carried out the THz harmonics nonlinear photonics. ACS Nano 15, 1145–1154 (2021). measurements.; A.P. developed the theoretical model for surface-state 13. Green, B. et al. High-field high-repetition-rate sources for the coherent THz cooling dynamics with input from K.J.T.; A.B. performed the RCWA control of matter. Sci. Rep. 6, 22256 (2016). simulations, supervised by K.J.T.; K.J.T. wrote the manuscript with input 14. Pogna, E. A. A. et al. Hot-carrier cooling in high-quality graphene is intrinsically from D.S.R., S.K., G.V.A., and all other authors. limited by optical phonons. ACS Nano 15, 11285–11295 (2021). 15. Mics, Z. et al. Thermodynamic picture of ultrafast charge transport in gra- Data availability pheme. Nat. Commun. 6, 7655 (2015). All data are available from the corresponding authors upon reasonable 16. Massicotte, M. et al. Hot carriers in graphene – fundamentals and applications. request. Supplementary information accompanies the manuscript on the Light: Nanoscale 13, 8376–8411 (2021). Science & Applications website (http://www.nature.com/lsa). 17. Hafez, H. A. et al. Terahertz nonlinear optics of graphene: from saturable absorption to high-harmonics generation. Adv. Optical Mater. 8, 1900771 Conflict of interest (2020). The authors declare no competing interests. 18. Tomadin, A. et al. The ultrafast dynamics and conductivity of photo- excited graphene at different Fermi energies. Sci. Adv. 4, eaar5313 (2018). 19. Mihnev, M. T. et al. Electronic cooling via interlayer Coulomb coupling in Supplementary information The online version contains supplementary multilayer epitaxial grapheme. Nat. Commun. 6, 8105 (2015). material available at https://doi.org/10.1038/s41377-022-01008-y. 20. Principi, A. & Tielrooij, K. J. Ultrafast electronic heat dissipation through surface- to-bulk Coulomb coupling in quantum materials. Phys. Rev. B 106, 115422, Received: 29 June 2022 Revised: 7 September 2022 Accepted: 8 October https://doi.org/10.48550/arXiv.2206.09119 (2022). 2022 21. Lang, M. R. et al. Revelation of topological surface states in Bi2Se3 thin films by in situ Al passivation. ACS Nano 6, 295–302 (2012). 22. Hsieh, D. et al. Selective probing of photoinduced charge and spin dynamics in the bulk and surface of a topological insulator. Phys. Rev. Lett. 107, 077401 (2011). References 23. Glinka, Y. D. et al. Ultrafast carrier dynamics in thin-films of the topological 1. Sirtori, C. et al. Giant, triply resonant, third-order nonlinear susceptibility χ(3)3ω in insulator Bi2Se3. Appl. Phys. Lett. 103, 151903 (2013). coupled quantum wells. Phys. Rev. Lett. 68, 1010–1013 (1992). 24. Zanotto, S. & Tredicucci, A. Universal lineshapes at the crossover between 2. Markelz, A. G. & Gwinn, E. G. Nonlinear response of quantum-confined elec- weak and strong critical coupling in Fano-resonant coupled oscillators. Sci. trons in nonparabolic subbands. J. Appl. Phys. 80, 2533–2535 (1996). Rep. 6, 24592 (2016). 3. Kono, J. et al. Resonant terahertz optical sideband generation from confined 25. Tarakina, N. V. et al. Suppressing twin formation in Bi2Se3 thin films. Adv. Mater. magnetoexcitons. Phys. Rev. Lett. 79, 1758–1761 (1997). Interfaces 1, 1400134 (2014). 4. Kozina, M. et al. Terahertz-driven phonon upconversion in SrTiO3. Nat. Phys. 15, 26. Orlita, M. et al. Magneto-optics of massive Dirac fermions in bulk Bi2Se3. Phys. 387–392 (2019). Rev. Lett. 114, 186401 (2015). 5. Nakamura, S. et al. Nonreciprocal terahertz second-harmonic generation in 27. Bonell, F. et al. Growth of twin-free and low-doped topological insulators on superconducting NbN under supercurrent injection. Phys. Rev. Lett. 125, BaF2(111). Cryst. Growth Des. 17, 4655–4660 (2017). 097004 (2020). 28. Zanotto, S. PPML - periodically patterned multi layer. https:// 6. Hafez, H. A. et al. Extremely efficient terahertz high-harmonic generation in it.mathworks.com/matlabcentral/fileexchange/55401-ppml-periodically- graphene by hot Dirac fermions. Nature 561, 507–511 (2018). patterned-multi-layer (2022). 7. Giorgianni, F. et al. Strong nonlinear terahertz response induced by 29. Ordal, M. A. et al. Optical properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Dirac surface states in Bi2Se3 topological insulator. Nat. Commun. 7, Pt, Ag, Ti, and W in the infrared and far infrared. Appl. Opt. 22, 1099–1119 11421 (2016). (1983). 8. Kovalev, S. et al. Terahertz signatures of ultrafast Dirac fermion relaxation at the 30. Zhai, M. et al. Characterization of nanoporous Al2O3 films at terahertz fre- surface of topological insulators. npj Quantum Mater. 6, 84 (2021). quencies. Opt. Lett. 45, 4092–4095 (2020). 9. Cheng, B. et al. Efficient terahertz harmonic generation with coherent accel- 31. Palik, E. D. Handbook of Optical Constants of Solids (Academic Press, 1985). eration of electrons in the Dirac semimetal Cd3As2. Phys. Rev. Lett. 124, 117402 32. Giuliani, G. & Vignale, G. Quantum Theory of the Electron Liquid (Cambridge (2020). University Press, 2005).