Dynamic construction of refractive index-dependent vibrations using surface plasmon-phonon polaritons

One of the fundamental hurdles in infrared spectroscopy is the failure of molecular identification when their infrared vibrational fingerprints overlap. Refractive index (RI) is another intrinsic property of molecules associated with electronic polarizability, but with limited contribution to molecular identification in mixed environments currently. Here, we investigate the coupling mode of localized surface plasmon and surface phonon polaritons for vibrational de-overlapping. The coupling mode is sensitive to the molecular refractive index, attributed to the RI-induced vibrational variations of surface phonon polaritons (SPhP) within the Reststrahlen band, referred to as RI-dependent SPhP vibrations. The RI-dependent SPhP vibrations are linked to molecular RI features. According to the deep-learning-augmented demonstration of bond-breaking-bond-making dynamic profiling in biological reaction, we substantiate that the RI-dependent SPhP vibrations effectively disentangle overlapping vibrational modes, achieving a 92% identification accuracy even for the strongly overlapping vibrational modes in the reaction. Our findings offer insights into the realm of light-matter interaction and provide a valuable toolkit for biomedicine applications.

The thickness of each layer of the stacked antenna also affects the stacked nanoantenna behavior (Supplementary Figure 3a).When the thicknesses of mental antennas are equal (t1=t3, Supplementary Figure 3b), an increase in the thickness of the spacer enlarges the splitting between the two modes.Because the increase in the spacer thickness reduces the strength of the magnetic response between the two metal antennas.When the spacer thickness t2 is fixed and t1≠t3, (Supplementary Figure 3c), splitting occurs at the configuration with larger values of t1 or smaller values of t3.It implies the sensitivity of the stacked nanoantenna to the direction of incidence of the light.The light incident from the thicker antenna to the thinner antenna is easy to excite the mode splitting.These interesting properties of the stacked nanoantenna can be used in the broadband design, that is, to find the starting point of splitting.
Additionally, there is no significant difference observed in the polarization response between the stacked antennas with stacked nanoantenna and the mono-layer nanorod (Supplementary Figure 4).Both the stacked antenna and the single-layer nanorods exhibit a polarization variation range of 35°, within which the resonance strength of the antenna exceeds half of that observed when the polarization angle is 0 (Supplementary Figure 4c,d).
Regarding the angle of incidence range, the stacked nanoantenna also does not affect the effective range of the antenna's angle of incidence (Supplementary Figure 5).Both the stacked antenna and the mono-layer nanorod have a 50° variation range, where the resonance strength is greater than half that at normal incidence (Supplementary Figure 5c,d).
In terms of novelty, we capitalize on the partial overlap of the two modes in the stacked antenna when they are in proximity to each other, thereby achieving broadband resonance.This method is effective and has not been reported in previously reported work.Current broadband strategies, such as supercells 1 and multi-pixel designs 2 , exploit in-plane space to achieve the goals.
Our method adopts an out-of-plane strategy, which efficiently avoids occupying the in-plane space.Consequently, our approach remains fully compatible with all existing methods while offering the potential to significantly enhance the bandwidth of these existing techniques.Supplementary Figure 1.Optical response of stacked mental antennas in a parallel plate configuration.a, Schematic of the parallel nanoantennas.Two resonant modes become excited: a lower-frequency bonding mode and a higher-frequency antibonding mode.b, Experimental results showing the broadband resonance of the parallel nanoantennas when the two resonant modes are in proximity to each other.c,d,e,f, The distribution of the z-axis electric and magnetic field components along the 2D xz plane.(c,e) is known as the antibonding mode, where two antennas have the same direction of (c) the electric field components and (e) the magnetic response between them is weak.(d,f) is known as the bonding mode, where (d) the direction of the electric field is opposite and (f) a significant magnetic response between the parallel antennas is observed.4. Polarization characteristics of stacked nanoantennas.a, Schematic of the stacked antenna and b, common mono-layer nanorod with a varying polarization angle.t1= t2= t3 = 50 nm.W= 300 nm.L =3 μm.c, Polarization angle versus spectrum wavenumber mapping of the stacked antenna.d, The polarization performance of the mon-layer nanorod.Both the stacked antenna and the mon-layer nanorod have a varying range of 35°, in which the resonance intensity is greater than that when the polarization angle is 0. Supplementary Note 2: The radiating oscillator model

