Surface Josephson plasma waves in a high-temperature superconductor

Electron density oscillations with acoustic dispersions and sustained at boundaries between different media provide information about surface and interface properties of hetero-structures. In ultra-thin metallic films these plasmonic excitations are heavily damped. Superconductivity is predicted to reduce dissipation allowing detection of these resonances. Emerging low-loss interface Cooper-pair waves have been studied before, however, the observation of surface-confined Josephson plasmons has remained elusive. Here, we report on generation and coupling to these excitations in an ultrathin single-crystal film of high-temperature superconductor La1.85Sr0.15CuO4. The film becomes brighter than Au below the critical temperature when probed with sub-gap THz photons. We show that the enhanced signal in the superconducting state, which can be visualized with a spatial resolution better than {\lambda}/3,000, originates from near-field coupling of light to surface Josephson plasmons. Our results open a path towards non-invasive investigation of enhanced superconductivity in artificial multilayers, buried interface states in topological hetero-structures, and non-linear phenomena in Josephson devices.


INTRODUCTION
Propagating sound-like collective modes in superfluids are converted to higher-energy plasma modes in superconductors 1 . Several decades ago, low-lying excitations with linear dispersions were detected in the microwave region, in aluminum films 2 . These are known as Carlson-Goldman modes. They consist of balanced oscillations of supercurrents and normal carriers 3 and are overdamped except in a narrow region close to Tc. In thin films and superlattices, a second type of superfluid acoustic mode should exist, where the electromagnetic field is confined to the interfaces between superconducting and dielectric regions 4 . Two main obstacles make the observation of these excitations difficult for traditional, far-field, optical techniques [5][6][7] . One is sensitivity, i.e. the capability to distinguish surface modes from bulk contributions. The other is generic to acoustic branches and entails overcoming the momentum mismatch between free-space photons and gapless plasmons. This can be realized in principle by coupling the evanescent waves through a prism, by periodic nano-fabricated arrays, or by bringing a nano-sized light source close to the sample 6,8,9 . Artificial periodic corrugation was indeed used to observe the emergence of microwave surface waves in superconductors [10][11][12] . These studies highlighted the role of reduced dissipation in the superfluid state and established these materials as candidates for temperature-tunable plasmonic structures 10,13 . In this work we use the second, scattering-based, approach, see Fig. 1, which has the benefit of increased spatial resolution due to the nanometer probe size.
There is also a practical interest in materials hosting stacks of conductive metallic sheets ordered with atomic scale precision as they can be used as building blocks for plasmonic devices. Along with doped semiconductor superlattices or modern two-dimensional (2D) materials, superconducting cuprates are examples of such layered systems [14][15][16] . Coupling mechanisms and the requirement of large momentum transfers for mapping dispersion branches gave electron energy loss spectroscopy the leading role in the study of plasmons, including pioneering work which enabled the extension of this concept to surfaces and interfaces 17 . Potential applications are met with substantial challenges because of electrical losses in conventional metallic structures 18 . The same problem is present and amplified in the normal state of superconductors with a high critical temperature (Tc). Topological protection was found to strongly suppress damping of acoustic plasmon excitations propagating on the surface of topological insulators 19,20 . Dissipation is also expected to be significantly reduced below Tc in superconductors 4 . Optical probing of these materials in the sub-gap regime are of interest because below Tc they are predicted to support low-loss plasmon waves and confinement of the photon field at deep sub-wavelength scales 13 . Near-field optics brings in the additional advantage of higher spatial resolution along with sub-surface sensitivity 21 . For high-Tc cuprates, the inherent non-linearity of the Josephson coupling between the superconducting CuO2 planes may enable photonic applications exploiting light-matter interactions in the THz range 22 , with superfluid surface plasmon excitations predicted to play a special role 7,23,24 . The main result of this work is the observation of these modes in strongly anisotropic single crystal superconductors, where their energy and dispersions depend crucially on the Josephson coupling between adjacent layers.
Systems of coupled metallic layers harbor a variety of plasmonic branches and here we describe their salient features (see also Supplementary Note 4 and Supplementary Figure 3). There is a very large energy cost associated with excitations of plasmons in bulk metals. If the same materials are confined to 2D, this energy becomes vanishingly small in the long-wavelength limit; the oscillation frequency scales as ω(q → 0) ∝ q 1/2 , where q is the in-plane momentum 25 . Quasi-2D surface excitations confined to planar interfaces between media with dielectric permittivities of opposite sign are also gapless but they display a sound-like dispersion ω(q → 0) ∝ q, see for example ref. 9 . Interestingly, a change from the expected square root to linear dispersion for 2D electronic states is also generated by bulk-crystal screening effects: examples are the acoustic surface plasmons in metals such as those in Be(0001), ref. 26 . A cross-over between these two regimes can be achieved as well in a configuration where the 2D electron gas is screened by the proximity of a conductive planar electrode/gate 12 , or by stacking 2D conductive and insulating layers in alternation on top of one another, ref. 27 . The latter geometry exhibits the (quasi-)2D to 3D panoply of gapped or gapless plasmon dispersions ω(q, kz) as dictated by kz, the momentum perpendicular to the planes 28 . In addition to these propagating guided modes, superlattices also display their characteristic surface-confined acoustic modes, where the electric and magnetic fields decay exponentially with the distance both inside and outside the layered material 6,29 .
In high-Tc copper oxides the normal state transport is extremely anisotropic: it is metallic in the ab-plane (parallel to the CuO2 layers) and insulating/incoherent in the perpendicular, c-axis, direction. Below Tc, copper oxides behave like a series array of intrinsic Josephson junctions 14 . Tunneling between CuO2 planes allows propagation of bulk plasmons only above ωJ, the screened Josephson plasma energy 6,8,23 . Josephson plasma edge can be probed with far-field techniques and its measured value ranges from ωJ ~50 GHz in Ba2Sr2CaCu2O8 (ref. 30 ) to ~1.5 THz in La2-xSrx-CuO4 (ref. 31 ) and even higher in YBa2Cu4O8 (ref. 32 ). Surface Josephson plasma waves (SJPWs) remain the only excitations with acoustic dispersion in these materials. An anisotropic two-fluid model was shown to capture the dielectric properties of the copper oxides 33 : Here 'S' and 'N' stand for the superfluid and normal components, ωp is the plasma frequency (for the c-axis direction ωJ = ωpS,c/√ε∞,c) and Γ is the scattering rate. The first term in Eq. (1) follows from the London model and the second is a Drude-model description of the uncondensed carriers.

