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 heterostructures. In ultrathin 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 in highly anisotropic superconductors 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 λ/3000, 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 heterostructures, 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 however 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 the superconducting temperature. 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 twodimensional (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 (T c ). 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 T c in superconductors 4 . Optical probing of these materials in the sub-gap regime are of interest because below T c they are predicted to support low-loss plasmon waves and confinement of the photon field at deep subwavelength scales 13 . Near-field optics brings in the additional advantage of higher spatial resolution along with sub-surface sensitivity 21 . For high-T c cuprates, the inherent non-linearity of the Josephson coupling between the superconducting CuO 2 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 Fig. 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 inplane 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) 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 27 . The latter geometry exhibits the (quasi-)2D to 3D panoply of gapped or gapless plasmon dispersions ω(q, k z ) as dictated by k z , 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-T c copper oxides the normal state transport is extremely anisotropic: it is metallic in the ab-plane (parallel to the CuO 2 layers) and insulating/incoherent in the perpendicular, c-axis, direction. Below T c , copper oxides behave like a series array of intrinsic Josephson junctions 14 . Tunneling between CuO 2 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~5 0 GHz in Ba 2 Sr 2 CaCu 2 O 8 30 to~1.5 THz in La 2−x Sr x CuO 4 31 and even higher in YBa 2 Cu 4 O 8 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. In optimally doped La 1.85 Sr 0 . 15 CuO 4 (LSCO) the ratio of the in-plane to the c-axis superfluid plasma frequencies is ω pS,c /ω pS,ab ≈ 25 at low temperatures 31,34 .

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 atomiclayer-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 LaSrAlO 4 substrates were polished perpendicular to the crystallographic [001] direction. The films are epitaxially locked to the substrate and pseudomorphic; the CuO 2 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 T s = 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. 1c 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 × 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 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 T c . 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 T c this ratio is fairly temperature independent, staying below unity around a valuẽ 0.95. There is an abrupt rise below T c 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 midinfrared photon energies ω ≈ 1000 cm −1 (30 THz or 125 meV), which are above the superconducting gap in LSCO.

Theoretical modeling 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 T c . 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 T c . Increased scattering due to larger far-field reflection coefficients in the superconducting state cannot account for the observations because of their temperature dependence 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 superconductivity-induced 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 nonmonotonic behavior of the relative LSCO/Au signal with reversal of intensity upon crossing T c as well as the energy scale which is consistent with the Josephson plasma energy ω J inferred from farfield 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 c-axis 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 Eq. (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 deepsubwavelength 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 r P (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[−2qz 0 ] where z 0 is the distance of minimum approach, chosen here to be equal to the tip radius r tip~2 0 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 r P (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 r P (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 Fig. 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 c-axis 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 k 2z in Supplementary Eq. (10), the vertical component of the wave momentum inside the slab, which makes r P (q ≈ 1/d film , ω e < ω P;c ) practically equal to the vacuum-film reflection coefficient r 12 . In other words, the pole in Im[r P (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 Figs. 3, 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. Figure 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 bulkcavity 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  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. Figure 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 Fig. 4, in our measurements we change the temperature and use the constant energy of our THz source. At a given temperature below T c , 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 nearfield 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 Fig. 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 T c , at optimal doping we obtain a critical current density of the c-axis junction J c~2 ×10 5 A/cm 2 and a Josephson penetration depth λ J = c/ω J~3 75 μ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.8T c ) 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 SrTiO 3 , 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 longwavelength 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·T c , 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 near-field environment can also enable characterization of confined modes in topological heterostructures 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.

METHODS
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.

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-to-peak) 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~5 Hz and~1 Hz 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.

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.