A. Modeling plasmon-phonon coupling
Our SP-PhP platform involves the coupling of surface plasmon polaritons and phonon polaritons.To describe this coupling system, we use a classical model comprising two coupled harmonic oscillators that are subject to external forces f1(t) and f2(t).The model includes a dark oscillator that represents the purely dissipative mode of the phonon.The model can be expressed as follows: The excitation, damping factor, and resonance frequency of the bright-mode resonators 1 and 2 are represented by p1,2, γ1,2, and ω1,2, respectively.The two resonators are linearly coupled with each other through a coupling strength μ•exp(iφ), where the phase shift φ between the two resonators is induced by the retardation effect.The dissipative phonon oscillator, which has a resonance frequency of ωm and damping factor γm, is excited by q(t), and its coupling strength with the plasmonic oscillator is represented by κ1,2.We make the assumption that f1(t)= f2(t)= f(t), exp(iφ1)= c1 and exp(iφ2)= c2.In this work, the phonon is coupled with the bonding mode of nanoantennas, so κ1 is equal to 0. Under steady-state conditions, p1(t), p2(t), q(t) and f(t) can be expressed as: where the amplitudes of the oscillator or radiation in the frequency domain are denoted by P1, P2, Q, and f0.By substituting equation (S4) into equations (S1), (S2) and (S3), we arrive at the following: Solving equation (S5) yields the resonant amplitudes of the modes in the coupled system as follows: where Since the dark oscillator lacks a dipole moment that matches the external field, it does not contribute to the surface current.Therefore, we can calculate the electric current density J as follows: 12 () According to equations (S18) and (S19), the surface conductivity can be determined as, We can calculate the scattering parameters, including reflection, transmission, and absorption, using equations (S9), (S10), and (S11).To illustrate, the reflection can be expressed as follows To explore phonon-plasmon coupling, it is essential to calculate the signal strength of the phonon in the coupled system.This can be achieved by subtracting the background spectrum of the plasmon from the spectrum that contains the phonon signal, as follows: where ω-ωm is defined as the detuning δ of the plasmon-phonon coupling system.
Sensitivity refers to the deviation from the ideal slope of a characteristic curve, which measures how much an output changes in response to a change in the input quantity it measures.
In the context of a phonon-plasmon coupling system, the output is the intensity of the reflection signal of the plasmonic resonance, while the input is the target molecules that cause a detuning between the plasmon and the phonon and thereby change the reflection intensity.Therefore, in our case, sensitivity can be calculated as follows: According to equation (S12), we can obtain Figure 2i in the main text.

B. Modeling the phonon polariton branches
The phonon polariton branches can be studied by analyzing two oscillators with eigenfrequencies ω+ and ω-, which are affected by the coupling strength between them.To observe how the vibration parameters are influenced by the coupling strength, we analyze the system of the two oscillators without any driving force.The resulting matrix equations describe the motion of the oscillators and can be expressed as follows: where we can assume |ω-ω2,m|<< ω when investigating the response of the system at or near resonance.To obtain a solution, we can solve the characteristic equation |H-ωI| = 0, where I represents the identity matrix.The dressed frequencies of the system can be expressed as follows: where g= ω+-ωis the frequency splitting.For a lossless (γ2=γm=0) system with m0   and g=κ ω0, we can plot the bare and dressed frequencies of the coupling system, as shown in Supplementary Figure 6.As the coupling strength increases, the split between the upper and lower branches becomes more pronounced.In the strong coupling regime, a clear Rabi split can be observed.According to equation (S14), we can obtain Figure 1c