Film synthesis and near-field data
We achieved coupling to SJPWs in a 13 nm thick LSCO film by using a customized cryogenic system based on Atomic Force Microscopy combined with Scanning Near-Field Optical Microscopy 35,36 (AFM-SNOM), see Fig. 1, Methods and Supplementary Note 1. For film synthesis, we have used an advanced Atomic-Layer-by-Layer Molecular Beam Epitaxy (ALL-MBE) system, a technique proven to provide the highest-quality, single-crystal LSCO films with atomically smooth surfaces and interfaces 37 . The surface morphology and crystalline structure of the film were monitored by reflection high-energy electron diffraction, in real time. A source of pure ozone was used to ensure sufficient oxidation under high-vacuum conditions. Single-crystal LaSrAlO4 substrates were polished perpendicular to the crystallographic [001] direction. The films are epitaxially locked to the substrate and pseudomorphic; the CuO2 layers are parallel to the LSCO film surface. During the growth, we kept the ozone partial pressure at p = 2 × 10 -5 Torr and substrate temperature at Ts = 650 ºC. To protect the film surface, we cover the films in situ with 10 nm thick layer of gold, deposited at room temperature. The device shown in the inset of Fig. 1b was patterned using photolithography and Argon-ion milling the LSCO film into a 10 mm long, 20 μm wide strip with 64 lateral leads to large contact pads. The protective gold layer was removed from the active LSCO strip by an appropriate gold etch.
The AFM-SNOM system was optimized for laser light access and large collection throughput, as well as for AFM tip positioning and device visualization 38 , see Fig. 1. All data shown here were acquired with photon energies ω = 26.7 cm -1 (0.8 THz or 3.52 meV), corresponding to a wavelength λ = 375 μm. Gold pads contacting Hall bar devices were used as the optical reference material. AFM-SNOM scans with typical size of 1.3 μm ⨉ 8 μm (more than 2 orders of magnitude smaller than the photon wavelength) were taken in the temperature range 20 K < T < 300 K. These images are shown in Fig. 2a,b. The near-field signal was demodulated up to the 4 th harmonic of the AFM tapping frequency, see Supplementary Note 1. We focus here on the data corresponding to the 3 rd harmonic, S3.
In Fig. 2b we compare the temperature dependence of the near-field signal with that of the device resistance, R(T). The corresponding temperatures are indicated to the right of each scan. LSCO regions are to the left and Au contact areas to the right in the panels. The dashed line between T = 35 K and T = 30 K scans separates the data collected at temperatures collected above and below Tc. The error bars are primarily determined by surface roughness and the imperfections at the edge of the sample due to lithography. Good reproducibility in the data is demonstrated by scans acquired at the same temperatures on cooling and warming during different runs for the same contact and also at different LSCO-Au contacts. The main experimental observation is that, when probed with energies of ω = 26.7 cm -1 (0.8 THz), the near-field intensity ratio between LSCO and Au changes substantially when the HTS sample enters the superconducting state, see Fig. 2c. Above Tc this ratio is fairly temperature independent, staying below unity around a value ~0.95. There is an abrupt rise below Tc where the relative intensity reaches a peak value close to 1.3 at T ≈ 25 K. The signal decreases slightly but remains above unity upon further cooling to T = 20 K. No such changes are observed if the sample is probed with mid-infrared photon energies ω ≈ 1,000 cm -1 (30 THz or 125 meV), which are above the superconducting gap in LSCO.

