Direct Visualization of Ultrastrong Coupling between Luttinger-Liquid Plasmons and Phonon Polaritons

Ultrastrong coupling of light and matter creates new opportunities to modify chemical reactions or develop novel nanoscale devices. One-dimensional Luttinger-liquid plasmons in metallic carbon nanotubes are long-lived excitations with extreme electromagnetic field confinement. They are promising candidates to realize strong or even ultrastrong coupling at infrared frequencies. We applied near-field polariton interferometry to examine the interaction between propagating Luttinger-liquid plasmons in individual carbon nanotubes and surface phonon polaritons of silica and hexagonal boron nitride. We extracted the dispersion relation of the hybrid Luttinger-liquid plasmon–phonon polaritons (LPPhPs) and explained the observed phenomena by the coupled harmonic oscillator model. The dispersion shows pronounced mode splitting, and the obtained value for the normalized coupling strength shows we reached the ultrastrong coupling regime with both native silica and hBN phonons. Our findings predict future applications to exploit the extraordinary properties of carbon nanotube plasmons, ranging from nanoscale plasmonic circuits to ultrasensitive molecular sensing.


Near-field amplitude and phase contrast of individual carbon nanotubes
In our study we only used near-field phase maps for the data processing. The near-field amplitude values measured at a specific pixel position depend sensitively on the tip-substrate separation. As the tip-sample separation increases, the amplitude value drops. Small non-resonant residues on the surface cause lower local amplitude values. This can alter the amplitude along the nanotube as well, yielding false, misleading artifacts. Fig. S1 presents the (a) AFM topography, (b) 4th harmonic amplitude map, (c) 4th harmonic phase map of the nanotube of interest after two days of measurements. Small dark spots in the amplitude map are well recognizable in the topography.
On the hexagonal boron nitride (hBN) flake (bottom part), we can also see a big blob of dirt collected by the tip causing a serious artifact in the amplitude map but having only a small effect on the phase measurements.
As mentioned in the main text, the contrast of the nanotube was normalized to the values of the underlying substrate both for the silica and hBN region. In order to make this contrast correction easier during data processing, the fast scan direction (horizontal axis) was set parallel to the material boundary.
The post-processing consisted of the following steps (in the case of amplitude correction a similar process can be used, replacing the subtraction operation by division): 1. Align rows: calculating a representative value (here median was used) of each line and subtracting it from all pixel values in the line. It is a very common scan line artifact correction method in SPM techniques. Because the sample is oriented in a way that the material boundary is horizontal, the median value of a line is representing the specific substrate level.
1 Figure S1: (a) AFM topography (b) Near-field amplitude (c) Near-field phase map of the area of interest. The optical images were taken at 920 cm -1 2. Plane leveling: the normalization then happened with a plane subtraction. It removes any remaining tilt and shifts the base value to 0.
3. The phase profiles were extracted along the nanotube starting from the most intense maximum of the interference pattern as demonstrated in Fig. 1 (c) in the main text. 4. The Fourier spectrum of the line profile was calculated by FFT giving the spatial frequency (f) of the oscillation pattern. The plasmon wavelength can be calculated by λ p = 4π/f . 5. The plasmon wavelength was acquired also by measuring the distances between each adjacent maxima through the whole line profile. The average and the standard deviation of these values were calculated.

Details on sample preparation
Our samples consist of horizontally aligned single-walled carbon nanotubes (SWCNTs) on silicon substrate with hBN flakes. The silicon wafer has a 2 nm native oxide layer that is naturally present on the surface. The samples were prepared in the following way: 1. hBN exfoliation: hBN flakes were mechanically exfoliated from a large piece of hBN crystal 2 (HQgraphene, h-BN, BN2A1) using blue nitto tape (Nitto Denko Co., SPV 224P). The tape with the exfoliated flakes was brought in contact with the silicon subtrate, then peeled off leaving a number of hBN flakes with a variety of size and thickness.
2. CNT growth: SWCNTs were grown by chemical vapor deposition technique on an r-cut quartz substrate using iron catalyst nanoparticles described in detail in [1]   Our instrument (NeaSNOM by NeaSpec GmbH) applies the so-called pseudo-heterodyne (PsHet) detection scheme. The illustration of such a setup can be seen in Fig. S2. This detection method allows to separate the pure near-field signal from the background scattering and additionally makes it possible to measure the amplitude and the phase of the scattered waves originating from the near-field interaction. The setup consists of the AFM part with focusing optics and an asymmetric Michelson interferometer as depicted below. The detailed derivation and description of this method can be found in the following references that created the basis to realize such an instrument: [2,3,4]. 3

Dielectric function of silica and hexagonal boron nitride
In our study we used a Lorentz dielectric function model to describe the optical response of silica given by: In the spectral region of interest we used only one oscillator, with ϵ ∞ = 1.88, ω 1 = 1071.9 cm -1 , γ 1 = 29.9 cm -1 . The parameters were found by fitting the dielectric function values in the SI of [5].
In case of hBN, since it is highly anisotropic, its dielectric function is different for the in-plane and out-of-plane direction: The values were taken from [6].

Fitting nanotube oscillator parameters
In order to apply the harmonic oscillator model, we have to acquire the parameters of each oscillator participating in the coupling. For the silica layer the parameters are straightforward to obtain as described in the main text. For the nanotube, it is more difficult because we do not know the dielectric tensor. However, we know the linear dispersion relation from which we calculate an oscillator's eigenfrequency. For the damping γ cnt we have to rely on the experimental dispersion map. To obtain the parameters, one could fit a vertical line cut from the dispersion map and fit the experimental data with the imaginary part of a Lorentz oscillator y = A 2 /(ω 2 i − ω 2 − iγω). As the experimental data suffer from large uncertainty, we used a vertical quasi-line cut. We take a narrow area as shown by the red square in Fig. S2 (a). By averaging the values in the q direction, we can enhance the signal-to-noise ratio a bit. The resulting quasi-line cut and the corresponding fit can be seen in Fig. S2 (b). The acquired damping value is γ cnt = 150 cm -1 . We note that this method overestimates the damping as the averaging broadens the experimental curve. This effect is, however, not really significant. By repeating the same procedure for a theoretical dispersion map we estimate the broadening and overestimation of γ cnt by only around 10% which does not alter our findings about the regime of coupling significantly.
Since one could argue that the oscillator parameter could not be the bare one because we are not far enough from zero detuning, we also calculated the nanotube plasmon oscillator damping in a different way.
We took a few vertical lines around zero detuning as shown in Fig. S3a by the red rectangular area. The resulting function of polariton amplitude versus frequency can be sen in Fig. S3b. One can clearly observe the two hybrid polariton modes with the splitting consistent with the values in the main text. We fitted this spectrum with the sum of two Lorentzian functions described above.
This yields a very good fit as plotted in Fig. S3 (b) by the red curve. In case of zero detuning, the damping of the hybrid polariton is γ pol = (γ cnt + γ SiO2 )/2. Utilizing this formula we can also  calculate the nanotube plasmon damping because the SiO 2 phonon damping is known. The fit provided γ pol = 90 cm -1 and using γ SiO2 = 30 cm -1 , we get the same value for the nanotube damping γ cnt = 150 cm -1 . The acquired splitting from the position of the fitted Lorentzians is also consistent with the Ω = 300 cm -1 value from the dispersion fitting.

6 Classical coupled oscillator model
We describe the coupling of nanotube plasmons and phonon polaritons via a classical model of coupled oscillators. For sake of simplicity we assume that only the nanotube plasmon is driven.
The equation of motion for each oscillator describing their coupling: x SiO2 + γ SiO2ẋSiO2 + ω 2 SiO2 x SiO2 + 2gẋ cnt = 0 (4) We seek the plane wave solution in the form x = Ae −iωt , the equations become where Γ j = ω 2 j − iω j − ω 2 , j stands for either nanotube (cnt) or silica (SiO 2 ) and K = 2igω. Solving for the steady state solution of the equation system, for x cnt the result becomes where f ∝ e/mE loc . From here the polarization induced in the nanotube is P = e · x. We then visualize the dispersion map as Im(P). As can be seen in Fig. 3. it truly describes the experimental map and the spectra calculated using this theoretical map as described in Section S7.

Spectrum calculation from the dispersion map
When polariton interferometry cannot be applied to measure experimental dispersion, the coupling can be extracted from the phase contrast spectrum of the nanotubes. At a specific illumination frequency, the tip launches every plasmon when its momentum is matching the available momentum range of the tip's near field. The phase contrast measured at a given illumination frequency contains the contributions from all these plasmon excitations. A quantity that is proportional to the measured phase contrast can be replicated by summing all Im(P calc (q, ω)) values calculated at a given frequency by taking into account the weights of the momentum components. We found that a Gaussian weight function W q (ω) applied to each row of the dispersion map replicates the experimental phase spectrum quite well (main text Fig. 4 and Fig. 5b). We set the position of the Gaussian weight function to q = 5.5 · 10 5 cm -1 , as this choice matches the experimental dispersion and replicates the contrast spectrum well. The Gaussian function applied to the theoretical dispersion map is shown in Fig. S5. The theoretical phase contrast spectrum C φ (ω) thus calculated  a change in dielectric environment. In the near-field phase maps one can notice a reflection phase shift (RPS), similarly to previous studies [7]. The first maximum of the interference fringe has lower intensity than the second one (see Fig. S1 (c)). This is a mark of reflection phase shift. Our analysis did not focus on RPS because the geometry here is not ideal to a precise study, however, we measured the value of RPS at two distinct frequencies (920 cm -1 , 1400 cm -1 ). We fit the phase profiles with a damped oscillator form of A · e −2πx/(Qλp cos (4πx/λ p + Φ), where x is the distance from the nanotube end, Q is the quality factor, λ p is the polariton wavelength, Φ is the reflection phase shift. The fit in both cases is plotted in Fig. S7. The acquired values of RPS 920 = 0.08π and RPS 1400 = 0.38π, which correlates with the previously reported values in the range of the accuracy of the reported study. 8 As it is shown in the main text the widely accepted criterion used for strong coupling is written where Ω = 2g (mode splitting) [8,9,10]. The obtained values of g SiO2 = 150 cm -1 leads to C SiO2 = 7.7, far beyond the limit of C > 1.
With the uncertainty of the obtained value of g, one could argue that, the criterion for strong coupling might not be satisfied at some point. Let us examine now the lower limit of strong coupling based on the criterion for C. We can say that the transition between weak and strong coupling is C = 1. From this statement, it follows for the value of g: where γ P h is the damping of the specific phonon mode for silica or hBN. This worst case scenario will lead to g SiO2 = 54.1 cm -1 . With these values we can replot the dispersion curves and the polariton spectrum to compare with the data and to visualize how far the strong coupling limit is from reality. The yellow dashed line in Fig. 10.1 (a) shows the dispersion curve for a system with the above provided coupling strength. As can be clearly seen from the figure, it does not matter if we consider the data low or high quality, these low g curves deviate considerably from the measurement. Furthermore, we also presented the splitting at zero-detuning both in the main text and in the supplementary information (Section S5). Here again, the yellow solid line in Fig.   10.1 (b) presents how the zero-detuning spectrum would look like in the case of g SiO2 = 54.1. It is obvious that g must have higher values, thus the coupling strength is unquestionably higher than it is required for the strong coupling (C > 1). We strongly believe that the presented explanation is scientifically sound, the data are clear, regardless of the noise and they prove our original statement unambiguously.
Regarding ultrastrong coupling (USC), it has stricter requirements and we will discuss this below.
In case of USC, the damping does not come into play at all, as it is defined by the normalized coupling strength η = g/ω g ≥ 0.1 where ω g is the midgap frequency [11]. Thus, the lower limit for USC is η = 0.1 which would translate to g SiO2 = 117.6 cm -1 . The orange lines in Fig. 10.1 (a) and (b) show the dispersion curve and the zero-detuning spectrum for these values of g. Already, the comparison of the experimental dispersion and the new low g curves point out that all the intense peaks of the data are above the curve which proves that the applied g is too low. If it is not convincing enough, the zero-detuning spectrum similarly, probably even more convincingly, shows that the splitting is higher than Ω = 2g = 235.2 cm -1 . The automated fitting provided the red curve with splitting Ω = 305 cm -1 , which we even rounded down to 300 cm -1 . Once again, the evidence is clear regardless of the relative data quality.
Here we comment on the measurements and results on the nanotube-hBN part. We show that the derived contrast spectrum can provide the value of g reliably enough to strongly suggest USC and undoubtedly SC. We repeat the approach presented in the previous section for silica but presenting the contrast spectrum. In case of the hBN phonon mode the USC criterion, η = 0.1, would mean in a spectrum significantly different from the measurements. Modeling only SC results in an even worse match between the model and experiment regardless of the signal-to-noise ratio.