Supplementary Note 3: More discussion on polariton-molecule interaction
Polaritons are hybrid quasiparticles formed due to the strong coupling of light to matter, including photon polaritons, exciton polaritons, and phonon polaritons 3 .In the IR region, the phonon polaritons and photon polaritons have received a great deal of attention, especially for sensing applications based on the interaction of polaritons and matter 4 .The large dipole moments and long lifetime of polaritons make them strongly interact with analytes in the surrounding environment.Therefore, polariton-based sensing sensors have the advantages of high sensitivity, low detection limit, nondestructiveness, and real-time detection.
The essence of the interaction between polaron and molecule is the interaction between polaron and refractive index and infrared vibration mode.In the case of the interaction with the refractive index, polaritons are sensitive to changes in the refractive index of the surrounding environment, which can be caused by the presence of analytes.When the refractive index of the environment changes, it affects the propagation of the electromagnetic field, and therefore the polariton properties, such as frequency shift and angle change.By measuring the changes in the polariton properties, we can detect the presence and concentration of analytes.In terms of the interaction with the vibrational modes of the molecules, the strong electric field component of the polariton.This interaction leads to changes in the vibrational properties of the molecules, such as changes in their intensity of vibration absorption.In the IR region, we know it as the Surface-Enhanced Infrared Absorption (SEIRA) effect.
Previous demonstrations have aimed to improve performance and push the technological frontier forward.In the case of polaronic refractive index-based sensors, the figure of merit (FOM) can be enhanced by reducing the bandwidth and increasing sensitivity.Methods to reduce bandwidth include utilizing bound states in the continuum (BIC) 5 and higher-order modes, among others 6 .Lowering the refractive index of the substrate is a way to increase sensitivity 7 .For SEIRA-based sensors, performance metrics such as LOD, sensitivity, and bandwidth can be improved by: 1) optimizing the design of the nanoantenna structure.For instance, reducing the distance between metamaterial pattern units can increase the near-field intensity, which is attributed to attractive electromagnetic force interactions 8 ; 2) constructing a perfect absorber to increase near-field confinement and enhancement.In a dipole absorber, the electric and magnetic fields generate electric and magnetic dipoles within and between metal pattern layers, respectively 9 ; 3) improving the spatial overlap between analytes and enhanced near-field by integrating microfluidic techniques with nanoantennas 10 ; 4) selecting appropriate plasmonic and dielectric materials.For instance, using Al as the plasmonic material can form a self-passivating oxide layer on the surface 11 ; and 5) employing machine learning algorithms to simplify complex design and increase the types of detectable targets 12 .
In our work, we develop a distinctive vibration, termed RI-dependent vibration, through the excitation of coupled surface plasmon-phonon polaritons.This unique vibration is a specific type of SPhP oscillation occurring within the reststrahlen band of polar dielectric crystals, exhibiting a notable sensitivity to the molecular refractive index.Notably, sensitivity to refractive index does not imply the simultaneous measurement of the real and imaginary parts of the composite refractive index in the same wavelength range.Our proposed method eliminates the need for any conversion using Kramers-Kronig relations and requires only 1 pixel to encompass refractive index information in our work.The physical mechanism behind the refractive index-dependent vibration response is the detuning of plasmonic resonance and SPhP mode in the reststrahlen band, as shown in Supplementary Figure 7.In our proposed stacked nanoantenna, the SPhP mode (Supplementary Figure 7a) and the plasmonic resonance mode (LSPP, Supplementary Figure 7b) The plasmonic resonance frequency is redshifted, while the SPhP mode frequency is fixed.This will cause a change in the detuning of plasmonic resonance and SPhP mode.The Asymmetric Least Squares algorithm is a potent tool for baseline correction.It is often used on a single spectrum, where the true baseline cannot be determined.In the case of multiple spectra, shared characteristics among them can be inferred to learn a gradually changing baseline, resulting in an accurate baseline.The particular technique involves estimating the baseline by penalizing the difference of the baseline correction signal, based on the resemblance between multiple spectra 13 .The optimization equation is based on the Whittaker Smoother and can be expressed as where y1, y2, ..., ym are spectral data with a length of m, and z1, z2, ..., zm are target baselines.The first term on the right represents the difference between the corrected spectra pair-wise, and optimizing the baseline entails ensuring that each corrected spectrum is in close proximity to the average corrected spectrum.The second term is an asymmetric fitting error that measures the data.
The final term imposes a smoothness constraint on the baseline.By introducing a relaxation factor au to each corrected spectrum, the corrected spectrum can deviate from the average corrected spectrum.Then, equation (S15) is extended to The solution to this equation is where E is the n×n identity matrix, and λ'=m λ, μk'=m μk.
After the initial value of each parameter is determined, the baseline can be fitted for multiple spectra according to equation (S17), and then the weight matrix is estimated as Then, we calculate the baseline using equation (S17) in Matlab software.Supplementary Figure 8a,b   the training data set, and 20% is for the testing data set.In this work, the DNN models for regression are developed using Keras, a high-level neural networks API written in Python that allows for easy and flexible model construction.The architecture of the DNN model consists of fully connected layers, also called dense layers.The activation function used in the layers is ReLu, which stands for Rectified Linear Unit, a popular choice for deep learning models because it is computationally efficient and allows for faster convergence during training.The DNN model is composed of input, hidden, and output layers.The input layer has 1714 nodes, corresponding to the wavenumber points of the spectrum, while the output layer has 5 nodes, matching the vibration numbers of the GER reactants.The model also includes two hidden layers, each with 20 nodes, designed to extract high-level features from the input data.These layers enable the model to learn complex relationships and make accurate predictions.
The loss function used in the model is mean square error (MSE), which measures the average squared difference between the predicted and actual values.The optimizer used in the final model is Adam, which is a powerful optimization algorithm that adapts the learning rate based on the gradient of the loss function, allowing for more efficient training of the model.During training, the model adjusts the weights and biases in the layers to minimize the loss function and improve the accuracy of the predictions.Once the training is complete, the model can be evaluated on a test dataset to assess its performance on new data.The result shows that the accuracy of the trained model reaches 92% (Supplementary Figure 9c,d).
For a machine learning model, its generality refers to the ability to generalize patterns and knowledge from the training data to make accurate predictions or decisions on unseen or new data.
In our case, we input a new measurement data set to predict the bioreaction process.Besides, during the training process, the data is split by a ratio of 4:1 (20% is for the testing data set).The results show that the identification accuracy is 92%.Therefore, from the perspective of the machine learning field, our model exhibits good generality.When the target molecules are changed to other molecules, the refractive index of these new molecules changes but remains detectable by the antenna.For the part of machine learning algorithms, when the target molecules are changed to other molecules, the structure of the algorithm is general and still applies.We only need to retrain the machine learning model using the new data from our nanophotonic device.Overall, our methodology and model are general and not limited to specific molecules.