Theoretical modelling and interpretation
The energy scale and the non-monotonic temperature dependence of the near-field signal across the superconducting transition provide strong constraints to possible interpretations of our data. Phonons do not display any abrupt changes at Tc. Moreover, the energies of optical and acoustic branches are either too high or too low to match the energy and range of momenta probed in our experiment 39 . Long wavelength and low energy magnetic excitations are already overdamped at much smaller Sr concentrations 40 . Local antiferromagnetic fluctuations persist at optimal doping 41 but they are at much higher energies, around 0.2 -0.5 eV, and they also do not display any abrupt changes with temperature upon crossing Tc. Increased scattering due to larger far-field reflection coefficients in the superconducting state cannot account for the observations because of their temperature dependence, refs. 31,42 , see also the discussion in Supplementary Note 5. We can also rule out electronic scenarios associated with both the critical behavior of the dielectric function around the percolation threshold 43,44 . A consideration of the energy scales, near-field coupling mechanisms, temperature dependence across the superconducting transition as well as of the parity and momentum selection rules can explain why scenarios invoking more exotic superconductivityinduced collective modes are very unlikely to account for our results 2,45-47 (a more detailed discussion in connection to these aspects can be found at the end of Supplementary Note 5).
In contrast, we show here that SJPWs explain observations quite naturally. This scenario accounts for all key features of our experiment: the coupling mechanism, the abrupt and non-monotonic behavior of the relative LSCO/Au signal with reversal of intensity upon crossing Tc as well as the energy scale which is consistent with the Josephson plasma energy ωJ inferred from far-field measurements. Coupling to plasmonic excitations is enabled by evanescent modes waves with high in-plane momenta in the proximity of the AFM tip and it is well documented in the literature 48 . The peaks in the near-field signal are correlated to energies of surface modes. The dispersion ω(q) for the propagating electromagnetic mode confined to the planar interface between an isotropic and an optically uniaxial material is implicitly given by, see Supplementary Note 3: where ε1 is the permittivity of the isotropic medium while εab(ω) and εc(ω) are the ab-plane and caxis dielectric constants of the anisotropic material. Note that if εab = εc, we recover the textbook dispersion formula of surface polariton in isotropic media, see ref. 9 and Supplementary Equation (7). AFM-SNOM has thus access to both in-and inter-plane charge dynamics, even if just one surface is accessible. Note also that far-field techniques require either large single crystals with surfaces parallel to the c-axis or, possibly, a grazing-incidence configuration for the case of thin films where the c-axis is perpendicular to the surface. Because of deep-subwavelength resolution, crystals or flakes with μm-length dimensions are sufficient for AFM-SNOM.
The energies and spectral weights associated with Josephson plasma modes in our film can be seen in Fig. 3, which shows the energy-momentum plot of the Fresnel reflection coefficient rP(q,ω). The AFM-SNOM momentum form factor is also shown in panel (a) of this figure by a white dashed line. Within the dipole model its functional form is given by q 2 ·exp[-2qz0] where z0 is the distance of minimum approach, chosen here to be equal to the tip radius rtip ~ 20 nm, see also Supplementary Note 6. The closely spaced branches above unity represent propagating modes inside the film, while the feature at ω/ωJ ~ 1, nearly dispersionless in this momentum range, correspond to SJPWs. The vertical axis in Fig. 3a represents energy in units of ωJ = 56.7 cm-1, which is the calculated Josephson plasma frequency for LSCO at T = 0 K. The values of rP(q,ω) were calculated using measurements of c-axis reflectance in bulk crystals 31 and in-plane THz complex conductivity data in single-crystal LSCO films grown by ALL-MBE 34 , see Supplementary Note 2. Our derivation of rP(q,ω) generalizes the characteristic matrix formalism for treating anisotropic materials in stratified media, see Supplementary Note 4. This approach also helps generalize Eq.
(2) for multilayers because even thin-film devices contain at least two interfaces: vacuum-LSCO and LSCO-substrate, in our case. In this configuration there are two SJPW branches whose energies at high momenta can be read off directly from the zeros of the denominator in Eq. (2), inserting for ε1 the values corresponding to vacuum and the substrate, respectively, see also Supplementary Figure 3. A related approach has been used in ref. 49 for obtaining surface mode dispersions in topological semiconductors. An important point for our interpretation is the realization that, in spite of the very different physical origin, the near-field signal tracks with a good approximation the screened Josephson plasma frequency. The reason is that the asymptotic values of the SJPWs at high wave-vectors are pushed very close to ωJ due the strong anisotropy of the ab-plane and caxis plasma frequencies in LSCO.
Thin films allow for the hybridization of surface modes that are either symmetric or anti-symmetric with respect to the reflection in horizontal symmetry plane. In the presence of a substrate these modes cease to have strict odd/even parity but crosstalk between them is in principle still allowed. In our experiment, i.e. for the range of momenta defined by our form-factor (see white dashed line in Fig. 3a) and energies right below the screened c-axis plasma, we are mostly sensitive to the surface mode confined to the top interface. This is because of the large and negative imaginary part of k2z in Supplementary Equation (10), the vertical component of the wave momentum inside the slab, which makes rP(q ≈ 1/dfilm, ω ⪝ ωP,c) practically equal to the vacuum-film reflection coefficient r12. In other words, the pole in Im[rP(q,ω)] associated with the mode propagating on the bottom interface has a vanishing strength. Further mathematical and visual explanations are given in Supplementary Notes 4 and 5 and Supplementary Figures 3 and 4.
Furthermore, the energy scale and the momentum distribution of our experiment (the white dashed line in Fig. 3a) also allow us to distinguish between bulk-cavity and surface modes. Propagating Josephson plasmons inside the film correspond to the higher energy branches in Fig. 3a while surface modes appear for our momentum range as a horizontal bright line. Fig. 3b shows that momentum integration across the full-width-at-half-maximum (FWHM) of the distribution does not affect the surface mode: the green and black dashed lines practically coincide at energies ω/ωJ ~ 1. However, it flattens out the peak structure of the bulk-cavity modes at higher energies. In conclusion, Fig. 3 reveals not only that modes propagating inside the film are at higher energy than the Josephson plasma edge, but also that momentum averaging completely smears out their spectral structure while leaving the surface modes intact because of their dispersionless nature at high momenta. We therefore find that our experiment is sensitive to SJPWs.
Last but not least, the scenario invoking surface Josephson plasmons is also able to explain the temperature dependence observed in Fig. 2b. The calculated temperature and frequency dependence of the near-field signal within the spherical dipole approximation is shown in Fig. 4, see also Supplementary Note 6. In this approximation the entire tip is replaced by a polarizable sphere whose radius is roughly given by the AFM tip apex curvature. The spherical-dipole approach underestimates the experimental contrast and has clear limitations because of its simplifying assumptions. While a quantitative agreement between the model and experiment is not anticipated, we nevertheless expect that the qualitative features of Fig. 2 are captured. Fig. 4 shows that this is indeed the case, attesting to the robustness of our interpretation in terms of SJPWs. The increase of the relative LSCO/Au signal above unity in the superconducting state as well as the non-monotonic temperature dependence are naturally understood in this scenario. As shown in Supplementary Figure 4, in our measurements we change the temperature and use the constant energy of our THz source. At a given temperature below Tc, the reflectance edge sweeps across our energy window and that is also the point where the near-field signal is enhanced. With further cooling, the superfluid density increases, the reflectance edge moves to higher energy, leading to a decrease in the near-field response, which is what we observe experimentally.

