Reshaping the phonon energy landscape of nanocrystals inside a terahertz plasmonic nanocavity

Phonons (quanta of collective vibrations) are a major source of energy dissipation and drive some of the most relevant properties of materials. In nanotechnology, phonons severely affect light emission and charge transport of nanodevices. While the phonon response is conventionally considered an inherent property of a nanomaterial, here we show that the dipole-active phonon resonance of semiconducting (CdS) nanocrystals can be drastically reshaped inside a terahertz plasmonic nanocavity, via the phonon strong coupling with the cavity vacuum electric field. Such quantum zero-point field can indeed reach extreme values in a plasmonic nanocavity, thanks to a mode volume well below λ3/107. Through Raman measurements, we find that the nanocrystals within a nanocavity exhibit two new “hybridized” phonon peaks, whose spectral separation increases with the number of nanocrystals. Our findings open exciting perspectives for engineering the optical phonon response of functional nanomaterials and for implementing a novel platform for nanoscale quantum optomechanics.


Supplementary Note 1: Fabrication of THz gold nanoantenna arrays
Matrices of nanoantenna chains were fabricated by means of the electron beam lithography (EBL) technique. After substrate cleaning in an ultrasonic bath of acetone, a 160-nm thick poly(methyl methacrylate) (PMMA) layer was spin-coated (1800 rpm) on a 500-µm thick, high resistivity (> 10 kΩ·cm) (100)-oriented silicon chip. Subsequently, annealing was performed at 180°C for 7 min. In order to prevent charging effects during electron beam exposure, a 10-nm thick Al layer was thermally evaporated on the PMMA surface. Electron beam direct-writing of the nanoantenna arrays was carried out using an ultra-high resolution Raith 150-Two e-beam lithography and imaging system. During exposure of the nanoantennas, a beam energy of 20 keV was employed, setting the exposure dose to 550 µC·cm -2 . After the Al conductive layer removal in a 1 M KOH solution, the exposed resist was developed in a solution of methyl isobutyl ketone (MIBK)/isopropanol (IPA) (1:3) for 30 s. Electron beam evaporation in a high vacuum chamber (base pressure 10 -7 mbar) was exploited to produce a 5-nm adhesion layer of titanium and a 50nm Au film, with a 0.3 Å·s -1 deposition rate. Finally, the unexposed resist was removed through a conventional lift-off process in hot acetone. To improve the lift-off efficiency, sonication for 2 minutes at 40 kHz frequency was employed. The substrate was then rinsed out in IPA for 30 s and dried in nitrogen. O 2 plasma ashing at 150 W for 10 minutes was used to remove residual PMMA resist and other organic residues. The fabricated two-dimensional arrays cover an area of 200 × 200 µm 2 and are composed of chains of 200-nm wide dipole antennas coupled end-to-end through a fixed gap G x = 30 nm along their long-axis direction and featuring a spacing G y = 8.5 µm along their short-axis direction.
Supplementary Figure 1 A schematic view of a THz nanoantenna array, together with two representative scanning electron microscope (SEM) images, is presented in Supplementary Figure 1. Arrays with different nanoantenna length L, ranging from 4 to 7 µm, were fabricated on the same silicon substrate in order to tune their plasmonic resonance in the THz frequency band comprising the Fröhlich (FR) optical phonon mode of the nanocrystals (NCs).

Methods
The NCs used in the present study were synthesized using a two-steps protocol, in which a CdS shell of desired thickness was grown on pre-synthesized CdS cores, following the procedure reported in Ref. [1] for CdSe/CdS giant-shell NCs.
Synthesis of CdS cores. CdS cores were obtained by mixing 3 g of TOPO, 0.280 g of ODPA and 0.060 g of CdO in a 25 mL tri-neck flask. The mixture was heated at 120°C under vacuum for 1 hour. Then, under nitrogen, the temperature was increased to 380°C to completely dissolve the CdO. When the CdO was dissolved and the solution became colourless, 1.5 mL of TOP was injected into the system and the temperature was allowed to recover to 380°C before the injection of TOP:S solution (0.020 g S+ 0.500 g TOP). The time of growth was set to 10 minutes to get cores of 6 nm in diameter. After the synthesis, the NCs were washed by precipitation and re-dissolution (3 times) in toluene and methanol. The final sample was dispersed in 3 mL of toluene.