B. DNN model for nanorod platform
The DNN model for the nanorod platform is developed using a similar method.Notably, there are no RI-determined SPhP vibrations to quantify the individual contributions of each analyte to refractive index change, as shown in Supplementary Figure 10.Hence, it becomes challenging to decouple analyte vibrational modes when they significantly overlap with each other.This is the reason for the poor performance of the nanorod platform in the demonstration of glucose enzymatic reaction (Figure 5e       The procedure started by cleaning a BaF2 wafer in acetone using ultrasonic treatment for 10 minutes.The wafer was then rinsed in isopropanol and dried with nitrogen before undergoing 5 minutes of oxygen plasma treatment.ii: A 400 nm thick layer of PMMA e-beam lithography resist (950 PMMA A5) was spin-coated at 3000 rpm, followed by the spin-coating of a commercial electron-conducting polymer.Then it was exposed using e-beam lithography and developed using deionized water, MIBK/IPA (1:3) mixture, and isopropanol.iii-vi: SiO2/Ti/Au/Ti/SiO2/ Ti/Au layers of desired thickness were sequentially deposited using e-beam evaporation.vii: The unexposed resist was removed using a lift-off process.viii: The final SP-PhP chip.b, The nanofabrication of the nanorod platform.i: Wafer Cleaning.ii: Spin-coating PMMA e-beam lithography resists and electron-conducting polymer, followed by e-beam lithography and developing.iii: Deposit Ti/Au layers using e-beam evaporation.iv: The final nanorod chip was obtained using a lift-off process.

Supplementary Figure 2 .Supplementary Figure 3 .
Geometry-dependent characteristics of stacked nanoantennas.a, Schematic of the plasmonic stacked nanoantenna.b, Sectional view of E-field distribution of bonding mode and c, antibonding mode.d, Antenna length L and e, width W versus spectrum wavenumber mapping.f, Spectra of the device with W = 0.2 μm and 0.6 μm.g, unit cell period in the x-axis Px and h, y-axis direction Py versus spectrum wavenumber mapping.i, Spectra of the device with W = 2.8 μm and 3.8 μm.Thickness-dependent characteristics of stacked nanoantennas.a, Sectional view of the plasmonic stacked nanoantenna.W= 300 nm.L =3 μm.b, Thickness versus spectrum wavenumber mapping when the stacked Au antenna thicknesses are equal t1=t3 varying from 10 nm to 150 nm and spacer thickness t2= 50 nm, 100 nm, and 150 nm.c, The spectral response mapping when t1≠t3.t1 or t3 changes from 10 nm to 100 nm, and t2 remains unchanged.

Supplementary Figure 5 .
Incident range characteristics of stacked nanoantennas.a, Schematic of the stacked antenna and b, common mon-layer nanorod with a varying incident angle.t1= t2= t3 = 50 nm.W= 300 nm.L =3 μm.c, Incident angle θ versus spectrum wavenumber mapping of the stacked antenna.d, The performance of the mon-layer nanorod.Both the stacked antenna and the mon-layer nanorod have a varying range of 50°, in which the resonance intensity is greater than that at normal incidence.

Supplementary Figure 7 .
are coupled to form an LSPP-SPhP mode hybridization.The spectra of the LSPP-SPhP nanoantenna are shown in Supplementary Figure 7c.where we can see a dip in the SPhP reststrahlen band (gray region).When loading a molecule, the real part of the molecule refractive index Δn will cause a redshift of the plasmonic resonance (LSPP) frequency.However, the SPhP mode frequency is fixed and not affected by Δn.This will cause a change in the detuning of plasmonic resonance and SPhP mode, as shown in Supplementary Figure 7d.Then, the change in the detuning further causes the change in the spectrum of SPhP mode in the reststrahlen band ΔI (because the coupling intensity of plasmonic resonance and SPhP mode is influenced).Supplementary Figure 7e shows the obtained spectrum of SPhP mode, which is determined by the real part of the molecule refractive index Δn.So we call it refractive index-dependent SPhP vibrations in the manuscript.Apparently, the frequency redshift-induced detuning of plasmonic resonance and SPhP mode is the key to this mechanism.Hopefully, our explanation has clarified the confusion for the reviewer.The physical mechanism behind the response of SPhP mode in the reststrahlen band to refractive index.a, SPhP mode is coupled with (b) plasmonic resonance (LSPP).c, The spectrum of plasmonic resonance coupled with the SPhP mode.d, The spectra when loading a molecule with the real part of the molecule refractive index Δn.e, The extracted spectrum showing the SPhP mode.Core physical mechanism behind the response of SPhP mode in the reststrahlen band to refractive index: the refractive index-induced detuning of plasmonic resonance and SPhP vibration Fig. 2c of manuscript

Supplementary Figure 8 .
illustrates the optimization process of baseline fitting by adjusting the weight parameter P and regularization parameter μ.It is observed that P can control the position of the baseline while μ represents the similarity constraint.By iteratively updating the baseline and relaxation factor based on the equation (S17), the final baseline can be derived.Supplementary Figure 8c demonstrates the application of this method for multispectral baseline fitting and extraction of IR fingerprints of bovine serum albumin (BSA).Asymmetric least square algorithm for multispectral baseline fitting and IR fingerprints extraction.a, The effect of the weight parameter p and b, the regularization parameter β on the fitting result.c, Baseline correction for extracting IR fingerprint signals of amide I and II of protein BSA using the asymmetric least square algorithm.Regularization parameter λ=5000.

Supplementary Figure 9 .
of the main text).The DNN model incorrectly assigns the vibration signal of H2O2 to GOD because it cannot distinguish between the amide vibration of GOD and the H-O vibration of H2O2.Training results of the DNN model.a, Real-time 3D plots of differential absorbance spectra when all analytes are detected by the SP-PhP platform at the same time.b, Corresponding regression curve.c, Loss (mean square error) and d, accuracy of the DNN model in the training and validation process with the increase in the number of epochs.