DISCUSSION
As expected, the spherical dipole model is not able to account quantitatively our data. For reasonable fitting parameters it renders values for the LSCO/Au intensity ratio that are too small. The shortcomings can be understood and are essentially related to the simplifying assumptions about effective tip polarizability, near-field distribution, field enhancement factors and the neglect of retardation effects. Indeed, increasing the tip polarizability or taking into account a more realistic shape of the AFM tip would lead to a larger LSCO/Au contrast than that shown in Fig. 4, see also Supplementary Figure 5 and the discussion in Supplementary Note 6. Substantial quantitative discrepancies related to the spherical dipole model were already remarked for experiments performed in the infrared range and more realistic models were proposed for the AFM tip geometry 50 as well as for the electrodynamic response of the sample-probe system 51,52 . Because in our experiment the wavelengths reach millimeter range, the shortcomings of this model are further exacerbated. Therefore, we believe that a meaningful quantitative agreement requires a full electrodynamical treatment taking accurately into account the shapes of AFM tips and the details of the scanning regime, which is outside the scope of this study.
Using the value of the Josephson plasma frequency ωJ = 0.8 THz at T = 0.8 Tc, at optimal doping we obtain a critical current density of the c-axis junction Jc ~ 2·10 5 A/cm 2 and a Josephson penetration depth λJ = c/ωJ ~ 375 μm, indicating that our 20 μm wide wire is in a quasi-1D long Josephson junction regime. For insulating barriers, ωJ tracks, in turn, the superconducting gap Δ.
Using the normal-state conductivity evaluated from far-infrared reflectance data 31 the gap can be estimated 53 to Δ(T = 0.8Tc) ~ 24 cm -1 (0.75 THz), see also Supplementary Note 2. This is a value which is close to the frequency of AFM-SNOM data and it is in very good agreement with recent measurements of the superconducting gap in LSCO-based tunnel junctions 54 . Note that the momentum distribution probed in AFM-SNOM experiments depend on the probe geometry, see Fig.  3a. A systematic exploration of the strong coupling (polaritonic), regime of the SJPW modes is important for probing superconducting nodes or analysis of dissipation channels. This task can be accomplished in several ways. One is by taking advantage of the fact that parameters such as the conical shapes or apex radii can be engineered for etched wire tips 55,56 . Other routes can be provided by 'tilting' the light line dispersion by using thin cover layers of high dielectric permittivity such as SrTiO3, or by using hyperspectral near-field imaging extended to the THz range 57 . It is important to note that Eq. (2), which deals with dispersion of excitations confined to planar interfaces, has to be amended when the wavelength of these modes is comparable to the sample width, which defines the Sommerfeld or Goubau regime 11,58,59 . A comparison with the momentum selective form-factor of the AFM tip, white dashed line in Fig.3, shows that in our configuration when the light couples to the sample only through the tip, the contribution of the surface modes in the long wavelength limit is vanishingly small. Consequently, Eq.(2) is applicable here for a broad momentum range which includes the relevant contributions for analysis. Note also that the extracted values for the Josephson penetration depth and superconducting gap are not obtained from direct measurements, but rather indirectly from the value of ωJ and theoretical considerations in ref. 53 . Measurements of samples with different widths combined with the use of engineered AFM tips with larger apex radii provide a clear path for checking experimentally the role of boundary conditions when exploring the strong light-matter regime.
The energy scales and the temperature dependence of the LSCO/Au near-field ratio constrain the interpretation sufficiently so that we can assign our observations to SJPWs. Other signatures compatible with this interpretation may exist in theory, for example real-space interference patterns of propagating plasmonic excitations 8 . We did not find observe such fringes experimentally. The reasons for their absence can be understood as the result of the d-wave gap and strong quasi-particle damping at temperatures T ≈ 0.8·Tc, which are close to the base temperature of our system. Customized AFM tip radii may be also needed to maximize the coupling efficiency. Plasma waves such as the ones reported here for LSCO have not been observed yet in Fe-based superconductors 60 , although they are layered materials as well. It would therefore be interesting to explore the coupling to these excitations as a function of the sample anisotropy and damping in engineered heterostructures as well as in other cuprate families.
In addition to deep sub-diffractional resolution, light demodulation at different harmonics of the AFM tapping frequency provides dielectric depth information 21 . Thus, cryogenic AFM-SNOM allows probing both the lateral and vertical spatial (in)homogeneity in superconductors, offering a path for mapping superfluid profiles confined to buried interfaces 15 . A variable-temperature nearfield environment can also enable characterization of confined modes in topological hetero-structures and quantification of length-scales associated to superconducting proximity effects in exposed surfaces or in buried interfaces 61 . Non-linear effects involving Josephson plasma in layered superconductors have been observed in LSCO by far-field methods 7,22 and they are predicted to be even stronger below ωJ. It is conceivable that electromagnetic field enhancement around the tip will allow further avenues for exploration and control of such photonic effects in the THz range.

