Introduction

Over the past few decades, metal halide perovskites have gained significant attention for potential use in solar cells1,2,3. However, their carrier mobilities are generally lower than those of conventional inorganic semiconductors, largely due to strong electron–phonon interactions4,5. This has led to the growing interest in phonon engineering within perovskites, as it can profoundly influence carrier mobility and, by extension, the energy conversion efficiency of devices.

The coherent manipulation of phonons using strong external laser fields has recently sparked considerable interest6,7,8,9,10. For instance, intense terahertz (THz) radiation can modify the band gap8 and photoluminescence spectra11 of perovskites. An alternative approach for controlling phonon properties involves cavity phonon-polaritons, which result from the coupling of phonons to the vacuum field of a cavity resonator12,13,14,15,16,17. This method bypasses the need for external light sources and the challenges associated with phonon anharmonicities. By exploiting vacuum fluctuations, the phonon properties can be engineered by adjusting the resonator geometry18. Particularly, deep subwavelength cavities make it possible to observe phonon-polaritons in nanoscale samples, with dimensions comparable to the carrier diffusion length, offering a promising strategy for mitigating carrier recombination in solar cells.

The coupling of phonons to THz subwavelength cavities also opens the door to the ultrastrong coupling (USC) regime of light–matter interaction19,20, where the coupling strength g becomes comparable to the bare mode frequencies. In this regime, the counter-rotating terms in the Hamiltonian result in the ground state becoming a squeezed vacuum21,22. Recent studies have shown that USC can give rise to intriguing phenomena, such as changes in electronic quantum transport23, tunable couplings between magnetic excitations24, and magnonic superradiant phase transitions25. A particularly intriguing aspect of the USC regime is the appearance of anomalous correlations in the polaritonic ground state, which contains both photon and matter excitations21.

A recent study on the impact of phonon-photon coupling on electron–phonon interactions in perovskites, using ultrafast pump-probe spectroscopy, showed that the mobility of photoexcited carriers remained unaffected by light–matter coupling17. This study was conducted in the strong coupling regime, where the ground state is a standard vacuum. Investigating this behavior under USC conditions remains an intriguing prospect.

Recent works have reported the observation of a single phonon mode in lead halide perovskites coupled to a THz resonator12,14,16. The possibility of coupling multiple phonon modes ultrastrongly to the resonator opens up effective avenues for modifying electron–phonon interactions in the material. In recent years, multimode light–matter coupling has gained increasing attention across various platforms26,27,28,29,30,31. Notably, multimode USC has been shown to induce ground-state correlations between different cavity modes31. However, the impact of multimode USC on intensity fluctuations in phonon-polariton systems has not yet been explored, to the best of our knowledge, despite the study of thermal photon statistics for a single two-level system in the USC regime32.

In this work, we report the observation of multimode USC between two optical phonon modes, with frequencies ω1 and ω2, in 3D MAPbI3 and 2D (BA)2MAPb2I7 lead halide perovskite crystals embedded in THz nanoslot cavities (Fig. 1a). By leveraging the small mode volume of the nanoslots, which significantly enhances the coupling strengths g1 and g2 of each phonon mode, we achieved normalized coupling strengths at resonance of gλ/ωλ ~ 0.3 (λ = 1, 2). By tuning the nanoslot resonator frequency, we observed three distinct phonon-polaritons, which exhibited two Rabi splittings in THz time-domain spectroscopy (THz-TDS) measurements. These experimental results were consistent with numerical simulations and calculations based on a Hopfield quantum model. In a polaritonic thermal state at room temperature, the model predicts the presence of “superthermal” phonon bunching in the off-resonance regime, where the nanoslot resonator frequency is much smaller than the bare phonon frequencies. For a resonator frequency ωc/(2π) = 0.1 THz, the intramode equal-time second-order correlation function \({g}_{\lambda,\lambda }^{(2)}(\tau=0)\), which quantifies the probability of simultaneous phonon emission in mode λ = 1, 2, exceeds the value \({g}_{\lambda,\lambda }^{(2)}(0)=2\) for bare phonons at thermal equilibrium and is found to be governed by the USC figure of merit gλ/ωλ. Moreover, while phonon emission in two distinct modes is uncorrelated without light–matter coupling (i.e., \({g}_{12}^{(2)}(0)=1\)), we show that multimode USC results in pronounced intermode bunching (\({g}_{12}^{(2)}(0)\approx 3\)), governed by the figure of merit g1 g2/ω1ω2.

Fig. 1: The perovskite–nanoslot hybrid system in the ultrastrong coupling regime.
figure 1

a Hybridization between a nanoslot-cavity mode, with frequency ωc, and two transverse optical phonon modes in perovskite materials, with frequencies ω1 and ω2, in the far-detuned, low-cavity-frequency regime, ωc ≪ ωλ (λ = 1, 2). The coupling strengths of these phonon modes are denoted as g1 and g2, respectively. Anomalous correlations between phonons are mediated by the cavity mode and governed by the coupling ratios gλ/ωλ. b Illustration of the perovskite–nanoslot hybrid system under illumination by terahertz light. Seven nanoslots of different lengths (l = 30–160 μm) were fabricated to tune the cavity resonance frequency. The inset shows a scanning electron microscope image showing a bare nanoslot (top view); Scale bar: 20 μm. c, Numerical simulation (COMSOL) showing an enhancement of the x component of the electric field (Ex) at resonance (0.77 THz) in a nanoslot filled with MAPbI3 perovskite. Left: top view (z = 0 plane); Right: cross-section (y = 0 plane). The white dotted lines outline the area filled with MAPbI3. The white solid lines outline the nanoslot area.

Full size image

Results

We fabricated an array of nanoslots (w = 950 nm) on quartz substrates with seven different lengths (l = 30, 40, 50, 60, 80, 120, and 160 μm) to tune the cavity mode frequency, given by \({\omega }_{{{{\rm{c}}}}}/(2\pi )={c}_{0}/(2l\sqrt{{\epsilon }_{{{{\rm{avg}}}}}})\), where c0 is the speed of light in vacuum and ϵavg = (ϵair + ϵsub)/2 represents the average dielectric constant of air and the quartz substrate (ϵsub = 2. 12)33; see Methods for sample preparation details. The resonance frequency is predominantly governed by the geometry of a single nanoslot rather than that of the periodic array34,35. Figure 1b illustrates the structure of our samples, where perovskite films (purple) are coated both on top of and within the slots.