2 Supplementary Figure 11 .
Calculation of refractive index using Kramers-Kronig relations.a, Original refractive index n and extinction k of glucose.b, Bare plasmonic device for SEIRAS.The resonance frequency between Device 1 and Device 2 is different.c, Simulated reflection result of SEIRAS Device 1 coupled with glucose molecules.d, Extracted differential reflectance spectra.e, Calculated n using Kramers-Kronig relations.f-h, Calculation of n using the data from another SEIRAS Device 2. i, Comparison of original n and calculated n using Kramers-Kronig relations.

Supplementary Figure 12 .Supplementary Figure 13 .
using Kramers-Kronig (device 1) n using Kramers-Kronig(device 2) Calculation of complex refractive index.a, Measured absorption spectra of glucose film.b, Real (black curve) and imaginary (red curve) parts of glucose calculated using the Drude permittivity model.c, Complex refractive index of glucose extracted from (b). d, FDTD simulation results using the calculated permittivity to verify the correctness of the calculated permittivity.Thickness-dependent transition of plasmon-phonon coupling.a, Schematic of the SP-PhP platform for investigating plasmon-phonon coupling.b, Experimental reflection spectra of the platform with SPhP layer thickness TSPhP = 30 nm and antenna length ranging from 2.2 to 4.2 μm, which corresponds to the plasmonic resonances from 900 to 1700 cm - 1 .The resonance range covers the reststrahlen band of the SPhP mode from ωTO to ωLO.c, The response of the platform when TSPhP = 70 nm and d, TSPhP = 100 nm.As observed, the strength of plasmon-phonon coupling can be controlled by changing the thickness of the SPhP antenna.

Supplementary Figure 17 .
Nanofabrication processes.a, The nanofabrication of the SP-PhP platform.i: