Terahertz nanospectroscopy of plasmon polaritons for the evaluation of doping in quantum devices

Terahertz (THz) waves are a highly sensitive probe of free carrier concentrations in semiconducting materials. However, most experiments operate in the far-field, which precludes the observation of nanoscale features that affect the material response. Here, we demonstrate the use of nanoscale THz plasmon polaritons as an indicator of surface quality in prototypical quantum devices properties. Using THz near-field hyperspectral measurements, we observe polaritonic features in doped silicon near a metal-semiconductor interface. The presence of the THz surface plasmon polariton indicates the existence of a thin film doped layer on the device. Using a multilayer extraction procedure utilising vector calibration, we quantitatively probe the doped surface layer and determine its thickness and complex permittivity. The recovered multilayer characteristics match the dielectric conditions necessary to support the THz surface plasmon polariton. Applying these findings to superconducting resonators, we show that etching of this doped layer leads to an increase of the quality factor as determined by cryogenic measurements. This study demonstrates that THz scattering-type scanning near-field optical microscopy (s-SNOM) is a promising diagnostic tool for characterization of surface dielectric properties of quantum devices.


Introduction
To achieve a quantum advantage over classical computers, the fabrication of high-fidelity qubits with long dephasing time and low coherence loss is required [1][2][3]. While superconducting qubits are one of the most mature quantum computing platforms, the realisation of their full potential is currently limited by material losses, many of which are introduced during fabrication. In practice, the primary decoherence mechanisms are concentrated on lossy dielectrics at device interfaces, particularly the substrate-air and metalair interfaces. These dielectric losses limit the low temperature performance of superconducting devices, stunting the progress of scalable quantum computing [3,4]. Silicon is the industry-standard substrate for superconducting and spin qubits [5,6]. However, to quantify and minimize losses introduced during fabrication, non-destructive diagnostic methods are needed for the nanoscale characterization of quantum relevant surfaces and interfaces.
Surface plasmon polaritons (SPPs) are a type of quasiparticles formed by the collective coupling of incident electromagnetic waves to free electrons and are bound to dielectric conductor interfaces [7,8]. Due to their strong confinement at interfaces, SPPs form the basis of highprecision sensors that can detect minute changes in the local dielectric environment. Therefore, a quantitative approach to evaluate surface dielectric properties via knowledge of surface plasmons is possible. To support SPPs, the real part of the complex permittivity must have the opposite sign across an interface under p-polarized electromagnetic radiation [9]. One strategy to meet this condition is doping, which allows control over the polariton dispersion and tune the supported frequency regimes for plasma frequencies within the THz regimes [10]. Therefore, dopants introduced during the fabrication of siliconbased superconducting quantum devices [2,11] can potentially support of plasmons under THz stimulus. To resolve evanescent light-matter responses, scattering-type scanning near-field optical microscopy (s-SNOM) has emerged as a key nanoscopy tool for long-wavelength radiation, such as THz waves [12][13][14][15][16]. It has been employed to characterize surface-confined electromagnetic responses in microstructures [17], study carriers plasmon and phonon polaritons in van der Waals materials [18][19][20], and reveal THz polaritons in topological materials [21,22] at sub-micron resolution.
To begin this study, we report the observation of THz plasmon polaritons in superconducting quantum devices. A multilayer model is considered to describe the fabricationinduced doping layer. To retrieve both the unknown thickness and complex permittivity of the doped layer without the prior assumptions on the form of a permittivity model, a multilayer extraction procedure for THz nanospectroscopy using vector calibration is demonstrated. The observed polaritonic features from THz near-field measurements are due to the fabrication-induced surface doping and thereby due to THz-induced surface plasmon polaritons. Finally, to evaluate the influence of doping on potential loss mechanism for superconducting qubits, low-temperature characterization measurements are performed on the as-prepared and etched quantum devices. Removing the doped layer increases the quality factor of the etched device, demonstrating the utility of the nanoscale characterization method to inform fabrication efforts. Therefore, THz surface plasmon polaritons serve as an indicator of semiconductor device surface quality. Figure 1(a) shows a three-dimensional rendering of the chip under investigation combined with a schematic representation of the s-SNOM geometry above a coplanar microwave resonator in Figure 1(b). The region of interest for our s-SNOM investigations is the etched domain defining the transmission line as shown in Figure 1(c) and this microscopic structure is repeatedly investigated with multiple tips and different samples in different device orientations. Tapping-mode AFM imaging in Figure 1(d) reveals that the device is dominated by a sharp 85-nm step demarcating the aluminium and silicon regions. Turning to the THz scattering signals, we observe a prominent increase in scattered intensity just within the silicon channel for the second-order harmonic (S 2 ) of the signal (Figure 1(f)). By extracting a sequence of line profiles (45 nm/pixel) across the boundary with averaging over 10 adjacent lines (center in the dashed line in Figure 1(d)), we observe a series of oscillations within S 2 signals, as shown in Figure 1(f).

THz near-field hyperspectral analysis
To better understand and identify the source of the observed spatial-varying oscillatory behaviour in the channel, we perform a high-resolution (20 nm/pixel) THz hyperspectral scan on a second identical device. This hyperspectral scan crosses the aluminium-silicon interface as shown in Figure 2. This dataset includes both the phase and amplitude information of s-SNOM scattering signals and allows analysis of the imaginary part of reflectivity (Im{r exp } = Im{re i }). To reconstruct the complex reflectivity (r exp = re i ) from the scattering signals (S n ), we perform a vector calibration (see Methods). Figure 2(a) shows the spectrally averaged line profile from the calibrated reflectivity, |r exp |, between 0.6 and 1 THz (see Figure S3(a)). Similar to the white-light (broadband scattered signals) line profile in Figure 1(g), a series of oscillations is observed well within the silicon channel. The broad transition in the scattered amplitude near the interface is indicative of  Figure S1). In-plane momenta extracted from fitting profiles (black lines) in (b) are displayed in (d).
the tip-induced edge darkening effect, commonly observed in s-SNOM, due to the tip radius (65 nm, from knife-edge tests in Figure S2). To eliminate the spurious feature contribution (edge darkening effect and edge mode) from the tip-edge convolution [23][24][25], we focus our attention to a region greater than 600 nm away from the interface (highlighted in Figure 2a), in light of the calibrated amplitude and phase contour in Figure S3(b) and (c) to get rid of potential edge-induced phase artefact, to study monochromatic responses on the fringes. Figure 2(b) shows the imaginary part of calibrated reflectivity, Im{r exp }, using the spectrally resolved (monochromatic) frequency response (see its contour in Figure S3(d)). The spectrally resolved scattering in this region shows a pronounced sinusoidal oscillation across the sampled spectral range (0.6-1.35 THz). We observe that the fringe spacing varies with the probing frequency and attribute it to the coherent interference between the tip-launched surface waves and the reflection from the interface (Figure 2c).
To better understand the fringes in the silicon channel, we consider a surface doped layer on a nominally undoped high-resistivity silicon substrate, consistent with earlier measurements of carrier concentration for etched silicon surfaces [11]. Theoretically, the in-plane momentum launched by the s-SNOM probe tip is described by a timeaveraged weight function for the duration of the probe's tapping motion, H(t), and follows a bell-shape distribution described by )⟩ t (see Figure S1) [26]. By focusing an incident THz field on the probe tip, an induced polarization appears on the probe which further polarizes the sample (see such a case in Figure S6(a)). Hence, an iterative polarization process occurs between the probesample pair, which can be regarded as a pair of electric dipoles [27] (see an example in Figure S6(c)). The oscillation motion of the tapping probe modulates the distance between this interacting dipole pair and yields re-radiated light accounting for the s-SNOM scattering response. We expect this scattering quantity to be proportional to the scattered probe-sample dipole moment [28], which may be described by finite dipole model (see Figure S6(b) and (d)). Additionally, this quantity contains a frequency dependent enhancement below a certain cut-off frequency originating from the antenna resonance effect in s-SNOM, depending on the probe geometry and the dipolar resonance between the tip-sample pair [29][30][31][32]. Therefore, the s-SNOM scattering response is convoluted with the material properties and instrumental responses.
The measured s-SNOM near-field response in our system is a weighted average of the in-plane (p-polarized) momentum-dependent reflection coefficient, The coherent-coupling mode in the near-field is characterized by the maxima of Im { r p ( , q )} at real-valued q [26]. In the interrogated multilayer system, the combination of thin-film thickness and surface doping effectively alters the plasma frequency to support the excitation of THz polaritons with the tip-launched in-plane momenta. Considering the nanostructure described in Figure 2(c) with frequency-dependent complex dielectric permittivity, its dispersion relation ranging from 0 to 2.5 THz is calculated and shown in Figure 2(d) (details in Methods). By combining the material dependent and tip dependent responses, we can then extract in-plane momenta q that coincide with the observed fringes in our experiment and is consistent with the calculated dispersion relation. These fits are shown in Figure 2(b) (solid lines) and the extracted momenta are shown in Figure 2(d).

Multilayer extraction with THz nanospectroscopy
As the substrate is high-resistivity silicon with a known permittivity at the probing frequency ( Si = 11.6), the fingerprints of the surface polariton in Figure 2(b) originate from the fabrication-induced surface doping layer reported previously [11]. Hence, a quantitative extraction of the thin-film thickness and permittivity is required to further validate the surface dielectric condition to support THz polaritons. We developed a multilayer extraction procedure which follows our previous work describing calibration of SNOM for direct permittivity measurements [33]. The details of this procedure are explained in detail below.
To quantitatively extract both the thickness and complex permittivity of the doping layer, we extend our previous s-SNOM calibration method to the characterization of multilayers. This process is described schematically in the flow diagram shown in Figure 3.
The steps to follow for the characterization of a multilayer structure using an s-SNOM are summarised in Figure 3  knowledge of the superstrate (air: 1 ) and substrate (high-resistivity Si: 3 ). -(4) obtain feasible solutions for the thin film after multilayer characterization and perform further analysis.
Here, we consider a multilayer structure with unknown thin-film thickness (d) and permittivity ( 2 ) under the quasistatic limit with the in-plane momentum (q). Employing our previous vector-nature calibration method, we are able to calibrate the system responses with three error terms (e D , e R , e S ) and further recover the calibrated reflection coefficient as shown in step 1 and 2 in Figure 3(a). The flow chart in Figure 3(b) describes the procedure to further extract the thin-film thickness and complex permittivity from the calibrated reflection coefficient ( exp ). To solve this inverse problem, first, a thin-film thickness (d) is selected. Then, an initial guess of complex permittivity ( 2 ) is used to compute the multilayer reflection coefficient ( theo ). To find the feasible solutions, the deviation between the measured and the computed reflection , is minimized. An optimization problem is solved subject to the knowledge of a typical material: (1) the material is lossy and (2) not actively emitting energy to the environment, as formulated below: This optimization problem is then solved repeatedly for each thickness in the predefined boundary for the thin-film . In this study, we choose 100 nm as the thin-film upper limit. Finally, the optimal ( 2 , d) pairs are returned as the feasible solutions of the thin film for the multilayer structure. For each monochromatic response, an optimal ( 2 , d) pair exists, including both the real and imaginary parts of the permittivity. Therefore, the final extracted results for a thin film would be two parameter planes as shown in Figure 3(c).
The power of this method is demonstrated by the extraction of the permittivity of the thin film layer for the as-prepared silicon device discussed above from multiple harmonics (S 1 -S 5 ) of the near-field scattering signals. Figure 4(a) shows the extracted THz permittivity (real: solid lines, imaginary: dashed lines) between 0.6 and 1 THz (for results above 1 THz, see Figure S4). For the as-prepared sample, the real part of the thin-film permittivity is negative, indicating it holds an opposite sign with the substrate (highresistivity silicon). This sign change corroborates with the dielectric conditions needed to support surface plasmon polaritons. We can then estimate the fabrication-induced doped layer thickness around 9-13 nm from Figure 4(c) according to the retrieved thickness from multiple orders (S 1 -S 5 ) of harmonic signals. As the effective probing depth decreases with the increasing of harmonic orders, the intersection of the retrieved thickness from several higher-order harmonics further confirms the existence of a doped thin film on the high-resistivity silicon substrate.
To quantitatively describe the doped surface, in Figure 4(c), we compute the equivalent conductivity ( ) from the imaginary part of the permittivity, where 0 is vacuum permittivity [34]. To better estimate the dopant concentration, we compute the uncertainty of the low-frequency THz responses (0.6-0.7 THz), as close to direct current conductivities, in our probing frequencies for extracted thin-film thickness. In addition, the thickness-dependent conductivity curve matches a typical exponential decay profile, where A is the initial amplitude and is the decay constant.
To further validate the retrieved information of the thin film, the dispersion curve is reconstructed from the experimental data of the imaginary part of the p-polarized reflection coefficient at real-valued q. A representative thickness (13 nm) and corresponding complex permittivity (extracted from S 3 ) are used to compute the curve in Figure 4(b). The feasible thin-film thickness range is determined by looking for solutions where the solutions converge and coexist across five harmonics as shown in Figure 4(c). Comparing with the theoretical dispersion curve in Figure 4(d), we notice that the retrieved in-plane momenta (circles) of nearfield fringes from hyperspectral measurements matches well with the retrieved dispersion curves from s-SNOM measurements. For more retrieved information from different harmonic signals and thin-film thickness, see Figure S5. We also validate on a known multilayer nanostructure with SiO 2 on Si (TGQ1, TipsNano Co, Estonia). The extracted values of SiO 2 (see Figure S7) agree with the literature and sample specification [35]. We anticipate an increase of extraction accuracy on multilayer nanostructures across multipleorder harmonic signals in the future with advancements in signal-to-noise ratio on THz s-SNOM. Cases for nano-FTIR and s-SNOM measurements on other spectral ranges are subject to the further study.

Device performances at low temperature
Motivated by the extracted conductive thin-film information, we move further to selectively etch the silicon channel of the quantum device to remove the fabricationinduced surface doping. Low-temperature (∼50 mK) transmission measurements, as shown in Figure 5, are performed to determine quality factors of the fabricated resonators with varying degrees of etching. We find three dips around 5 GHz, 6 GHz and 7 GHz, which represent the response of those resonators. These three dips are then fit to obtain resonance frequencies and quality factors (internal (Q i ), external (Q c ), and loaded (Q l ) quality factors, see Methods).  The measured resonator quality factors are shown in Table. 1. Internal quality factor (Q i ), shown in Figure 5(b), is of interest since it quantifies internal losses caused by coupling to two level fluctuations, which may occur due to impurities in the substrate [2,36]. An increase of internal quality factor is observed after the surface etching in the silicon channel, suggesting an effective removal of substrate dissipation loss dependent on the leakage conduction current [37].

Conclusions
In this article, we observe THz plasmon polaritons in prototypical superconducting quantum devices through THz s-SNOM hyperspectral nanospectroscopy and simulations. Our observations are supported by the existence of a 9-13 nm thick doped layer at the surface of the highresistivity silicon substrate, thereby providing the necessary dielectric change of sign required to support polaritons. To verify the thickness and dielectric properties of the doped layer, we introduce a multilayer extraction procedure and quantitatively retrieve both the thickness and permittivity of the doped thin film layer using multiple s-SNOM harmonics. This multilayer extraction procedure does not require the prior knowledge of the permittivity model or dopants identity and is applicable to other devices structures. The convergence of multiple harmonics as well as the agreement between the hyperspectral experiment measurements and theoretical calculations corroborate the existence of nanoscale THz surface plasmons due to the doped thin film on top of the bulk silicon substrate. By selectively etching the surface of the silicon channel to mitigate the fabrication-induced surface doping, an increased quality factor is observed on the etched when compared to the unetched one. This further confirms the existence of fabrication-induced surface doping in our device. This study demonstrates that THz s-SNOM is an ideal tool for the nondestructive near-surface characterization of film thickness and dopant concentration and can be extended broadly to other nanodevices [38][39][40].

Coplanar waveguide fabrication
High resistivity Si(100) (Topsil, floating zone grown, > 10kΩ cm) was prepared by solvent cleaning, Piranha etching, and then etching in 5% buffer oxide etchant for 20 s to generate a hydrogen-passivated surface. The etched substrate was transferred to a Plassys MEB 550S electron beam evaporator and 100 nm of aluminium (1 nm/s by quartz microbalance) was deposited at room temperature. The sample is spincoated with AZ1512-HS resist and patterned using direct-write lithography (Heidelberg Instruments μPG 101). Once the pattern is defined and developed (AZ726), aluminium is selectively etched from exposed regions using an etchant containing 21% deionized (DI) water, 73% H 3 PO 4 3% acetic acid, and 3% HNO 3 by volume. To selectively etch the silicon channel, reactive ion etching (RIE) is performed. The chip is sent to Oxford PlasmaPro80 to selectively etch the surface of exposed silicon substrate. AZ1512-HS resist is used to protect the underlying aluminium film from bombardment of ions. Flow of SF5 (20 ccm) and CHF3 (35 sccm) are introduced in the chamber with a pressure maintained at 10 mTorr. Ionized fluorine atoms react with Si under RF electric field (82 W). An additional piece of cleaned silicon wafer is used to cover the target area of the chip from its top side. When the etching is complete, the remaining photoresist is removed by submerging for 2 min in 60degree VLSI acetone, followed by a 15 s rinse in VLSI isopropanol. The wafer is then dried with nitrogen gas.

THz s-SNOM measurements
Scattering-type near-field optical microscope (Neaspec GmbH, Germany) with a broadband THz time-domain spectroscopy (0.6-1.6 THz) was employed to interrogate the resonator chip. Commercially available atomic force microscopy probes (25PtIr200B-H, Rocky Mountain Nanotechnology LLC, U.S.A.) are used in this study. A 10-ps THz scattering field is collected from each spatial position (step size = 20 nm) by sweeping the optical scanning delay line. Probe tapping amplitude for THz near-field hyperspectral measurements is ∼290 nm.

Momentum-dependent p-polarized reflection coefficient
The momentum-dependent multilayer reflection coefficient where the normal component of the wave vector is momentum-dependent, k z = iq, under the electrostatic limit. Since the scattering signal measured by the receiver is momentum-averaged sampled through multiple periods of probe oscillation, the near-field reflection coefficient is given by [41][42][43]: where qe −2qz is the tip-coupling weight function, z = H(t) + eff expresses the distance between the effective monopole of the probe to the top layer of a multilayer nanostructure probed by an s-SNOM, cos( t)) describes the tip temporal oscillation during the near-field scattering measurements, is the probe oscillation frequency, a is the probe tapping amplitude, and eff describes the distance between the effective monopole and the tip apex. (See Figure S6

Fringe fitting procedure
To retrieve the polariton momenta from fringes in Figure 2(b), we use a Gaussian-type function as a fitting function: where A, B, q ′ and q ′′ are the fitting parameters. In the fitting procedure, x 0 are kept as the constant for all fringes to determine the starting position for the fitting procedure. To avoid the edge artefact and considering the tip radius estimated from the edge-knife test (∼65 nm, see Extended Data, Figure S2), we choose to fit the fringe data starting ∼200 nm away (spatial position = ∼700 nm in Figure S3(d)) from the aluminium-silicon interface (spatial position = ∼500 nm in Figure S3(d)) till to the far end in the Si channel.

s-SNOM system response calibration
To reconstruct the reflectivity from the scattering amplitude, we perform a calibration using a vector method as described by Guo et al. [33] The calibration considers both, the intrinsic errors introduced via the SNOM itself and extrinsic errors from the environment (such as loss due to humidity or shifts in ambient temperature). It also acts as a converter taking the scattering near field signal and transforming it into the equivalent far-field reflectivity. The s-SNOM system response is calibrated using reference signals from a gold mirror, high-resistivity Si and air. These three standards allow us to extract the spatial-dependent complex THz reflectivity from hyperspectral data (see Extended Data, Figure S3). This calibration also allows direct extraction of complex pairs to describe materials under test (e.g., complex reflectivity or permittivity) without the requirement to fit the approach curve for calibrating tip-sample responses in the case of probe geometries are complex, which may be described by two probe geometry response factors in finite dipole model (see f 0,1 in Figure S6(d)) for mid-infrared cases.

Low-temperature characterization of quantum devices
For the low-temperature characterization of these devices, we packaged and mounted this chip onto the mixing chamber of a dilution refrigerator. Then the mixing chamber was cooled down to around 50 mK and the S 21 parameters were measured. We observed three dips around 5 GHz, 6 GHz and 7 GHz, representing the characteristic resonance frequencies of the three difference resonators. The microwave source was attenuated such that the average photon number circulating in the resonators was in the low-power regime (<1000 circulating photons). The three resonances were then fitted to obtain resonance frequencies and internal (Q i ), external (Q c ) and loaded (Q l ) quality factors using algorithms described in Probst et al.      The incident field E in is focused on the tip apex along with a contribution from the sample's surface reflection % E in due to a wavelength-dependent beam focus (which is larger than the tip radius). The scattering signal is proportional to the overall momentum of the probe and consists of induced momentum p probe and % p probe due to the near-field interaction and surface reflection. % is transverse magnetic (TM or p) polarised Fresnel coefficient and = exp(− Δϕ) accounts for possible phase retardation between the direct radiation to/from the probe and reflection from the sample (which is usually regarded as 1 in s-SNOM analyses). (c) Point dipole model: the probe is approximated as a dielectric sphere with radius and effective momentum due to the near-field interaction. The investigated sample is at a distance ( ) below the probe and mimicked as a mirror image with momentum + and the effective tip-probe distance is . (d) Finite dipole model: the probe is approximated as an elongated spheroid (shank length 2 , radius ) with the initial polarisation , induced by the incident field in and the near-field polarisationinduced by the probesample coupling. It assumes only monopoles , andnear the sample surface participate in the near-field interaction. ,,-denotes the distance between monopole ,,-and probe apex, , + and -+ are the corresponding mirror images and ( ) is the realistic tip-sample distance. (e, f) Multilayer treatments on near-field reflections: multiple reflections from the layered system (permittivity ϵ / and thickness d 0 ) are approximated as an effective monopole eff at the distance eff above the surface. eff is the effective distance between eff and the tip apex.