These films exhibit two distinct optical phonon modes in free space, labeled as λ = 1 and λ = 2, corresponding to the rocking and stretching of Pb–I bonds, respectively. Due to the orientational disorder of methylammonium molecules, which breaks the lattice space-group symmetry, these phonons acquire a mixed transverse-optical (TO) and longitudinal-optical (LO) character36. As a result, they not only exhibit strong infrared absorption36,37 but also interact with lattice electrons, as recently observed5. Moreover, low-frequency phonons are particularly beneficial for achieving USC, as the normalized coupling strength g/ω increases with decreasing phonon frequency. This study, therefore, focuses on the phonons that are most relevant to strong interactions with both photonic and electronic degrees of freedom. Since electron–phonon interactions dictate charge mobility and recombination through long-range Coulomb forces, phonon-polariton formation involving low-frequency hybrid TO/LO phonons could offer effective pathways to engineer charge transport in lead halide perovskites.

Nanoslot resonators provide significant electric field enhancement due to strong optical confinement within and around the slots38,39. Since the phonon-photon coupling strength \(g\propto \sqrt{N/V}\), where N is the number of unit cells in the crystal and V is the resonator mode volume, the small mode volume of nanoslot resonators enables ultrastrong light–matter interaction regimes even with small perovskite crystals.

The in-plane spatial distribution of the cavity mode, computed using COMSOL for a perovskite-filled nanoslot, is shown in Fig. 1c (left panel). The field profile follows a sinusoidal pattern16 along the y-axis, with an electric field enhancement factor of 20 relative to transmission through a bare quartz substrate. The strong confinement of the x-component of the electric field (Ex) along the x- and z-axes results in a nearly uniform electric field within the perovskite region (Fig. 1c, right panel). Although the nanoslot thickness is 130 nm, the cavity mode extends beyond the nanoslot into the surrounding MAPbI3 layer, as depicted in Fig. 1c. The electric field strength above the nanoslot remains comparable to that inside, indicating that the perovskite film covering the slot also contributes to light–matter coupling. Notably, when t becomes comparable to the mode’s spatial extent along z, where t is the perovskite film thickness, g saturates at its maximum value; see Supplementary Note 4 and Supplementary Fig. 8.

We characterized the perovskite-nanoslot hybrid system using THz-TDS at room temperature. A normal-incident THz beam was linearly polarized along the x-axis. In free space, a 200-nm-thick MAPbI3 film exhibits transmittance dips at ω1/(2π) = 0.96 THz and ω2/(2π) = 1.9 THz, corresponding to the two phonon modes λ = 1 and λ = 2, respectively (Fig. 2a). The bare cavity resonance appears as a single peak in the transmission spectrum. Figure 2b displays the cavity resonance frequency as a function of cavity length l. By adjusting l, the cavity mode can be brought into resonance with either the λ = 1 mode (ωc = ω1) or the λ = 2 mode (ωc = ω2).

Fig. 2: Terahertz transmission spectra.
figure 2

a Transmission spectra for bare cavities (nanoslots) with different lengths l (green curves) showing a single cavity mode. The green dashed line shows the simulated transmittance through the nanoslot (l = 80 μm). Transmission spectrum for a 200-nm-thick bare MAPbI3 film (black curve) showing two transmission dips due to the two optical phonon modes (λ = 1 and λ = 2) with angular frequencies ω1 and ω2, respectively. b The bare cavity resonance frequency as a function of nanoslot length l in the reciprocal axis (green circles). The linear fit (green dashed line) shows good agreement with the experimental data. The λ = 1–cavity and λ = 2–cavity resonances occur with an 80-μm-long slot and with a 50-μm-long slot, respectively, when the cavity mode frequency coincides with the phonon frequencies (red and blue dashed lines). c Transmission spectra for the MAPbI3–nanoslots hybrid system showing three polariton branches. UP upper polariton, MP middle polariton, and LP lower polariton. The dashed lines indicate the two phonon frequencies. The spectra are vertically offset by 0.2 for clarity. d Numerical simulation (COMSOL) of the transmission as a function of cavity frequency (color map). Each spectrum has been normalized by its maximum transmittance to clearly show the three polariton branches; the black solid circles are the experimental results.

Full size image

Figure 2 c presents the transmission spectra of MAPbI3–nanoslot structures for different l values. As the nanoslots predominantly reflect incoming radiation, the observed polariton modes appear as transmission peaks. The spectra exhibit three distinct polariton branches: lower (LP), middle (MP), and upper (UP) polartions. These branches are separated by the uncoupled phonon modes λ = 1 and λ = 2 (dashed lines). As l decreases, the LP branch shifts toward λ = 1, the MP branch moves away from λ = 1 and approaches λ = 2, while the UP branch shifts away from λ = 2. Two anticrossings are observed at l = 80 μm and l = 50 μm, corresponding to ωc ≈ ω1 and ωc ≈ ω2, respectively (Fig. 2b). Due to the larger oscillator strength of the λ = 2 mode (Fig. 2a), the second Rabi splitting at l = 50 μm exceeds the first at l = 80 μm.

We carried out finite element simulations (COMSOL) to validate our experimental results, using conductivity values extracted from THz-TDS measurements (Supplementary Fig. 1) as input parameters. The simulated transmission spectra (Fig. 2d, colormap) closely match the experimental data, with black solid circles marking the resonance frequencies obtained from Fig. 2c via Lorentzian fitting. Minor discrepancies in the UP frequencies are attributed to slight shifts in the bare cavity mode (dashed green line, Fig. 2a) and additional coupling with a 3.8 THz phonon mode in the z-cut quartz substrate.