Synthesis of CdS@CdS giant shell NCs.
For the synthesis of core@shell CdS@CdS NCs, cadmium and sulfur precursors were prepared by respectively dissolving 0.640 g of CdO in 5 mL of OA and 0.160 g S powder in 5 mL of TOP. Both precursors were then diluted with ODE to get a 0.5 M Cd oleate solution in ODE (solution 1) and a 0.5 M TOP:S solution in ODE (solution 2).
In a 25 mL tri-neck flask a mixture of 3 mL of ODE and 1 mL of CdS cores solution was first degassed at 80°C under vacuum to remove water. Then, under nitrogen, the solution was heated to 300°C and a mixture of solution 1 and solution 2 (0.65 mL each) was injected into the system drop by drop in 1 hour, targeting a shell thickness of 2.35 nm, corresponding to the growth of 7 CdS monolayers (each layer having a thickness of 0.3358 nm, which is equivalent to half the lattice parameter along the crystallographic c axis).
The final sample was washed by centrifugation and re-dissolution with toluene and 2propanol (3 times) and finally suspended in 3 mL of toluene.

Transmission Electron Microscopy (TEM).
Conventional TEM images were acquired on a JEOL JEM-1011 microscope equipped with a thermionic gun at 100 kV accelerating voltage.
Powder X-ray Diffraction (XRD) Analysis. XRD patterns were obtained using a PANalytical Empyrean X-ray diffractometer equipped with a 1.8 kW Cu Kα ceramic X-ray tube, PIXcel3D 2×2 area detector and operating at 45 kV and 40 mA. The diffraction patterns were collected in air at room temperature using a Parallel-Beam (PB) geometry and symmetric reflection mode. All XRD samples were prepared by drop casting a concentrated solution on a zero diffraction silicon wafer.

Inductively Coupled Plasma Optical Emission Spectroscopy (ICP-OES).
ICP-OES analysis performed on an iCAP 6000 spectrometer (ThermoScientific) was used to estimate the concentrations of CdS cores and CdS@CdS Giant Shell NCs solutions. Prior to measurements, the samples were decomposed in aqua regia (HCl/HNO 3 = 3/1 (v/v)) overnight.
Supplementary Figure 2 shows TEM and XRD results obtained for the CdS cores (panels a and b) and the CdS@CdS Giant Shell NCs (panels c and d). As evident from TEM characterizations, the two samples are quite monodisperse, with an average size of (5.6 ± 0.4) nm and (10.2 ± 0.6) nm for the CdS cores and the CdS@CdS Giant Shell NCs, respectively. XRD patterns confirm the crystallinity of both samples, showing an excellent match with reference card ICSD 98-015-4186, corresponding to the hexagonal wurtzite structure with a = b = 4.1365 Å and c = 6.7160 Å lattice parameters. Moreover, the absence of other peaks related to different phases and to precursors' residues certifies the high purity of the used samples. The concentration of the NCs solutions obtained following the previously described protocols was retrieved via ICP-OES. The CdS cores concentration was 34.24 µM, while the concentration of CdS@CdS Giant Shell NCs in the final sample solution was 10.95 µM.

Deposition of CdS@CdS Giant Shell NCs on THz gold nanoantenna arrays
The deposition of the CdS@CdS Giant Shell NCs on the THz gold nanoantenna arrays was carried out via spin coating under ambient conditions and at room temperature. Specifically, the solution obtained as described in the Methods subsection was diluted to a final NCs concentration of 1 µM. Then, 20 µL of the diluted solution were spin cast onto the patterned substrate at 2000 rpm for 60 s (acceleration of 1000 rpm·s -1 ). The procedure resulted in the deposition of a uniform monolayer of CdS@CdS Giant Shell NCs, as shown in Figure 1d of the main manuscript.
Thicker layers were obtained through multi-step layer-by-layer spin coating deposition; to avoid complete re-dissolution of the already deposited NCs, spin coating steps were interspersed with exposure of the sample to 20 µL of TBAI solution in methanol (10 mg·mL -1 ) for 30 s, followed by three rinse-spin steps with methanol, as described in Ref. [2].