A. La 1.85 Sr 0.15 CuO 4 (LSCO) film synthesis and device fabrication.
Optimally-doped, LSCO films with a thickness of 13 nm were grown on LSAO substrates using ALL-MBE. The stoichiometry was controlled before and during deposition by a scanning quartz crystal monitor and by a calibrated custom-built 16-channel atomic absorption spectroscopy system, respectively. Film thickness was monitored in-situ by counting the reflection high-energy electron diffraction oscillations and post growth by AFM, profilometry and the finite-thickness oscillations in X-ray reflectance and diffraction. The latter also attest to the film being atomically smooth. The substrate temperature during the growth was 650 ºC and the ozone partial pressure was kept at p = 2 × 10 -5 Torr. The device, part of which is shown in Fig. 1b was patterned using photolithography into a narrow strip with 32 pairs of lateral contact pads, as described in the main text. The data shown in Fig. 2 were acquired from an area in the vicinity of one such lateral contact (marked by a red rectangle in Fig. 1b) and the same results were reproduced for other contacts.

B. AFM -SNOM measurements.
Near-field and AFM data were acquired simultaneously in a custom-built AFM-SNOM setup. The sample chamber was optimized for laser light access to the AFM tip, sample visualization and a large numerical aperture in the collection path. AFM tips are reproducibly obtained by etching 100 μm diameter W wires in a 2M-NaOH solution and using a two-ring geometry. Voltages (peak-topeak) in the 2 -10 V range and frequencies between 5 Hz and 800 Hz can be used for shaping the AFM tip shaft and apex radius. The tip radius used here is ~20 nm as imaged by Scanning Electron Microscopy. The tips were glued to a piezo-actuated quartz tuning fork (TF) resonator operated at its resonance frequency of 28.38 kHz (the resonance frequency without the mounted tip was 2 15 Hz = 32.768 kHz). AFM data were taken in amplitude-modulation mode with the feedback based on the TF piezocurrent. The resonance FWHM was ~ 5Hz and ~ 1Hz at room and low temperatures, respectively. Q-control was used to broaden the TF resonance when the system was cold. The AFM oscillation amplitude was A ~ 80 nm as determined from both AFM approach curves as well as from precise interferometric measurements. Light from a Backward Wave Oscillator (Microtech Instruments) with a wavelength λ = 375 μm (800 GHz = 26.7 cm -1 ) was focused on the tip with the help of an off-axis Au-coated copper paraboloid. The SNOM signal was focused on a LN2 -LHe cooled hot electron bolometer and demodulated in real-time up to the 4th harmonic of the AFM tapping frequency with the help of a home-made data acquisition system. We analyzed the SNOM data within the dipole model, which roughly approximates the metallic AFM tip with a sphere whose polarizability is modulated by the interaction with the LSCO sample. Details about the model can be found in the Supplementary Note 6.

C. AFM -SNOM Data Acquisition (DAQ).
A custom made DAQ system was built in order to analyze the SNOM signal in-real time, simultaneously with AFM topography data. The hardware consists of a PC computer, a National Instruments high-speed digitizer and a microcontroller. Custom software provides the graphical user interface and synchronizes the operation of the system. The reference clock of the digitizer, reduced to its sampling frequency of a few MHz, is sent to the microcontroller which serves as a digital counter. The microcontroller also triggers the digitizer to start recording time-series data upon receiving pixel pulses from the ASC500 unit, which is the Attocube scanning probe microscope (SPM) controller. At each pulse, the microcontroller sends the current value of its counter to the PC to precisely identify each pixel within the digitizer's data streams. Two data streams are digitized and sent to the PC simultaneously: a reference sine wave output at the tuning fork's resonance frequency (ΩTF) from an output of the ASC500 controller and the actual detector signal. The reference has a dual purpose: (i) it allows for clock-independent measurements of ΩTF for subsequent Fourier analysis of the detector signal at each pixel and (ii) it allows the SNOM-DAQ to operate the SPM in both amplitude and frequency modulation modes. A lock-in method allows demodulation of the SNOM signal at desired multiples of ΩTF. The calculated values for each pixel are then used to construct 2D dielectric SNOM images of the material surface.

Data availability
The authors declare that all data supporting the findings of this study are available within the manuscript and the supplementary information. They are available from the corresponding author upon reasonable request.

Code availability
All numerical codes in this paper are available upon request to the authors. 33       . This is the frequency at which the data in Fig. 2 were acquired. The grey dashed line marks equal near-field intensity from LSCO and Au, as in Fig. 2b. Our calculations reproduce well the experimental observations: a non-monotonic temperature dependence of the normalized S 3 (T,ω) signal and the reversal of contrast in SNOM signal between LSCO and Au upon entering the superconducting state. This is due to the resonant excitation of SJPWs in the LSCO film.

Supplementary Information for Supplementary Note 1. Additional information about AFM-SNOM data: near-field interaction and spatial resolution of our measurements
Supplementary Figure 1 demonstrates that our optical signal originates in the true near-field interaction between the AFM tip and the sample. The near-field interaction is non-linear as a function of the tip-sample distance 1,2,3-5 . As a result, the reflected signal obtained when the AFM is operated in tapping mode contains harmonics of the AFM resonance. Supplementary Figure 1a shows a typical AFM-SNOM sample approach signal demodulated up to the 4 th harmonic of this frequency. All harmonics display very good suppression of the background component, attesting thus for the true near-field character of our optical measurements. Panel (b) in Supplementary Figure 1 shows the ratio of the near-field signal corresponding to the 3 rd and 2 nd harmonics while the AFM tip is in intermittent contact to the sample. This ratio is around 0.25 for both Au and LSCO and it is only very weakly temperature dependent. This value can be used in principle for refining models for the tip-sample interaction, e.g. to help fix adjustable parameters in various theoretical approaches 4,6 . For the Au reference, we used the data from refs. 14 and 15 . The value εAu = -738752 + i·2077715 was chosen based on the plasma frequency ωp,Au = 8.5 eV and the scattering rate ΓAu ≈ 75 cm -1 , for T ≈ 77 K. At room temperature, ΓAu is a factor of ≈ 5 larger, as indicated by transport measurements. Reasonably large variations in the real and imaginary parts of εAu were found to have no impact on the overall results. The dielectric function for the etched Tungsten tip 16 and LSAO substrate 17 used for the dielectric function at ω ~ 26.7 cm -1 are: εW = -12377 + i·12319 and εLSAO = 19.

Supplementary
These optical parameters also allowed us to estimate the c-axis normal state resistivity. Together with the value of the screened Josephson plasma edge ωJ, which allows determination of the critical current along this axis, one could try to estimate the value of the superconducting gap Δ using the Ambegaokar-Baratoff formula, see the discussion section in the main text and ref. 18 . The set of equations we used are: The critical current density jc was estimated from Supplementary Equation (1a) using Ic = jc·S and C = εrε0S/d, where C is the capacitance between two CuO2 planes and d = 6.62 Å is the distance between two CuO2 planes. The relative permittivity εr was determined from the far-IR reflectance data in ref. 7

Supplementary Note 3. Dispersion of surface waves at the interface between an isotropic and uniaxial material
Here we derive Eq. (2) in the main text, which is a generalization of the textbook relation for surface plasmons, see for example ref. 19 . The geometry of the problem is similar to that of Supplementary Figure 3a, but with only the top interface present. Maxwell's equations for time harmonic fields read: where the free carrier contribution is included in the complex dielectric constant shown in Equation (1) of the main text. The currents and surface charges can be obtained from j = σ·E and σsurf (x) = ε0·[E1z(x) -E2z(x)] where the conductivity σ is a tensor with components given by εα = ε∞,α + iσα/ε0ω where α is the ab-plane or c-axis index. It also has two components corresponding to the normal and superfluid parts, respectively. The equations and boundary conditions allow for separation of the transverse electric and magnetic modes. We confine ourselves to the latter as surface modes can be excited only with p-polarized light, i.e. Ey Bx Bz = 0. We look for propagating solutions of Maxwell equations, with the electric and magnetic fields decaying exponentially away from the interface: where 'j' is an index for the two media. Gauss' laws and the boundary conditions give: Supplementary Equation (4) is a set of four homogeneous equations with four unknowns: E1x, E1z, E2x and E2z. The condition for non-trivial solutions and the bulk dispersion relations for the isotropic and uniaxial medium read: Supplementary Equation (5) is used for deriving the dispersion relations for the surface modes.
The first equality in this equation is equivalent to: which is Equation (2) of the main text. If both materials are isotropic, i.e. εab = εc = ε2 one obtains the standard textbook formula for the dispersion of the surface plasmons: The asymptotic value of the dispersion ω(q→∞) in Supplementary Equation (6), i.e. the electrostatic limit, is given by the condition of vanishing denominator  (6), i.e. the asymptotic energy of the surface mode for high momenta, is given with a very good approximation by ωP,c/√ε∞,c. This is marks the energy of the reflectivity edge, which in the superconducting state is the Josephson plasma energy ωJ. This is why the green dashed lines in the anisotropic case of Supplementary Figure 3c are pushed closer to ωJ compared to the plots in panels (a) and (b) of the same figure. This is also the reason why the near-field signal has a peak which is very close to the reflectivity edge, see Supplementary Figure 4.

Supplementary Note 4. Dispersions of bulk and surface plasma waves in films
The JPWs are non-radiative modes, i.e. for these excitations the electromagnetic fields are exponentially damped outside of the film. They are also true normal modes of the system. Here we provide derivations for the basic formulas in connection to guided and surface JPWs in films, see Fig. 3 in the main text. They appear in various forms in previous publications [20][21][22][23][24][25][26][27][28][29][30][31] . Typically they are obtained by writing down the field solutions along with the corresponding boundary conditions and by looking for non-trivial solutions of the system of homogeneous equations, i.e. a similar approach as the one in Supplementary Equations (3) and (4).
In this work we generalize the characteristic matrix formalism 32 for extraction of the reflection coefficient rP(q,ω) to the case of anisotropic media. This approach is useful because it can be applied easily to larger number of layers, which corresponds to the typical experimental situation.
In this section we also illustrate qualitatively the evolution of the dispersions of the two surfaces modes as a function the dielectric properties of the surrounding media. The characteristic matrix for an anisotropic medium such as the film in Supplementary Figure 3a is given by: The transmitted / reflected fields are at z = 0 and z = d, respectively. Using the boundary conditions along with Fresnel equations, Supplementary Equation (8)  As discussed at the end of Supplementary Note 3, anisotropy in the plasma frequencies pushes the asymptotic energies of both these branches close to ωJ. In Supplementary Figure 3d, a relatively small ratio ωp,ab / ωp,c = 3 was chosen, but the effect is nevertheless substantial. For LSCO this ratio is ~ 25, almost one order of magnitude higher, and this is why the energies probed in our experiment are so close to ωJ, see Fig. 3 in the main text and Supplementary Figure 4. Inter-layer electron tunneling gaps the dispersions of the propagating modes in a superlattice 33 , which otherwise obey ω(q→0, kz ≠ 0) → 0. Furthermore, these guided modes appear only when anisotropy is present, see Supplementary Figure 3d, because propagation inside the film along the z-axis is only possible for frequencies where the ab-plane and c-axis dielectric constants have opposite signs, which is the so-called hyperbolic regime. This can be seen by inspection of the second equality in Supplementary Equation (11).

Supplementary Note 5. Additional information about the interpretation in terms of Surface Josephson Plasma Waves (SJPWs)
The interpretation of our data in terms of SJPWs rests on three pillars: (1) the intimate connection to superconductivity: the changes in the near-field signal are activated by the material entering the superconducting state; (2) the energy scales and (3) the non-monotonic temperature dependence of the relative near-field signal of LSCO compared to Au. The LSCO/Au ratio goes above unity upon crossing Tc and then it decreases again displaying a non-monotonic behavior. We found that these observations provide sufficient constraints for the interpretation. In the main text we mentioned other scenarios that could be, in principle, invoked to explain our results. In this section we look at these aspects in more detail. We start by showing why surface Josephson plasma waves can explain the experimental findings in a straightforward way. We then discuss and rule out the possible role played by far-field reflection coefficients, inhomogeneous behavior close to Tc as well as scenarios involving more exotic superconductivity-induced collective modes In the interpretation section of the main text we also argued that SJPWs provide a natural explanation for the data in Fig. 2 of the main text. Here we explain why this is so. The basis for our discussion is Supplementary Figure 4. Panel 4b of this figure shows that the near field signal at a given temperature is peaked very close to the reflectivity edge. In that sense one can say that SNOM signal tracks approximately the frequency of the plasma edge. However, in spite of similar energies it is crucial to realize that the origin of these two signals is very different. The SNOM signal is peaked at energies corresponding to the poles of the Fresnel reflection coefficient rP(q,ω), i.e. at energies of surface modes. The reflectance edge is given by energy of the screened plasma frequency, ωJ = ωP/√ε∞, which is a very different quantity. The reason these two values are so close to each other is the very anisotropic nature of the sample (ωP,ab/ωP,c ≈ 25 at low temperatures for LSCO) and it is explained mathematically at the end of Supplementary Note 3 and graphically in Supplementary Figure 3d. The idea is that the anisotropy makes the asymptotic value of surface modes energy at high momenta, a pole of rP(q,ω), to be very close to the screened c-axis plasma frequency.
With this discussion we can see how Supplementary Figure 4a explains the temperature dependence of the near-field signal in Fig. 2 of the main text. In our measurements we change the temperature and use the constant frequency of our THz source, the black dashed line in panel (a). At some temperature below Tc the reflectance edge sweeps across our energy window and that is the point where the SNOM signal is enhanced. As the superfluid plasma frequency increases further with cooling, the reflectance edge moves to higher energies and the SNOM signal decreases again, showing thus non-monotonic behavior.
Increased scattering due to larger far-field reflection coefficients in the superconducting state cannot account for the observations. Similarly to a normal metal, the temperature-induced changes of the far-infrared in-plane reflectivity in optimally doped LSCO are very small: by less than 2% on cooling from room temperature to T = 5 K and by less than 1% below Tc, ref. 12 . Such changes are very gradual and also barely detectable by conventional techniques, so they cannot be responsible for the observed effect. More dramatic changes occur along the c-axis due to Josephson tunneling 7 . However, an explanation in terms of increased far-field c-axis reflectance is incompatible with the non-monotonic temperature dependence of the near-field signal below Tc. Lowering the temperature will only increase the brightness of the sample which becomes monotonically more reflective with cooling. Furthermore, because our wavelength (λ = 375 μm) is much larger than our scanning areas it is difficult to understand how far-field reflections from a relatively large diffraction limited spot can produce effects that will not divide out in the LSCO/Au ratio as we move the AFM tip only slightly from one material to another. In fact, the approach curves shown in Supplementary  Figure 1 clearly demonstrate that the sample is probed in the near-field, which also means a distribution of finite momenta for the excitations as shown in Fig.3a of the main text. In conclusion, far-field effects can be ruled out.
We conclude this Supplementary Note by discussing other electronic scenarios compatible, in principle, with near-field coupling. Associating our results with the critical behavior of the dielectric function around the percolation threshold 34 is not a viable option. The fluctuation peaks seen in the 10's of GHz range 35 do not survive in the THz region 11 . Consistent with these expectations, in spite of our deep sub-diffraction capabilities, see Supplementary Figure 2, we did not detect experimentally inhomogeneous behavior in our optimally doped sample. The only other conceivable candidates in the electronic channel that are activated by the superconducting transition are those associated either with the c-axis Josephson tunneling, discussed in the next paragraph, or with the opening of the d-wave gap and the appearance of superconductivity-induced collective modes. Examples are Carlson-Goldman excitations 36 , amplitude oscillations of the order parameter (Higgs modes) or Bardasis-Schrieffer excitonic modes 37 corresponding to sub-dominant pairing. A theoretical review of THz near-field coupling to these excitations can be found in ref. 38 .
Here we briefly comment why these more exotic candidates are unlikely to explain our data. The coupling to Carlson-Goldman modes is expected to be weak because of almost complete charge neutrality of this excitation. Furthermore, disorder is also expected to heavily damp this response 39 .
There is no direct optical coupling to pair-breaking excitations or to the Higgs mode of the order parameter 38 . In addition, their energy scale, set by twice the superconducting gap, is larger than our photon energy. The requirements for near-field coupling to excitonic modes also preclude their observation in our experiment. Parity selection rules would forbid an expected s-to d-transition of the pair wave function. More importantly, Cooper pairs in cuprates are three orders of magnitude smaller than the AFM tip radius used here, so there is a large mismatch in the momentum transfer required to excite these modes. In conclusion, an assessment of near-field coupling mechanisms and energy scales shows that superconductivity-induced collective modes are unlikely to explain our data.

Supplementary Note 6 The spherical dipole model used for AFM-SNOM data analysis
This section provides details about the model used to calculate the AFM-SNOM signal and the momentum dependence of the probe-sample coupling in Fig. 3 of the main text. We analyzed the tip-sample interaction within the point dipole mode, i.e. by approximating the AFM tip with a sphere whose polarizability is modulated by the presence of the sample [1][2][3][4]40 , see also Supplementary Figure 5. The effective tip polarizability μ is given by: where α = a 3 4πε0(εT -1)/(εT + 2) is the bare spherical tip polarizability, εT the dielectric constant of the tip (Tungsten in our case; εT = -2144 +i·1015) and 'a' its radius. The tip-sample distance zt = z0 + A[1 + cos(Ωt)] in the tapping mode oscillates with the amplitude 'A' and has the distance of minimum approach z0 ≈ a. The Fresnel reflection coefficient rP(q,ω) was discussed in the previous section. Within the dipole model the near-field signal demodulated at the m-th harmonic of the tapping frequency is given by: Integration over time of the first term in the Taylor expansion of the denominator in Supplementary Equation (14) for the m-th harmonic generates a momentum dependent function Fm (q, z0, A) given by: The AFM-SNOM momentum form-factor shown with the white dashed line in Fig. 3a of the main text is given by Supplementary Equation (15) for m = 3 and for a distance of minimum approach, z0, taken to be equal to the tip radius a ≈ 20 nm.
The main approximation of the spherical dipole model is the replacement of the actual tip with a spherical particle. The two adjustable parameters of the model are the effective tip polarizability ∝ a 3 , where a is roughly given by the AFM apex radius, and z0, which is the distance of closest approach of the sphere center. The latter parameter should be also of the order of the tip radius.
Considering the typical shapes of AFM tips obtained from etched metallic wires, see for example refs. 41,42 , as well as the large value of the radiation wavelengths in the THz range, it is clear that the spherical dipole model has its limitations. These shortcomings have quantitative impact on the effective tip polarizability, near-field distribution, field enhancement factors and on effects of retardation. We believe they are at the origin of the quantitative discrepancy between the experimental near-field contrast in Fig. 2 and the model calculation in Fig. 4 of the main text.
These quantitative limitations were also observed in the infrared range. It is expected that the THz regime amplifies them because wavelengths are in the millimeter range. A more realistic tip-shape is going to improve of the agreement with the experiment. Indeed, increasing the effective tip polarizability in the spherical model or the elongation of a prolate spheroid in the quasi-static approximation 43 was found to lead to an increased near-field contrast. These considerations help us understand that the our model can be improved by a more adequate probe geometry. However, because the reasons described above, we believe that a meaningful quantitative agreement can only be achieved by a full electrodynamical treatment taking accurately into account AFM tip shapes and details of the scanning regime.