We also investigated a 2D perovskite material composed of metal halide layers separated by organic molecules, which enhances stability compared to 3D perovskites40,41,42 and holds promise for solar cell applications. Unlike 3D MAPbI3 (Fig. 3a), the presence of BA cations (CH3(CH2)3NH3) reduces the number of Pb-I bonds per unit volume, weakening the phonon mode oscillator strength. The layered structure, (BA)2(MA)n−1PbnI3n+1 (with n = 2)41, is shown in Fig. 3b. Here, n denotes the number of PbI6 octahedral layers between the BA spacer layers. The phonon modes λ = 1 and λ = 2 are slightly blueshifted compared to MAPbI3, with dips in the transmittance of a bare (BA)2MAPb2I7 200-nm-thick film at ω1/(2π) = 1.09 THz and ω2/(2π) = 2 THz, respectively (Supplementary Fig. 2). The transmission spectra of 2D perovskites embedded in nanoslot resonators resemble those of their 3D counterparts, with a larger Rabi splitting for λ = 2 due to its higher oscillator strength.

Fig. 3: Phonon-polariton properties in perovskite–nanoslot hybrid systems.
figure 3

Top: MAPbI3 films (3D perovskite). Bottom: (BA)2MAPb2I7 (2D perovskite) films. a, b Crystal structures of MAPbI3 and (BA)2MAPb2I7. BA: \({{{{\rm{CH}}}}}_{3}{({{{{\rm{CH}}}}}_{2})}_{3}{{{{\rm{NH}}}}}_{3}^{+}\), MA: \({{{{\rm{CH}}}}}_{3}{{{{\rm{NH}}}}}_{3}^{+}\). c, f Polariton dispersion as a function of cavity frequency; UP upper polariton, MP middle polariton, LP lower polariton. Solid circles: Peak frequencies extracted from the experimental transmission spectra. Solid lines: Fit of the extracted peak frequencies using the microscopic Hopfield model. The dashed lines indicate the λ = 1 and λ = 2 phonon modes and the cavity resonance. The two polariton gaps (see text) are denoted as Δ1 and Δ2. d, g Phonon Hopfield coefficients (H.C.) of the LP as a function of cavity frequency, showing a divergence in the low cavity frequency limit. e, h Theoretical predictions: Equal-time second-order phonon–phonon correlation functions \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(\tau=0)\) for a polariton thermal state at room temperature as a function of cavity frequency. The inset in (e) shows \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)\) as a function of temperature T for a cavity frequency of 0.1 THz.

Full size image

While classical electrodynamics simulations accurately reproduce the transmission spectra, we now adopt a complementary approach by utilizing a multimode Hopfield model43 to gain a deeper understanding of the ultrastrong light–matter coupling in our system and investigate its potential implications. The microscopic Hamiltonian is given by (see “Methods”)

$$\hat{H}= \hslash {\omega }_{{{{\rm{c}}}}}{\hat{a}}^{{{\dagger}} }\hat{a}+\sum\limits_{\lambda }\hslash {\omega }_{\lambda }{\hat{b}}_{\lambda }^{{{\dagger}} }{\hat{b}}_{\lambda }-i\sum\limits_{\lambda }\hslash {g}_{\lambda }\left({\hat{b}}_{\lambda }^{{{\dagger}} }-{\hat{b}}_{\lambda }\right)\left(\hat{a}+{\hat{a}}^{{{\dagger}} }\right) \\ +\sum\limits_{\lambda }\frac{\hslash {g}_{\lambda }^{2}}{{\omega }_{\lambda }}{\left(\hat{a}+{\hat{a}}^{{{\dagger}} }\right)}^{2},$$
(1)

where \({\hat{a}}^{{{\dagger}} }\) (\(\hat{a}\)) represents the creation (annihilation) operator of a cavity photon, while \({\hat{b}}_{\lambda }^{{{\dagger}} }\) (\({\hat{b}}_{\lambda }\)) denotes the creation (annihilation) operator of a phonon in the mode λ. The first two terms correspond to the bare photon and phonon Hamiltonians, respectively. The third term describes the light–matter interaction, with a coupling strength given by \({g}_{\lambda }=\frac{{\nu }_{\lambda }}{2}\sqrt{\frac{{\omega }_{\lambda }}{{\omega }_{{{{\rm{c}}}}}}}\), which is proportional to the effective ion plasma frequency νλ. The fourth term, known as the A2-term, induces a blueshift in the cavity mode frequency. The effective ion plasma frequency is determined by the effective charges associated with Pb2+ and I− ions; see “Methods” and Supplementary Note 2 for details.

The eigenfrequencies and eigenvectors of the Hamiltonian Eq. (1) are obtained via the Hopfield transformation: \({\hat{p}}_{\alpha }={\sum }_{\lambda }{X}_{\lambda,\alpha }{\hat{b}}_{\lambda }+{\sum }_{\lambda }{\widetilde{X}}_{\lambda,\alpha }{\hat{b}}_{\lambda }^{{{\dagger}} }+{Y}_{\alpha }\hat{a}+{\widetilde{Y}}_{\alpha }{\hat{a}}^{{{\dagger}} }\), where \({\hat{p}}_{\alpha }\) is the annihilation operator of a polariton in the mode α = LP, MP, UP, with frequency ωα. Up to a constant term, Eq. (1) can then be expressed in its diagonal form as \(\hat{H}={\sum }_{\alpha }\hslash {\omega }_{\alpha }{\hat{p}}_{\alpha }^{{{\dagger}} }{\hat{p}}_{\alpha }\). The system enters the USC regime when the normalized coupling strength at resonance satisfies gλ/ωλ = νλ/2ωλ ≳ 0.1. In this regime, the counter-rotating terms \(\propto {\hat{b}}_{\lambda }\hat{a},{\hat{b}}_{\lambda }^{{{\dagger}} }{\hat{a}}^{{{\dagger}} }\) in the Hamiltonian Eq. (1), along with the anomalous Hopfield coefficients \({\widetilde{Y}}_{\alpha }\) and \({\widetilde{X}}_{\lambda,\alpha }\), play a significant role.