Supplementary Note 3: Estimation of the number of NC layers
The estimation of the average number of NC layers over the sample surface was performed by directly counting in the SEM images the number of NCs within the nanocavities (where the THz electric field mainly concentrates over the sample surface, thus representing the key loci of the electromagnetic interaction) and averaging these values for around twenty-thirty nanogaps.
An example of this procedure, shown for a nanocavity covered with a single layer of NCs, is presented in Supplementary Figure 3 For each nanocavity, several SEM images taken with different tilt angles were considered, in order to confirm the positioning of the NCs within the nanocavity and also evaluate their stacking in the case of multiple layers. From the average number of NCs obtained for the reference case of the monolayer n avg 1 , we then estimated the number of layers N EST using the following formula N EST = n avg / n avg 1 , where n avg is the average number of NCs extracted for the three different sample coverage conditions investigated in our work. The obtained values for N EST are summarized in Supplementary Table 1 below.
To further corroborate these estimates, we also performed Electron Dispersive X-Ray (EDX) measurements, by employing a FEI Helios Nanolab 650 microscope working at an accelerating voltage V = 5 kV, current I = 0.4 nA and acquisition time t = 120 s. SEM images of the nanogaps were collected with 250.000× magnification, in order to include in the images the CdS NCs as well as the edges of the gold nanocavities. The relative atomic fractions of Cd and Au extracted from the EDX measurements taken in such a configuration were then used to estimate the equivalent number of NC layers, by fixing as a reference (N EDX = 1) the case of a monolayer of NCs. As can be seen in Supplementary As for the case of the Raman measurements, the number of layers N associated with each experimental data point refers to the actual number of NCs n contained in the specific nanogap under measurement, estimated by directly counting the NCs in the SEM image of that specific nanocavity ( N = n / n avg 1 ).

Supplementary Note 4: THz characterization
THz transmission spectra were collected at the high-brightness synchrotron infrared source SISSI@Elettra. THz radiation is generated at the bending magnet 9.1 and propagated in ultrahigh-vacuum to the SISSI laboratory, where the radiation enters a Bruker 70v FTIR spectrometer equipped with a Si beam-splitter, and is then coupled to a Hyperion infrared microscope. The microscope is provided with an external detector port allowing to install a He-cooled Si bolometer. Thanks to the high brightness of synchrotron light, this set-up allows performing polarized THz microscopy.
Supplementary Figure 4: Schematic layout of the experimental laboratory at the SISSI-mat beamline: a Bruker 70v FTIR spectrometer is directly coupled to the synchrotron source, and then a Hyperion infrared microscope is used to perform THz transmission measurements.
A grating polyethylene polarizer was positioned at the entrance of the microscope, oriented along the main linear polarization axis of synchrotron radiation, corresponding to the electron's orbit plane. Polarized transmittance measurements were performed by orienting the sample mounted on a goniometer, in such a way that the long axis of the nanoantennas was set parallel to the main polarization axis of the synchrotron beam. A Cassegrain objective (15×, mean incidence angle: 18°) was used to focus THz light onto the sample, while an aperture was employed to limit the focal spot size to a diameter of approximately 260 µm. Data were acquired by co-adding 512 scans with a spectral resolution of 4 cm -1 (spectral point spacing of 2 cm -1 ).

Supplementary Note 5: Raman characterization Experimental details of Raman measurements
Raman investigation was performed in a backscattering configuration using a microspectrometer system (Renishaw inVia) equipped with a 150× LEICA PL APO objective (numerical aperture NA = 0.95) and a thermo-electrically cooled CCD as detector (working temperature -60°C). Spectra were collected exciting the system at λ = 632.8 nm with a He:Ne laser; the wavelength was chosen after optimizing the experimental conditions in terms of scattering efficiency, suppression of fluorescence emission from the CdS NCs (bandgap ~2.4 eV), and beam size at the diffraction limit. The laser power was fixed at 1.7 mW with an integration time of 10 sec. The scans along the nanoantennas were performed using an XYZ motorized stage with a nominal 0.01 μm step. The polarization of the exciting laser source was set by means of a half waveplate. Supplementary Figure 6 shows Raman spectra acquired on a resonant plasmonic nanocavity using different laser illumination powers (100% power level corresponds to 1.7 mW on the sample surface). As can be seen, apart from differences in the signal-to-noise ratio, the spectral positions of the two hybridized peaks do not show any power dependence. For power levels below 85%, the two hybridized peaks cannot be distinguished anymore in the Raman spectra.

Evaluation of the Rabi splitting from Raman spectra
Raman measurements were acquired illuminating a sampling area of ~1.39 µm 2 . Rabi splitting values were calculated from the measured spectra using the following procedure.
A spectrum of interest (collected in a nanocavity region) was first corrected through linear baseline removal, and then by the subtraction of a reference spectrum (acquired on a spot just outside the cavity region). In this way, the contributions deriving from the transverse optical (TO), longitudinal optical (LO) and uncoupled FR phonon modes, as well as from the gold antenna luminescence background, were removed from the spectral data.
On the resulting signals (Supplementary Figure 7), a deconvolution in the region 200-300 cm -1 was operated, by fitting the peaks with two Lorentzian curves, thus extracting the hybridized resonance positions. Supplementary Figure 7: Examples of the fitting procedure described above, for the Raman spectra reported in Figure 3h.

Polarization properties of Raman spectra
As can be seen in Figure 2f of the main text (by comparing the spectra taken on the nanoantennas with the one measured on the silicon substrate), no enhancement of the intensity of Raman scattering was observed for the "non-hybridized" phonon modes during the experiments. This behaviour can be explained considering that terahertz nanoantennas are far out of resonance for the exciting visible radiation (λ = 632.8 nm), thus not contributing to antenna-assisted SERS enhancement.
In order to further corroborate this result, additional Raman measurements were taken to compare spectra where the polarization of the exciting visible laser was set either parallel (red line in Supplementary Figure 8 As shown in Supplementary Figure 8, in both cases no sign of SERS enhancement was retrieved. Furthermore, no significant differences in the hybridized peak positions and amplitudes were found, confirming that the phonon hybridization observed in Raman measurements does not depend on the polarization properties of the exciting visible laser. Besides this, additional Raman analysis was performed on resonant plasmonic nanocavities, by placing an output polarizer (analyser) in the path of the scattered light collected from the sample. As can be seen in Supplementary Figure 9, by increasing the angle formed by the main axis of the output analyser relative to the one of the input polarizer (the latter being set either along the long - Supplementary Figure 9a -or short - Supplementary Figure 9b -axis of the nanoantennas in these measurements) a progressive reduction of the intensity of the Raman response is observed for both the main LO peak and the hybridized FR peaks. All these peaks eventually disappear when the main axis of the output analyser is set perpendicular to the input polarization state. Indeed, the LO phonon resonance of CdS NCs is known to show a negligible degree of depolarization in Raman measurements, since the symmetric A1(LO) phonon mode is dominant [3]. These results thus seem to evidence, at least within our experimental configuration, a similar phonon symmetry between the hybridized FR response and the LO mode of the NCs.

Estimation of the Raman signal enhancement for the hybridized resonance
To extract an estimate of the enhancement related to the plasmon-phonon hybridization in our Raman measurements, we employed the following procedure: We can consider the measured Raman signal where 0.006 A = µm 2 is the area of a nanocavity and 1.39 B = µm 2 is the spot size of the micro-Raman excitation laser on the sample. Using Eq. (1), we can then estimate the Raman signal enhancement of the hybridized Raman peaks K ν ± as: The table below summarizes the enhancement values (spanning the range ∼ 100 -400), as extracted from the data reported in Figure 3h of the main manuscript.
Supplementary  [4]. To extract the dielectric properties of the layer of CdS NCs, we employed the following procedure: we prepared a thick layer (2.5 µm) of NCs over a silicon substrate and characterized its THz transmission (Supplementary Figure 10a, blue circles). We then fitted this measurement (red solid line), considering an effective permittivity eff ε for the NC layer derived from the Maxwell-Garnett approximation for a dielectric mixture [5]: where b ε and CdS ε are the permittivities of the background (air in our case, b 1 ε = ) and filling (CdS NCs) materials, respectively, while f is a dimensionless quantity representing the filling ratio of the mixture. For CdS ε , we considered the common formula used to describe the permittivity of a material in the reststrahlen region [6]: where γ is a damping constant, ε ∞ the high-frequency dielectric permittivity, and LO ω , TO ω the LO and TO phonon frequencies, respectively.
In the fitting procedure of the experimental data of Supplementary Figure 10a, we used the following literature values relative to bulk crystalline CdS: ε ∞ = 5.3, TO ω = 240 cm -1 , LO ω = 302 cm -1 [7]. We then tuned the filling factor to match the experimental resonance position at around 7.9 THz, obtaining 0.75 f = , a value that is consistent with the densely packed arrangement of our CdS layer. Finally, we tuned the damping constant to well reproduce the width of the resonance, which led to a value γ = 15 cm -1 .
The NC self-assembled layers under investigation were then simulated in our numerical analysis by employing an equivalent layer with effective permittivity as in Eq. (2) and appropriate thickness. This allowed us to drastically simplify the meshing requirements in the simulation domain and speed up the calculations. Supplementary Figure 10b shows the absorption of a single layer of NCs, obtained adopting this procedure.
Supplementary Figure 10: a, Relative transmittance spectra (blue circles: experimental measurement, red line: fit with Eq. (2)) of a 2.5-µm-thick CdS NC layer. The relative transmittance is obtained normalizing the sample transmission with the one of a bare silicon substrate. b, Absorption spectrum of a monolayer of CdS NCs, as obtained using Eq. (2) and the parameters extracted from the fit in a.

Calculation domain and illumination conditions
In our simulations, we considered nanoantenna arrays on a silicon substrate with spacing G x = 30 nm along the nanoantenna long axis (x) and G y = 8.5 µm in the perpendicular direction (y). The nanoantenna width was set to 200 nm, the height to 55 nm, while the nanoantenna length was varied in the range 4 -7 µm, to tune the array resonance peak position. All the nanoantenna sharp edges were rounded with a radius of curvature of 20 nm, except the ones perpendicular to the substrate in correspondence of the nanocavities, where a curvature of 40 nm was introduced to better resemble the fabricated structures.
The response of an infinite nanoantenna array was obtained by imposing appropriate periodic boundary conditions in the plane of the array. Two perfectly matched layers (PMLs) were instead set as boundary conditions at the top (air) and bottom (silicon substrate) of the calculation domain, to avoid spurious reflections.
For the input illumination conditions in the simulations, we employed a linearly polarized plane wave. The plane of incidence was tilted by 45° with respect to the nanoantenna long axis, while the angle of incidence in this plane was set to 18° (in accordance with the mean incidence angle of the Cassegrain objective, see Supplementary Note 4). The polarization state of the wave was then set taking into account the input polarization (parallel to the nanoantenna long axis) impinging on the Cassegrain objective. We found that this type of illumination well reproduces the transmission response of the arrays under the Cassegrain objective as well as the experimental position of the lattice mode at around 9.8 THz (Supplementary Figure 11). Supplementary Figure 11a shows the nanoantenna extinction efficiency (the ratio of the extinction cross section to the nanoantenna geometrical area) for arrays with different nanoantenna lengths as a function of frequency. Supplementary Figure 11b shows instead the dispersion of the near field enhancement value (the ratio of the local electric field to the background field) in the centre of the nanocavities. In both figures, the presence of a lattice mode at around 9.8 THz is clearly visible.
The lattice modes of a plasmonic array [8,9] arise from the mutual electromagnetic coupling of the individual plasmonic scatterers, associated with a grazing (diffracted) wave in the plane of the array. When a lattice mode is close to the plasmonic resonance of the array elements, the energy transfer between these two excitations can reshape the overall spectral response of the array, also leading to an increase of its resonance quality factor [10][11][12].
Lattice resonances are known to be found in correspondence of the so-called Wood-Rayleigh anomalies [8,9]. For our specific configuration, the lowest-energy anomaly appears at the frequency that satisfies the following equation: The 2D maps of the extinction as a function of frequency and nanoantenna length are shown in Supplementary Figure 13  Finally, the simulated transmission spectra of the bare nanoantenna arrays, as well as of the arrays covered with 1, 2, and 3 CdS NC layers are plotted in Supplementary Figure 14a Figure 15 shows the electric field distribution in the nanocavity area at resonance, on the planes perpendicular to the main simulation axes. Starting from the distribution of the local electric field in the 3D space, the mode volume of the plasmonic nanocavity can be estimated as follows [13]: where the numerator can be read as the full domain integral of the energy density, and the denominator as the maximum value of the energy density extracted in the same domain.
Considering the case of the array with nanoantenna length of 5.75 µm (corresponding to a resonance peak at around 8 THz, well aligned to the NC phonon resonance) we obtain a mode volume at resonance as small as V mod = 1.43 × 10 6 nm 3 . Compared to the geometrical volume of the nanocavity (V geo ≈ 200 × 30 × 55 nm 3 = 0.33 × 10 6 nm 3 ), the mode volume is thus just ~4.3 times larger. A graphical comparison of these two volumes is given in Supplementary Figure 16, in which the effective squeezing of the radiation in the nanocavity region becomes evident.
Supplementary Figure 16: Size comparison between the geometrical nanocavity volume and the mode volume.
Having an estimate of the mode volume, we can also evaluate the THz vacuum electric field within the cavity [14]: where ν res is the resonance frequency, h the Planck's constant, ε and 0 ε the relative and vacuum permittivities, respectively. Using Eq. (6), we find

Supplementary Note 7: Three-coupled-oscillator model
As shown in the main text, due to the presence of the lattice mode at around 9.8 THz, our overall system in THz transmission measurements is better described by a three-coupled oscillator model [15,16]. The following matrix equation can be used to retrieve the characteristic frequency dispersion of the investigated interaction: where pl ν , ph ν , and lat ν are the frequencies of the noninteracting plasmon, phonon, and lattice modes, respectively, while ν are the eigenfrequencies of the coupled system. pl-ph V ( pl-lat V ) is the coupling constant between the plasmon and the phonon (lattice) modes. Finally, α , β and γ represent the coefficients of the basis functions of the bare plasmon, phonon and lattice modes, respectively. By diagonalizing this matrix, the dispersion of the hybrid eigenfrequencies can be readily evaluated. While this model does not include damping, it is worth underlining that such an idealized picture has proven to be successful in describing experimental systems with finite linewidths, as shown in a variety of studies, for both two- [17,18] and three- [15,16] coupledoscillator interactions. From an operative point of view, a real system with finite resonance linewidths can be considered as strongly coupled whenever the two hybridized resonances are clearly visible in the response spectrum (see, for example, [19,20]). This is evidently the case of our plasmon-phonon interaction, since the double peak formation can be distinctly seen not only in the direct THz extinction response, but also in the Raman measurements, with no THz illumination. Such splitting is indeed an intrinsic characteristic of the coupling between the matter transition and the cavity vacuum-field, and thus it manifests itself in any optical property of the system.
In order to compare the above model to our numerical and experimental results, the uncoupled plasmon resonance frequency needs to be estimated. The resonance wavelength Res λ of nanorod-shaped plasmonic nanoantennas typically obeys the following relation: Res eff 2n L λ ≈ [21,22], where eff n is the effective index of the nanorod surface mode. While this equation well describes the behaviour of an isolated nanoantenna in the THz range, an additional parameter needs to be added in the case of an array of interacting nanoantennas: Res eff 2 , n L S λ = + to also take into account the expected red-shift of the resonance induced by the near-field coupling between the nanostructures [23]. In particular, by fitting the resonance wavelength dependence on the nanoantenna length for the case of bare nanoantenna arrays with the above formula (far from the lattice mode), we then obtain 0 eff 2.375 n = and 10.2 S = µm (see red circles and related fit in Supplementary Figure 17). The uncoupled plasmon resonance frequency in our system can be thus written as: pl eff , 2 c n L S ν = + (8) with c being the speed of light in vacuum.
Supplementary Figure 17: Resonance wavelength as a function of the antenna length for bare arrays.
The above relation (Eq. (8)) has been used to fit the numerical and experimental polariton dispersion curves with the model expressed in Eq. (7). Supplementary Figure 18 shows the result of this procedure on the numerical data regarding a monolayer of NCs.
Supplementary Figure 18: Polariton trace fit of the simulation results for the case of 1 NC layer (compared to Figure 2d in the main text, here the axes are extended to frequency values above 10 THz).
From these fits, we can also extract the quantities 2 α , 2 β and 2 γ , which can be used to estimate the relative contribution of the plasmon, phonon, and lattice modes to the overall hybridized state. As can be seen from Supplementary Figure 19 below, at the FR resonance frequency an essentially complete hybridization of the plasmon-phonon modes occurs (i.e.,