Figure 3c presents the polariton dispersion for the MAPbI3–nanoslots system. The coupling strengths gλ are extracted by fitting the peak frequencies (solid circles) of the transmission spectra to the calculated eigenfrequencies ωα (solid lines). When the nanoslot resonator is resonant with the phonon modes λ = 1 and λ = 2, we obtain normalized coupling strengths of g1/ω1 = 0.28 (ωc = ω1) and g2/ω2 = 0.3 (ωc = ω2), respectively. These values confirm that both phonon modes are in the USC regime with the nanoslot resonator. The corresponding Rabi splittings at the two resonances are 0.45 THz and 1.13 THz. Notably, while the Rabi splitting equals exactly 2g for a single matter and cavity mode in the strong coupling regime, this relation breaks down in the USC regime due to counter-rotating terms. In our case, the inclusion of two phonon modes leads to further deviations from the 2g value. Importantly, the polariton dispersion should be understood as the result of the simultaneous coupling of both phonon modes to the cavity mode, with all three degrees of freedom treated on equal footing. This becomes evident when examining the contribution of the two phonon modes to the MP mode at around the resonance between the λ = 1 phonon and the cavity mode. As shown in Supplementary Fig. 4b, the contribution \({W}_{\lambda }^{{{{\rm{MP}}}}}=| {X}_{\lambda,{{{\rm{MP}}}}}{| }^{2}-| {\widetilde{X}}_{\lambda,{{{\rm{MP}}}}}{| }^{2}\) of the two phonon modes (λ = 1, 2) to the MP mode is indeed of comparable magnitude, indicating that the MP branch involves significant hybridization with both phonon modes.

In the USC regime, distinctive features appear not only near resonance, but also when the resonator frequency is much lower than the phonon frequencies, ωc ≪ ωλ. Unlike in the strong coupling regime, the polariton modes do not converge to the uncoupled mode frequencies when ωc ≪ ωλ. This behavior defines the so-called polariton gaps19,44, given by \({\Delta }_{1}={\lim }_{{\omega }_{{{{\rm{c}}}}}\to 0}{\omega }_{{{{\rm{MP}}}}}-{\omega }_{1}\) and \({\Delta }_{2}={\lim }_{{\omega }_{{{{\rm{c}}}}}\to 0}{\omega }_{{{{\rm{UP}}}}}-{\omega }_{2}\), as shown in Fig. 3c. As detailed in Supplementary Note 3, the frequencies of the MP and UP modes, ωMP and ωUP, asymptotically approach \({\widetilde{\omega }}_{1}=\sqrt{{\omega }_{1}^{2}+{\nu }_{1}^{2}}\) and \({\widetilde{\omega }}_{2}=\sqrt{{\omega }_{2}^{2}+{\nu }_{2}^{2}}\), respectively, in the low-resonator frequency limit ωc → 0.

This unconventional behavior is linked to strong light–matter hybridization, which persists even in this far-detuned, low-resonator-frequency regime. This is reflected in the divergence of the light–matter coupling strength, \({g}_{\lambda }\propto \sqrt{1/{\omega }_{{{{\rm{c}}}}}}\), and the A2-term, which scales as \({g}_{\lambda }^{2}\), as ωc → 0. In this regime, the MP (UP) mode is mainly a hybrid of the phonon mode λ = 1 (λ = 2) and cavity photons. The corresponding normal (Yα) and anomalous (\({\widetilde{Y}}_{\alpha }\)) Hopfield coefficients become large and comparable, scaling as \({Y}_{{{{\rm{MP}}}}} \sim {\widetilde{Y}}_{{{{\rm{MP}}}}} \sim {\nu }_{1}/\sqrt{{\omega }_{{{{\rm{c}}}}}{\widetilde{\omega }}_{1}}\) and \({Y}_{{{{\rm{UP}}}}} \sim {\widetilde{Y}}_{{{{\rm{UP}}}}} \sim {\nu }_{2}/\sqrt{{\omega }_{{{{\rm{c}}}}}{\widetilde{\omega }}_{2}}\). In contrast, the LP mode mixes the cavity field with both phonon modes. The phonon contributions to this polariton, quantified by the coefficients Xλ,LP and \({\widetilde{X}}_{\lambda,{{{\rm{LP}}}}}\), also grow large and comparable, with \({X}_{\lambda,{{{\rm{LP}}}}} \sim {\widetilde{X}}_{\lambda,{{{\rm{LP}}}}} \sim {\nu }_{\lambda }/\sqrt{{\omega }_{{{{\rm{c}}}}}{\omega }_{\lambda }}\). These coefficients are shown in Fig. 3d, while the other Hopfield coefficients are provided in Supplementary Figs. 3 and 5.

Due to the large anomalous Hopfield coefficients, the polaritonic ground state \(\left\vert G\right\rangle\), defined by \({\prod }_{\alpha }{\hat{p}}_{\alpha }\left\vert G\right\rangle=0\), takes the form of a multimode squeezed vacuum in the low resonator frequency regime. This state contains correlated photon pairs, contributed by the MP and UP, as well as intermode and intramode phonon pairs originating from the LP. This multimode squeezed vacuum exhibits strong entanglement between the two phonon modes, as discussed below.

It is important to highlight that while both the normal and anomalous Hopfield coefficients diverge in the limit ωc → 0, the total phonon and photon weights remain finite due to the normalization of the Hopfield coefficients (see Supplementary Figs. 4 and 6). Moreover, we stress that the low-resonator-frequency regime (long cavity) does not imply the absence of a cavity, as its transverse confinement remains deeply subwavelength.

A distinctive feature of the multimode USC regime is the presence of anomalous correlations between the phonon modes. By inverting the Hopfield transformation, one obtains the correlation functions:

$$\langle {\hat{b}}_{\lambda }^{{{\dagger}} }{\hat{b}}_{{\lambda }^{{\prime} }}\rangle=\sum\limits_{\alpha }{\left({\widetilde{X}}_{\lambda }^{\alpha }\right)}^{*}{\widetilde{X}}_{{\lambda }^{{\prime} }}^{\alpha }(1+{n}_{\alpha })+\sum\limits_{\alpha }{X}_{\lambda }^{\alpha }{\left({X}_{{\lambda }^{{\prime} }}^{\alpha }\right)}^{*}{n}_{\alpha },$$
(2)
$$\langle {\hat{b}}_{\lambda }{\hat{b}}_{{\lambda }^{{\prime} }}\rangle=-\sum\limits_{\alpha }{\left({X}_{\lambda }^{\alpha }\right)}^{*}{\widetilde{X}}_{{\lambda }^{{\prime} }}^{\alpha }(1+{n}_{\alpha })-\sum\limits_{\alpha }{\left({X}_{{\lambda }^{{\prime} }}^{\alpha }\right)}^{*}{\widetilde{X}}_{\lambda }^{\alpha }{n}_{\alpha }.$$
(3)

Here, \({n}_{\alpha }=\langle {\hat{p}}_{\alpha }^{{{\dagger}} }{\hat{p}}_{\alpha }\rangle\) represents the population in the polariton mode α. In the polaritonic ground state (nα = 0), Eqs. (2) and (3) show that such correlations arise only when the anomalous Hopfield coefficients \({\widetilde{X}}_{\lambda }^{\alpha }\) are nonzero. Moreover, these correlations are further enhanced in excited polariton states where nα ≠ 0. To explore the impact of multimode USC on correlated phonon emission, we consider the second-order correlation function45, which quantifies the joint probability of a phonon being emitted in the mode \({\lambda }^{{\prime} }\) at time t + τ given that a phonon was emitted in the mode λ at time t:

$${g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(\tau )=\frac{\langle {\hat{b}}_{{\lambda }^{{\prime} }}(t+\tau ){\hat{b}}_{\lambda }(t){\hat{b}}_{\lambda }^{{{\dagger}} }(t){\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }(t+\tau )\rangle }{\langle {\hat{b}}_{\lambda }(t){\hat{b}}_{\lambda }^{{{\dagger}} }(t)\rangle \langle {\hat{b}}_{{\lambda }^{{\prime} }}(t+\tau ){\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }(t+\tau )\rangle }.$$
(4)

Although classical wave-based methods could, in principle, be used to study intensity correlations, such calculations are technically challenging and non-trivial to implement. For this reason, we focus in this work on quantum predictions for the second-order correlation functions. Assuming a thermal polariton state at temperature T, with \({n}_{\alpha }={\left({e}^{\hslash {\omega }_{\alpha }/{k}_{{{{\rm{B}}}}}T}-1\right)}^{-1}\) and kB the Boltzmann constant, the equal-time intramode (\(\lambda={\lambda }^{{\prime} }\)) and intermode (\(\lambda \, \ne \, {\lambda }^{{\prime} }\)) correlation functions are given by

$${g}_{\lambda,\lambda }^{(2)}(0)=2+\frac{\langle {\hat{b}}_{\lambda }{\hat{b}}_{\lambda }\rangle \langle {\hat{b}}_{\lambda }^{{{\dagger}} }{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle }{{\langle {\hat{b}}_{\lambda }{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle }^{2}},$$
(5)
$${g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)=1+\frac{\langle {\hat{b}}_{\lambda }{\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }\rangle \langle {\hat{b}}_{{\lambda }^{{\prime} }}{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle }{\langle {\hat{b}}_{\lambda }{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle \langle {\hat{b}}_{{\lambda }^{{\prime} }}{\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }\rangle }+\frac{\langle {\hat{b}}_{\lambda }{\hat{b}}_{{\lambda }^{{\prime} }}\rangle \langle {\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle }{\langle {\hat{b}}_{\lambda }{\hat{b}}_{\lambda }^{{{\dagger}} }\rangle \langle {\hat{b}}_{{\lambda }^{{\prime} }}{\hat{b}}_{{\lambda }^{{\prime} }}^{{{\dagger}} }\rangle },$$
(6)

respectively. For vanishing anomalous Hopfield coefficients (\({\widetilde{X}}_{\lambda }^{\alpha }=0\,\forall \lambda\)) or in the absence of phonon–photon coupling (νλ = 0), Eq. (5) simplifies to \({g}_{\lambda,\lambda }^{(2)}(0)=2\), which corresponds to intramode phonon bunching-a hallmark of thermal states. In contrast, the intermode correlation function satisfies \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)=1\) for \(\lambda \, \ne \, {\lambda }^{{\prime} }\), indicating that intermode phonon emission remains uncorrelated and follows Poissonian statistics46.

In the multimode USC regime, we predict a significant modification of equal-time phonon-phonon second-order correlations. As shown in Fig. 3e, the various contributions \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)\) to a value of 3 in the limit of a vanishing resonator frequency (ωc → 0). As ωc increases, \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)\) decreases monotonically, approaching 2 for intramode correlations (\(\lambda={\lambda }^{{\prime} }\)) and 1 for intermode correlations (\(\lambda \ne {\lambda }^{{\prime} }\)). These limiting values correspond to the correlations expected for bare phonons, which are recovered in the high resonator frequency regime (ωc ≫ ωλ), where the LP and MP asymptotically approach the uncoupled phonon frequencies ω1 and ω2, respectively.

For a detuned cavity with ωc/(2π) = 0.1 THz at room temperature (T = 300 K), our theoretical model predicts \({g}_{1,1}^{(2)}(0)\approx 2.86\), \({g}_{2,2}^{(2)}(0)\approx 2.96\), and \({g}_{1,2}^{(2)}(0)\approx 2.82\). These results indicate that multimode USC should lead to strong phonon bunching in both intramode and intermode correlations. This effect primarily arises from the LP, which exhibits large normal and anomalous phonon Hopfield coefficients (\({X}_{\lambda }^{\alpha },{\widetilde{X}}_{\lambda }^{\alpha }\)), as discussed earlier, along with a significant population in the low resonator frequency regime (nLP ≈ 80 for ωc/(2π) = 0.1 THz and T = 300 K). The inset of Fig. 3e illustrates the temperature dependence of the calculated second-order phonon correlations, showing that at T = 0 K, intramode and intermode correlations are enhanced by approximately 10% and 40%, respectively, compared to the bare phonon case. Notably, at room temperature, \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)\) remains in the saturation regime.

Figure 3f presents the extracted peak frequencies for the 2D perovskite (BA)2MAPb2I7–nanoslots system, alongside theoretical predictions (solid lines). In this case, the λ = 1 mode exhibits a very small polaritonic gap, which is consistent with the normalized coupling strength \({g}_{1}^{{\prime} }/{\omega }_{1}=0.13\) extracted from the fit. This value suggests that λ = 1 is on the verge of the USC regime, leading to reduced Hopfield coefficients X1,LP and \({\widetilde{X}}_{{{{\rm{1,LP}}}}}\) compared to the MAPbI3–nanoslots system (see Fig. 3g). Conversely, the polaritonic gap of the λ = 2 mode remains similar to that observed in the MAPbI3–nanoslots system, consistent with the large coupling ratio \({g}_{2}^{{\prime} }/{\omega }_{2}=0.23\) extracted from the fit.

At ωc/(2π) = 0.1 THz and room temperature, the calculated phonon-phonon correlations for \(\lambda={\lambda }^{{\prime} }=1\) and \(\lambda=1,{\lambda }^{{\prime} }=2\) are significantly lower than in the MAPbI3-nanoslots system (\({g}_{1,1}^{(2)}(0)\approx 2.61\), \({g}_{1,2}^{(2)}(0)\approx 2.52\)), in line with the weaker coupling strength of λ = 1. However, \({g}_{2,2}^{(2)}(0)\approx 2.94\) remains nearly unchanged compared to the 3D perovskite system, as shown in Fig. 3h.

Using a perturbative expansion valid for ωc/ωλ ≪ 1 and νλ/ωλ ≪ 1, we show in Supplementary Note 3 that the second-order correlation functions can be approximated as

$${g}_{1,1}^{(2)}(0)\approx 2+{\left(\frac{{g}_{1}}{{\omega }_{1}}\right)}^{4}{\left(\frac{1+2{n}_{{{{\rm{LP}}}}}}{1+{n}_{{{{\rm{MP}}}}}}\right)}^{2}$$
(7a)
$${g}_{2,2}^{(2)}(0)\approx 2+{\left(\frac{{g}_{2}}{{\omega }_{2}}\right)}^{4}{\left(\frac{1+2{n}_{{{{\rm{LP}}}}}}{1+{n}_{{{{\rm{UP}}}}}}\right)}^{2}$$
(7b)
$${g}_{1,2}^{(2)}(0)\approx 1+2{\left(\frac{{g}_{1}}{{\omega }_{1}}\right)}^{2}{\left(\frac{{g}_{2}}{{\omega }_{2}}\right)}^{2}\frac{{(1+2{n}_{{{{\rm{LP}}}}})}^{2}}{(1+{n}_{{{{\rm{MP}}}}})(1+{n}_{{{{\rm{UP}}}}})}.$$
(7c)

These results show that the intramode correlation functions are primarily controlled by the standard USC figure of merit, gλ/ωλ. In contrast, intermode correlations are governed by the product g1g2/ω1ω2, which becomes an important figure of merit for multimode USC. For instance, we find g1g2/ω1ω2 = 0.084 in the 3D perovskite-nanoslot system and g1g2/ω1ω2 = 0.03 in the 2D system. The specific form g1g2/ω1ω2 suggests that intermode correlations arise from the effective coupling between phonons mediated by the far-detuned cavity, where ωc ≪ ωλ.

Discussion

We report the observation of cavity phonon-polaritons in the multimode USC regime. The light–matter coupling strength was controlled by tuning the number of PbI6 octahedral layers between the BA spacer layers of the perovskite, directly influencing the phonon oscillator strengths. Unlike recent studies on multimode USC, which have primarily focused on engineering photonic properties in THz cavities through coupling with inorganic semiconductor (e.g., GaAs) quantum wells29,30,31, our complementary approach leverages a deep-subwavelength cavity resonator to mediate effective coupling between matter excitations. This approach has the potential to modify fundamental material properties, such as charge carrier mobilities. Given the relevance to solar cell applications, we focus on multimode USC of phonons in lead halide perovskite thin films, which are known to exhibit strong electron–phonon interactions5,47,48. Our cavity-mediated phonon–phonon coupling mechanism provides an effective route for controlling phonon–phonon correlations at thermal equilibrium, without requiring external driving fields or phonon anharmonicities.

The small mode volume of the nanoslots enabled USC with the highest resonant coupling strengths reported in cavity phonon-polariton systems. The use of deep-subwavelength resonators filled with lead halide perovskite films of a few hundred nanometers in thickness—comparable to the carrier diffusion length—is fully compatible with solar cell applications49. In the off-resonance regime, where the cavity frequency is much lower than the phonon frequencies, the coupling strength scales as \({g}_{\lambda }\propto 1/\sqrt{{\omega }_{{{{\rm{c}}}}}}\), allowing access to a unique regime where counter-rotating terms in the Hamiltonian become as significant as the rotating-wave terms. This leads to anomalous correlations governed by the USC figure of merit gλ/ωλ at thermal equilibrium, even in the absence of nonlinear interactions. We demonstrate theoretically that in this regime, the cavity mode mediates an effective interaction between the two phonon modes λ and \({\lambda }^{{\prime} }\), resulting in superthermal intermode phonon bunching \(\propto {g}_{\lambda }{g}_{{\lambda }^{{\prime} }}/\sqrt{{\omega }_{\lambda }{\omega }_{{\lambda }^{{\prime} }}}\). This corresponds to the correlated emission of phonons, characterized by an equal-time second-order phonon–phonon correlation function \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0) > 2\). In contrast, for bare phonons in thin films without a cavity, phonon emission in different modes remains uncorrelated, i.e., \({g}_{\lambda,{\lambda }^{{\prime} }}^{(2)}(0)=1\).

Although directly measuring phonon–phonon correlations remains challenging, indirect evidence can be obtained through quantum optics techniques that measure the second-order photon–photon correlation function. As mentioned earlier, this function also exhibits superthermal bunching and can be directly measured using femtosecond noise correlation spectroscopy, a method previously used to study mode fluctuations in quantum fields50. Additionally, a recent study on magnons demonstrated that statistical correlations between two probe pulses can be used to extract fluctuations of collective excitations51, a technique that can be adapted to phonons. Finally, nonlinear spectroscopic techniques, such as two-dimensional THz spectroscopy, can provide further indirect insights into phonon–phonon correlations.

Compared to recent studies on single-mode phonon-polaritons in similar systems12,14,16, our multimode scenario presents opportunities for controlling electron–phonon interactions in lead halide perovskites, with implications for light-harvesting and light-emitting devices. This advantage stems from the unique nature of the phonons studied here: they are strongly coupled both to the cavity and to charge carriers due to their mixed TO/LO character. Furthermore, we show that by employing sufficiently long resonators, it is possible to achieve a high phonon-to-cavity frequency ratio, ω/ωc, which effectively compensates for the small coupling ratio, ν/ω, of higher-energy phonons with a predominantly LO character and consequently low oscillator strength. Previous studies have demonstrated that charge carrier scattering, mediated by the Fröhlich interaction with LO phonons near 3 THz, dominates electron–phonon coupling in these materials at room temperature4. In long nanoslots, this compensation effect would enhance the coupling strength of these LO phonons, \(g/\omega=(\nu /\omega )\sqrt{\omega /{\omega }_{{{{\rm{c}}}}}}\), potentially yielding a large multimode USC figure of merit and, therefore, strong intermode phonon bunching with the low-frequency phonons λ = 1, 2.

Our approach thus enables control over high-frequency LO phonons via USC coupling to low-frequency IR-active phonons in long cavities, potentially leading to substantial modifications in electron–phonon scattering. This motivates further pump-probe photoconductivity experiments under multimode USC conditions to explore the modulation of photoexcited carrier mobility through electron–phonon interactions in perovskite solar cells.

More broadly, our findings open effective directions for phonon-based quantum technologies52,53,54 in nonequilibrium scenarios, with applications ranging from the control of superconductivity55 and multimode entanglement56 to the generation of coherent THz sources57 and enhanced energy transfer58 in solid-state systems. Notably, while optical phonons typically do not contribute to heat transfer due to their vanishing group velocity, in our system, the LP acquires a finite group velocity and retains a substantial (~20%) phonon weight in the low-cavity-frequency regime. This suggests the intriguing possibility that superthermal phonon bunching in the multimode USC regime may influence heat transport in perovskite materials.

Methods

Sample preparation

Dimethyl Formamide (DMF), Dimethyl Sulfoxide (DMSO), Lead Oxide (PbO), Butylamine (BA), Hydriodic acid (HI), Hypophosphorous acid (H3PO2), and Diethyl Ether were purchased from Sigma Aldrich and used without any further treatment. Methylammonium Iodide (MAI) and Methylammonium Chloride (MACl) were purchased from Greatcell Solar. Lead Iodide (PbI2) was purchased from TCI Chemicals.

For the synthesis of MAPbI3 (3D) films, the precursor solution was made by dissolving 95.4 mg MAI, 276.6 mg PbI2 in 638 μl DMF and 71 μl DMSO, and stirred on the hotplate at 70 °C for 3 hours. 4 mg MACl was added to improve the film crystallinity. 70 μl solution was then spin-coated on the nanoslots at 5000 rpm and 3500 rpm/s acceleration for 30 seconds. 600 μl of Diethyl ether was dripped at 10 seconds from the start. The films were then annealed at 100 °C for 10 minutes.

(BA)2MAPb2I7 (2D) crystals (in the form of small plates) were prepared by the following procedure59. The parent crystals were dissolved in DMF at a concentration of 0.2 M (30 mg/100 μl) and stirred on the hotplate at 70 °C for 2 hours in an argon glovebox. (BA)2MAPb2I7 (2D) films were made with the phase-selective method60. The solution was then transferred to a different glovebox where it was spin-coated on the nanoslots at 5000 rpm with 3500 rpm/s acceleration for 30 seconds. The films turn red-brown during the spin-coating process and are annealed at 100 °C for 5 minutes.

To fabricate the nanoslots, we utilized a standard photolithography technique to pattern photoresists to form an array of rods (950 nm by l), followed by Au deposition (150 nm) by an electron beam evaporator. Then, we performed an Ar beam ion milling on the samples to facilitate a lift-off process. In this process, the thickness of the Au films decreased to 130 nm by the ion milling. After the lift-off process with acetone, we obtained an array of nanoslots. An array of fabricated bare nanoslots is presented in a scanning electron microscope image (top view) in the inset of Fig. 1b. Then, perovskite polycrystalline films ( ~ 200 nm thick) were coated on the nanoslots.

THz time-domain spectroscopy

We performed THz-TDS transmission measurements in a dry air environment at room temperature. The total measurement time for each sample was less than 15 minutes to avoid the degradation of perovskite films. To access high-frequency THz emission (up to 3 THz), we utilized InGaAs photoconductive antennas for both emitter and detector, which are fiber-coupled with an Er-fiber laser (80 MHz, 1.5 μm). Electric field amplitudes were low enough to avoid any field strength-dependent nonlinear effects. The emitted THz waves are guided to be sequentially focused on the samples and the detector by four 90°-off off-axis parabolic mirrors. The THz beam size at the focal point was about 1 mm. To obtain transmission spectra of samples \(\tilde{T}={\left\vert {E}_{{{{\rm{sample}}}}}(\omega )/{E}_{{{{\rm{ref}}}}}(\omega )\right\vert }^{2}\), we first measured transmitted electric fields Esample(t) of a sample and those of a bare quartz substrate Eref(t) as a reference, where t is a delay time. Then, we performed the Fourier transformation to obtain Esample(ω) and Eref(ω).

Microscopic model

The microscopic phonon-polariton model is an extension of the Hopfield model to the multimode regime. Its physical validity is supported by a rigorous comparison between two equivalent formulations of cavity quantum electrodynamics, as detailed below.

The ionic vibrations within a unit cell of the perovskite material are modeled as effective spherical ions, each characterized by a reduced mass Mj and an effective charge Zj, where j = 1, 2,⋯   labels the effective ions in the unit cell of volume a3. The reduced mass Mj captures the relative motion of all atoms involved in a given optical phonon mode, thus forming the effective ion j61. The perovskite is assumed to occupy a volume V = w(h + t)l. The electromagnetic field confined within the nanoslots is considered to be spatially uniform along the x and z directions, resulting in a mode volume equal to V.

Within the minimal coupling framework, the system Hamiltonian is expressed as

$$\hat{H}= \, \frac{{\epsilon }_{0}}{2}\int\,d{{{\bf{r}}}}\left[\epsilon {\hat{{{{\bf{E}}}}}}^{2}({{{\bf{r}}}})+{c}^{2}{\left({{{\boldsymbol{\nabla }}}}\times \hat{{{{\bf{A}}}}}({{{\bf{r}}}})\right)}^{2}\right]\\ +\sum\limits_{j}\frac{1}{2{M}_{j}}\int\frac{d{{{\bf{r}}}}}{{a}^{3}}{\hat{{{{\boldsymbol{{{{\mathcal{P}}}}}}}}}}_{j}^{2}({{{\bf{r}}}})+\frac{{M}_{j}{\omega }_{\lambda }^{2}}{2}\int\frac{d{{{\bf{r}}}}}{{a}^{3}}{\hat{{{{\boldsymbol{{{{\mathcal{Q}}}}}}}}}}_{j}^{2}({{{\bf{r}}}})\\ -\sum\limits_{j}\frac{{Z}_{j}}{{M}_{j}}\int\frac{d{{{\bf{r}}}}}{{a}^{3}}{\hat{{{{\boldsymbol{{{{\mathcal{P}}}}}}}}}}_{j}({{{\bf{r}}}})\cdot \hat{{{{\bf{A}}}}}({{{\bf{r}}}})+\sum\limits_{j}\frac{{Z}_{j}^{2}}{2{M}_{j}}\int\frac{d{{{\bf{r}}}}}{{a}^{3}}{\hat{{{{\bf{A}}}}}}^{2}({{{\bf{r}}}}).$$
(8)

Restricting the model to the fundamental cavity mode, the vector potential operator can be written as

$$\hat{{{{\bf{A}}}}}({{{\bf{r}}}})=\sqrt{\frac{\hslash }{{\epsilon }_{0}\epsilon {\omega }_{{{{\rm{c}}}}}V}}\sin \left(\frac{\pi y}{l}\right)\left(\hat{a}+{\hat{a}}^{{{\dagger}} }\right){{{\bf{v}}}},$$
(9)

with an analogous expansion for the electric field \(\hat{{{{\bf{E}}}}}\). Here, ϵ0 is the vacuum permittivity, ϵ the background dielectric constant of the perovskite, ωc the frequency of the fundamental cavity mode, and v the unit polarization vector. The field is quantized via the bosonic operators \(\hat{a}\) and \({\hat{a}}^{{{\dagger}} }\).

Imposing strict boundary conditions in all spatial directions, the displacement and momentum fields of ion j, \({\hat{{{{\boldsymbol{{{{\mathcal{Q}}}}}}}}}}_{j}({{{\bf{r}}}})\) and \({\hat{{{{\boldsymbol{{{{\mathcal{P}}}}}}}}}}_{j}({{{\bf{r}}}})\), are expanded in Fourier series as follows:

$${\hat{{{{\boldsymbol{{{{\mathcal{Q}}}}}}}}}}_{j}({{{\bf{r}}}})=\sum\limits_{\lambda,n,q}\sqrt{\frac{4\hslash {a}^{3}}{{M}_{j}{\omega }_{\lambda }V}}\, {\phi }_{n,q}({{{\bf{r}}}})\left({\hat{b}}_{nq,\lambda }^{{{\dagger}} }+{\hat{b}}_{nq,\lambda }\right){{{{\bf{u}}}}}_{\lambda,j}$$
(10)
$${\hat{{{{\boldsymbol{{{{\mathcal{P}}}}}}}}}}_{j}({{{\bf{r}}}})=i\sum\limits_{\lambda,n,q}\sqrt{\frac{4\hslash {M}_{j}{\omega }_{\lambda }{a}^{3}}{V}}\, {\phi }_{n,q}({{{\bf{r}}}})\left({\hat{b}}_{nq,\lambda }^{{{\dagger}} }-{\hat{b}}_{nq,\lambda }\right){{{{\bf{u}}}}}_{\lambda,j},$$
(11)

where \({\phi }_{n,q}({{{\bf{r}}}})=\sin \left(\frac{\pi nx}{w}\right)\sin \left(\frac{\pi y}{l}\right)\sin \left(\frac{\pi qz}{h}\right)\), and \(n,q\in {\mathbb{N}}\). The operators \({\hat{b}}_{nq,\lambda }\) and \({\hat{b}}_{nq,\lambda }^{{{\dagger}} }\) annihilate and create a TO phonon in branch λ, of frequency ωλ. The polarization vectors uλ,j are real, mode-independent, and the phonons are considered dispersionless.

Introducing collective phonon operators \({\hat{b}}_{\lambda }\) and \({\hat{b}}_{\lambda }^{{{\dagger}} }\) as coherent superpositions of phonon modes (see Supplementary Note 2), the minimal-coupling Hamiltonian in Eq. (8) is obtained as:

$$\hat{H}= \, \hslash {\omega }_{{{{\rm{c}}}}}{\hat{a}}^{{{\dagger}} }\hat{a}+\sum\limits_{\lambda }\hslash {\omega }_{\lambda }{\hat{b}}_{\lambda }^{{{\dagger}} }{\hat{b}}_{\lambda }-i\sum\limits_{\lambda }\frac{\hslash {\nu }_{\lambda }}{2}\sqrt{\frac{{\omega }_{\lambda }}{{\omega }_{{{{\rm{c}}}}}}}\left({\hat{b}}_{\lambda }^{{{\dagger}} }-{\hat{b}}_{\lambda }\right)\left(\hat{a}+{\hat{a}}^{{{\dagger}} }\right)\\ +\frac{{\sum }_{j}\hslash {f}_{j}^{2}}{4{\omega }_{{{{\rm{c}}}}}}{\left(\hat{a}+{\hat{a}}^{{{\dagger}} }\right)}^{2},$$
(12)

where \({f}_{j}=\sqrt{{Z}_{j}^{2}/({\epsilon }_{0}\epsilon {M}_{j}{a}^{3})}\) denotes the ionic plasma frequency of ion j. The coefficient of the A2-term can be directly related to the effective plasma frequency νλ that governs the light–matter coupling strength. This relation is obtained by reformulating the Hamiltonian in the Power–Zienau–Woolley (PZW) gauge and comparing it with the transformed minimal coupling Hamiltonian, as described in the Supplementary Note 2. The equivalence \({\sum }_{j}{f}_{j}^{2}={\sum }_{\lambda }{\nu }_{\lambda }^{2}\) results from this procedure. Substituting this identity into Eq. (12) yields Hamiltonian in Eq. (1), which is diagonalized using a Hopfield–Bogoliubov transformation, as explained in the main text.

The approximate second-order correlation functions defined in Eq. (7) are computed by applying a Schrieffer–Wolff transformation to the PZW Hamiltonian, which enables the adiabatic elimination of the cavity field. This yields an effective phonon–phonon interaction mediated by virtual photons. The derivation is presented in detail in the Supplementary Note 3.