Disentangling collective coupling in vibrational polaritons with double quantum coherence spectroscopy

Abstract

Vibrational polaritons are formed by strong coupling of molecular vibrations and photon modes in an optical cavity. Experiments have demonstrated that vibrational strong coupling can change molecular properties and even affect chemical reactivity. However, the interactions in a molecular ensemble are complex, and the exact mechanisms that lead to modifications are not fully understood yet. We simulate two-dimensional infrared spectra of molecular vibrational polaritons based on the double quantum coherence technique to gain further insight into the complex many-body structure of these hybrid light-matter states. Double quantum coherence uniquely resolves the excitation of hybrid light-matter polaritons and allows to directly probe the anharmonicities of the resulting states. By combining the cavity Born-Oppenheimer Hartree-Fock ansatz with a full quantum dynamics simulation of the corresponding eigenstates, we go beyond simplified model systems. This allows us to study the influence of self-polarization and the response of the electronic structure to the cavity interaction on the spectral features even beyond the single-molecule case.

I Introduction

Light-matter coupling between optical resonances of a cavity and molecular transitions can result in the formation of molecular polaritons [1, 2, 3, 4]. When coupling overcomes the dissipative processes, the strong coupling regime is reached, and these new hybridized states with mixed photon-matter character can be observed spectroscopically. A pair of characteristic peaks, called the lower polariton (LP) and upper polariton (UP), can be observed shifted above and below the field-free transition frequency. By controlling of the photonic environment it becomes possible to selectively couple the cavity photon modes to vibrational or electronic transitions in molecules, called vibrational-strong coupling (VSC) and electronic-strong coupling (ESC), respectively. Both types of strong coupling are discussed as effective tools not only for modifying photophysics and photochemistry [5, 6, 7, 8] but also for altering electronic ground-state reaction rates and product branching ratios [9, 10, 11]. However, the current theoretical description of VSC in particular is far from complete and many open questions remain. For example, it is not clear how the polariton states contribute to the experimentally observed modification of reaction rates. The formation of polaritons in an ensemble is inherently a collective phenomenon in which the excitation is delocalized over the whole system. However, a chemical reaction is thought to occur locally in individual molecules.

Consequently, both linear [12, 13, 14] and non-linear spectroscopy [15, 16, 17, 18, 19, 20, 21, 22] are methods that can be used to advance the understanding of molecular polaritonics. In particular, coherent multidimensional infrared spectroscopy is a powerful tool that can probe anharmonicities and vibrational energy relaxation pathways, opening the possibility of further insight into processes under VSC. Among the multipulse non-linear spectroscopic techniques available [23] we focus on the double quantum coherence (DQC) technique [24, 25, 26, 27, 22] in this work. The study of DQC for molecular vibrational polaritons in optical cavities allows one to directly probe the double excitation manifold without being convolved with single exciton resonances. Furthermore, the effect of strong couplings on the vibrational anharmonicities in the systems can be observed.

In this manuscript we simulate DQC spectra for a single HF molecule and a pair of HF molecules resonantly coupled to a single-photon mode of an optical cavity in the infrared. Our simulation is based on full-dimensional cavity potential energy surfaces (cPESs) calculated on the cavity Born-Oppenheimer Hartree-Fock (CBO-HF) [28, 29] level of theory. The CBO-HF ansatz allows the molecular system, e.g., the electronic structure, to respond to the cavity field due to the self-consistent field (SCF) procedure, and the inclusion of dipole self-energy (DSE) terms allows self-polarization [30, 31, 28, 29] of the coupled cavity-molecular system. In addition to a detailed analysis of the obtained DQC spectra, we also investigated the influence of the DSE, e.g. self-polarization, and the SCF procedure, e.g. response of the electronic structure, on the spectral features. Finally, we study the changes introduced by going beyond a single molecule and the possible influence of cavity-mediated intermolecular interactions on the DQC signal.

II Theory and models

II.1 Cavity Born-Oppenheimer approximation

In our recent work [28, 32] we have introduced the CBO-HF approach, which represents a formulation of the well-known Hartree-Fock ansatz in the context of cavity Born-Oppenheimer approximation (CBOA). Within CBOA, the cavity field modes are grouped with the nuclei in a generalized Born-Huang expansion [33, 34], and then one can subsequently solve the quantum problem of the electrons and then of the nuclei and photons. The electronic CBOA Hamiltonian for a single cavity field mode takes the form of

$$\hat{H}_{CBO}=\hat{H}_{el}+\frac{1}{2}\omega_{c}^{2}q_{c}^{2}-\omega_{c}q_{c}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)+\frac{1}{2}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)^{2}\,,$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_q start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ( bold_italic_λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( bold_italic_λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (1)

where $\bm{\hat{\mu}}$μみゅー end_ARG represents the molecular dipole operator and $\hat{H}_{el}$ is the Hamiltonian for the field-free many-electron system. The second term is a harmonic potential introduced by the photon displacement field, with the classic photon displacement coordinate $q_{c}$ and $\omega_{c}$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT being the frequency of the cavity mode. The third term describes the dipole coupling between the molecular system and the photon displacement field, which is characterized by the coupling strength $\bm{\lambda}_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT. The last term is the DSE operator [35, 30, 31], which describes the self-polarization of the molecule-cavity system. The cavity mode-specific coupling parameter $\bm{\lambda}_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT for a cavity with effective mode volume $V_{c}$ is defined as follows:

$$\bm{\lambda}_{c}=\bm{e}_{c}\lambda_{c}=\bm{e}_{c}\sqrt{\frac{4\pi}{V_{c}}}\,.$$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = bold_italic_e start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = bold_italic_e start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT square-root start_ARG divide start_ARG 4 italic_πぱい end_ARG start_ARG italic_V start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG end_ARG . (2)

The unit vector $\bm{e}_{c}$ denotes the polarization axis of the cavity mode. In this work, we make use of two variations of the electronic CBOA Hamiltonian $\hat{H}_{CBO}$, describing the many-electron problem, and solving both cases in the SCF CBO-HF ansatz [28, 32]. In the first case labeled linear CBO-HF we neglected the DSE operator and the resulting energy $E_{CBO}^{lin}$ reads:

$$\bigl{\langle}E_{CBO}^{lin}\bigr{\rangle}=E_{el}+E_{lin}+E_{dis}\quad\text{ with}\quad E_{dis}=\frac{1}{2}\omega_{c}^{2}q_{c}^{2}\,.$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (3)

For the second case, the standard CBO-HF Hamiltonian including the DSE operator is used to self-consistently determine the energy $E_{CBO}$:

$$\bigl{\langle}E_{CBO}\bigr{\rangle}=E_{el}+E_{lin}+E_{dse}+E_{dis}\,.$$ (4)

In addition to considering the full DSE operator, the SCF treatment itself has proved crucial to capture relevant aspects in the description of strongly coupled molecules [29, 28, 32]. To study this aspect, we used the field-free Hartree-Fock energy and expectation values of the dipole moment and the DSE operator to construct a extended Tavis-Cummings (ETC) model [36, 37, 38, 39, 40]. The corresponding energy expectation value that defines the cPESs has the following form

$$\bigl{\langle}E_{ETC}\bigr{\rangle}=E_{el}-\omega_{c}q_{c}\bigl{\langle}\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\bigr{\rangle}_{0}+\frac{1}{2}\bigl{\langle}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)^{2}\bigr{\rangle}_{0}+\frac{1}{2}\omega_{c}^{2}q_{c}^{2}$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_q start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⟨ bold_italic_λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ⟩ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ⟨ ( bold_italic_λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ⟩ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (5)

Both the linear light-matter interaction term, as well as a quadratic DSE term are calculated with the corresponding field-free expectation values $\bigl{\langle}\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\bigr{\rangle}_{0}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ⟩ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and $\bigl{\langle}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)^{2}\bigr{\rangle}_{0}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⋅ overbold_^ start_ARG bold_italic_μみゅー end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ⟩ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT respectively.

II.2 Double Quantum Coherence Spectroscopy

The DQC technique is a 2D-IR method performed with four temporally well-separated laser pulses. The first three pulses with wavevectors $\bm{k}_{1}$, $\bm{k}_{2}$, and $\bm{k}_{3}$ generate a non-linear polarization in the molecular ensemble, which is heterodyne detected with the fourth pulse. We focus on the signal generated along the phase matching direction $\bm{k}_{III}=\bm{k}_{1}+\bm{k}_{2}-\bm{k}_{3}$ [23, 26, 27]. The signal recorded against the three delay times between the pulses $t_{1}$, $t_{2}$, and $t_{3}$ is denoted by $\mathcal{S}(t_{3},t_{2},t_{1})$. By invoking the rotating wave approximation, only the dominant contributions to $\mathcal{S}$, where all interactions are resonant, are retained. Consequently, there are only two contributions to the $\bm{k}_{III}$ signal for our coupled cavity molecule systems represented by the double-sided Feynman diagrams shown in Fig. 1. The two diagrams represent the same evolution during the first two time intervals: during $t_{1}$ the density matrix oscillates with the frequency ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{eg}={ϵ}_e-{ϵ}_g$Ωおめが𝑒𝑔subscriptitalic-ϵ𝑒subscriptitalic-ϵ𝑔\Omega_{eg}=\epsilon_{e}-\epsilon_{g}roman_Ωおめが start_POSTSUBSCRIPT italic_e italic_g end_POSTSUBSCRIPT = italic_ϵ start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT - italic_ϵ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT, and during $t_{2}$ the density matrix oscillates with the frequency ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{fg}={ϵ}_f-{ϵ}_g$Ωおめが𝑓𝑔subscriptitalic-ϵ𝑓subscriptitalic-ϵ𝑔\Omega_{fg}=\epsilon_{f}-\epsilon_{g}roman_Ωおめが start_POSTSUBSCRIPT italic_f italic_g end_POSTSUBSCRIPT = italic_ϵ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_ϵ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. Here $\epsilon_{g}$, $\epsilon_{e}$ and $\epsilon_{f}$ are the vibrational/polaritonic eigenenergies of the vibrational/polaritonic ground state $g$ (green), the first-excitation manifold $e$ (yellow) and the second-excitation manifold $f$ (red). During $t_{3}$ the diagrams yield an oscillation frequency of either ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{e^{\prime }g}$Ωおめがsuperscript𝑒′𝑔\Omega_{e^{\prime}g}roman_Ωおめが start_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_g end_POSTSUBSCRIPT, or ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{f{e}^{\prime }}$Ωおめが𝑓superscript𝑒′\Omega_{fe^{\prime}}roman_Ωおめが start_POSTSUBSCRIPT italic_f italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT (see Figs. 1 a) and  1 b) respectively). Thus, transitions that form a harmonic ladder, that is, ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{e^{\prime }g}={\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{f{e}^{\prime }}$Ωおめがsuperscript𝑒′𝑔subscriptΩおめが𝑓superscript𝑒′\Omega_{e^{\prime}g}=\Omega_{fe^{\prime}}roman_Ωおめが start_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_g end_POSTSUBSCRIPT = roman_Ωおめが start_POSTSUBSCRIPT italic_f italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT, cancel, and such transitions vanish in the DQC signal. Consequently, one of the main features of the DQC technique is its sensitivity to anharmonicities in the molecular systems studied or, in our case, polaritonic systems. Another one is the possibility of having direct access to information on not only the first but also on the second excitation manifold. The obtained 3D signal $\mathcal{S}(t_{3},t_{2},t_{1})$ can be written as double Fourier transform with respect to $t_{3}$ and $t_{2}$:

$$𝒮({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3,{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2,t_1)={\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}𝑑t_3𝑑t_2{e}^{i\left({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3{t}_3+{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2{t}_2\right)}𝒮(t_3,t_2,t_1).$$Ωおめが3subscriptΩおめが2subscript𝑡1superscriptsubscript0superscriptsubscript0differential-dsubscript𝑡3differential-dsubscript𝑡2superscript𝑒𝑖subscriptΩおめが3subscript𝑡3subscriptΩおめが2subscript𝑡2𝒮subscript𝑡3subscript𝑡2subscript𝑡1\mathcal{S}\left(\Omega_{3},\Omega_{2},t_{1}\right)=\int_{0}^{\infty}\int_{0}^% {\infty}dt_{3}dt_{2}e^{i\left(\Omega_{3}t_{3}+\Omega_{2}t_{2}\right)}\mathcal{% S}\left(t_{3},t_{2},t_{1}\right)\,.caligraphic_S ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT , roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_d italic_t start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT italic_d italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_i ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT + roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) end_POSTSUPERSCRIPT caligraphic_S ( italic_t start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) . (6)

Expanding the expression $𝒮({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3,{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2,t_1)$Ωおめが3subscriptΩおめが2subscript𝑡1\mathcal{S}\left(\Omega_{3},\Omega_{2},t_{1}\right)caligraphic_S ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT , roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) in the basis of the polaritonic eigenstates, for $t_{1}=0$, the signal becomes

$$𝒮\left({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3,{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2,t_1=0\right)=\sum _{e,e^{\prime },f}\frac{1}{\left({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2-{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{fg}+i\mathrm{<ruby><rb>\gamma </rb><rt>がんま</rt></ruby>}\right)}\times \left[\frac{{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{g{e}^{\prime }}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{e^{\prime }f}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{fe}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{eg}}{\left({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3-{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{e^{\prime }g}+i\mathrm{<ruby><rb>\gamma </rb><rt>がんま</rt></ruby>}\right)}-\frac{{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{g{e}^{\prime }}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{e^{\prime }f}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{fe}{\mathrm{<ruby><rb>\mu </rb><rt>みゅー</rt></ruby>}}_{eg}}{\left({\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3-{\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_{f{e}^{\prime }}+i\mathrm{<ruby><rb>\gamma </rb><rt>がんま</rt></ruby>}\right)}\right],$$Ωおめが3subscriptΩおめが2subscript𝑡10subscript𝑒superscript𝑒′𝑓1subscriptΩおめが2subscriptΩおめが𝑓𝑔𝑖𝛾delimited-[]subscript𝜇𝑔superscript𝑒′subscript𝜇superscript𝑒′𝑓subscript𝜇𝑓𝑒subscript𝜇𝑒𝑔subscriptΩおめが3subscriptΩおめがsuperscript𝑒′𝑔𝑖𝛾subscript𝜇𝑔superscript𝑒′subscript𝜇superscript𝑒′𝑓subscript𝜇𝑓𝑒subscript𝜇𝑒𝑔subscriptΩおめが3subscriptΩおめが𝑓superscript𝑒′𝑖𝛾\mathcal{S}\left(\Omega_{3},\Omega_{2},t_{1}=0\right)=\sum_{e,e^{\prime},f}% \frac{1}{\left(\Omega_{2}-\Omega_{fg}+i\gamma\right)}\times\left[\frac{\mu_{ge% ^{\prime}}\mu_{e^{\prime}f}\mu_{fe}\mu_{eg}}{\left(\Omega_{3}-\Omega_{e^{% \prime}g}+i\gamma\right)}-\frac{\mu_{ge^{\prime}}\mu_{e^{\prime}f}\mu_{fe}\mu_% {eg}}{\left(\Omega_{3}-\Omega_{fe^{\prime}}+i\gamma\right)}\right]\,,caligraphic_S ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT , roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0 ) = ∑ start_POSTSUBSCRIPT italic_e , italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_f end_POSTSUBSCRIPT divide start_ARG 1 end_ARG start_ARG ( roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - roman_Ωおめが start_POSTSUBSCRIPT italic_f italic_g end_POSTSUBSCRIPT + italic_i italic_γがんま ) end_ARG × [ divide start_ARG italic_μみゅー start_POSTSUBSCRIPT italic_g italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_f end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_f italic_e end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_e italic_g end_POSTSUBSCRIPT end_ARG start_ARG ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT - roman_Ωおめが start_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_g end_POSTSUBSCRIPT + italic_i italic_γがんま ) end_ARG - divide start_ARG italic_μみゅー start_POSTSUBSCRIPT italic_g italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_f end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_f italic_e end_POSTSUBSCRIPT italic_μみゅー start_POSTSUBSCRIPT italic_e italic_g end_POSTSUBSCRIPT end_ARG start_ARG ( roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT - roman_Ωおめが start_POSTSUBSCRIPT italic_f italic_e start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT + italic_i italic_γがんま ) end_ARG ] , (7)

where $\mu_{ij}$μみゅー start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT are the transition dipole moments between polaritonic states $i\rightarrow j$ and $\gamma$γがんま is an empirical dephasing rate.

Figure 1: The two double sided Feynman diagrams that contribute to the DQC signal $\mathcal{S}(t_{3},t_{2},t_{1})$ in the phase-matching direction ($\bm{k}_{III}=\bm{k}_{1}+\bm{k}_{2}-\bm{k}_{3}$) in the rotating wave approximation. The ground state is colored green, the single excited states yellow, and the double excited states are red.

II.3 Computational Details

The vibrational eigenenergies and the transition dipole moments needed to calculate the DQC spectra according to Eq. 7 are obtained using the CBO-HF ansatz, in the Psi4NumPy environment [41], which is an extension of the PSI4 [42] electronic structure package. All calculations were performed using the aug-cc-pVDZ basis set [43] and the geometry of the isolated single HF molecule was optimized at the Hartree-Fock level of theory. In this study, we will consider the interaction between a single lossless mode cavity and an ensemble of $N_{mol}=1$ and $N_{mol}=2$ HF molecules. We assume that the molecules are aligned with the polarization of the cavity mode and we consider a uniform electric field within the cavity. In order to compare the cases $N_{mol}=1$ and $N_{mol}=2$, we apply a scaling factor $1/\sqrt{N_{mol}}$ on the collective coupling strength $\bm{\lambda}_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT to obtain a fixed Rabi splitting.

$$\bm{\lambda}_{c}=\frac{\lambda_{0}}{\sqrt{N_{mol}}}\bm{e}_{c}$$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = divide start_ARG italic_λらむだ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG square-root start_ARG italic_N start_POSTSUBSCRIPT italic_m italic_o italic_l end_POSTSUBSCRIPT end_ARG end_ARG bold_italic_e start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT (8)

Here $\lambda_{0}$λらむだ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is equivalent to $\lambda_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT in Eq. (2) in the single-molecule case. As a result, we increase the mode volume $V_{c}$ of the cavity, but by including more molecules, we keep the average density of molecules $N_{mol}/V_{c}$ fixed. We use a coupling strength $\lambda_{0}$λらむだ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT of $0.03\text{\,}\mathrm{\text{au}}$, which corresponds to a cavity electric field strength of $1.5\text{\,}\mathrm{V}\text{\,}{\mathrm{nm}}^{-1}$ in a Fabry-Pérot-like setup.

By scanning along the bond length of each HF molecule $r_{HF}$ and the photon displacement coordinate $q_{c}$, we construct the $(N_{mol}+1)$-dimensional cPES together with the corresponding dipole moment surfaces. The potential energy surfaces for the cases $N_{mol}=1$ and $N_{mol}=2$ are interpolated to an equally spaced grid of 128$\times$64 ($r_{HF}\times q_{c}$) grid points and 128$\times$128$\times$64 ($r_{HF}\times r_{HF}\times q_{c}$) grid points, respectively. A Gaussian-shaped trial function is numerically propagated in imaginary time  [44] (time step $0.1\text{\,}\mathrm{\text{au}}$ and 70000 time steps) on cPES with the Arnoldi propagation scheme [45] to obtain the first ten nuclear-photonic eigenfunctions and to determine the corresponding transition dipole moments. All quantum dynamics simulations are performed with the open source quantum dynamics code QDng [46]. All calculations were performed in a reproducible environment using the Nix package manager together with NixOS-QChem [47] (commit f5dad404) and Nixpkgs (nixpkgs, 22.11, commit 594ef126).

III Results

In the following, we present DQC spectra simulated for a single HF molecule and a pair of HF molecules resonantly coupled to a single-photon mode of an optical cavity. Self-consistent treatment within the CBO-HF ansatz allows the electronic structure of the molecular ensemble to respond to the cavity field. Moreover, the inclusion of DSE terms allows for a proper description of the self-polarization of the coupled cavity-molecular system [31, 28, 29].

III.1 Vibrational polaritons in hydrogen fluoride molecules

Figure 2 shows schematic energy level diagrams for a single HF molecule and two HF molecules without a cavity in a) and c) and resonantly coupled with a single cavity mode ($\omega_{c}=\omega_{1}$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = italic_ωおめが start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT) in b) and d) respectively. Apart from the expected energetic differences between the results obtained with the three different energy expectation values, the schematic energy-level diagrams shown are the same in all three cases. An analysis of the discussed polaritonic states for both the single-molecule case and the two-molecule case in terms of uncoupled bare states $|v,n\rangle$ can be found in the Supporting Information section S1. In the notation used for the bare states, $v$ describes the molecular vibrational excitation of the uncoupled system and $n$ is the uncoupled photon number.

Figure 2: Schematic energy level diagrams for a single HF molecule and a pair of HF molecules without a cavity in a) and c) and resonantly coupled with a single cavity mode in b) and d). The cavity frequency $\omega_{c}$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is resonant with $\omega_{1}$ωおめが start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. The ground state is colored green, the single excited state is colored yellow, and the double excited states are colored red. The optically dark states originating from the ground state to the single excited manifold and from the single excited manifold to the double excited manifold are shown in gray. Full circles indicate states of predominantly matter character and half-filled circles indicate states with a mixed matter and photon contribution.

In case of a single uncoupled HF molecule, the energy diagram shown in Fig. 2 a), consists of the vibrational ground state $g$, a single excited state $e$ and a double excited state $f$, which correspond to vibrational quantum numbers $v=0,1,2$ respectively. These states are separated by $\omega_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG and $\omega_{2}=$4108\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 4108 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, respectively. The corresponding molecular anharmonicity $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}={\mathrm{<ruby><rb>\omega </rb><rt>おめが</rt></ruby>}}_1-{\mathrm{<ruby><rb>\omega </rb><rt>おめが</rt></ruby>}}_2$Δでるたsubscript𝜔1subscript𝜔2\Delta=\omega_{1}-\omega_{2}roman_Δでるた = italic_ωおめが start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT for HF is $173\text{\,}{\mathrm{cm}}^{-1}$. When a single-cavity photon mode strongly couples to a single HF molecule aligned with the cavity polarization axis, two single excited states and three double excited states in the system are formed, see Fig. 2 b). The cavity frequency $\omega_{c}$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is chosen to be resonant with the fundamental transition $\omega_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG and the coupling strength $\lambda_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is $0.03\text{\,}\mathrm{\text{au}}$. A complementary analysis for the case $\omega_{c}=\omega_{2}$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = italic_ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, corresponding to the first hot transition, can be found in the Supporting Information section S2. Within the single excitation manifold, this resonant interaction leads to the expected formation of a pair of hybrid polaritonic states $UP^{(1)}$ and $LP^{(1)}$ with a Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT of $60\text{\,}{\mathrm{cm}}^{-1}$ for the chosen coupling strength. The lowest energy state $f$ in the second excitation manifold of the coupled cavity-molecule system is characterized mainly by a double excitation in the molecular part ($v=2$, bare state $|2,0\rangle$). Compared to its cavity-free counterpart, this state is stabilized by approximately $20\text{\,}{\mathrm{cm}}^{-1}$. The remaining two states are on average $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}$Δでるた\Deltaroman_Δでるた higher in energy than $f$ and form a second pair of hybrid polaritonic states $UP^{(2)}$ and $LP^{(2)}$ with a Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(2)}$Ωおめが𝑅2\Omega_{R}^{(2)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT of $90\text{\,}{\mathrm{cm}}^{-1}$. As expected, the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(2)}$Ωおめが𝑅2\Omega_{R}^{(2)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT in the second excitation manifold is increasing by approximately a factor of $\sqrt{2}$ compared to ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT. $UP^{(2)}$ and $LP^{(2)}$ are formed by the bare state $|1,1\rangle$, in which the molecule and the cavity are both singly excited, and the bare state $|0,2\rangle$, in which the cavity is doubly excited. None of the states in this coupled cavity-molecule system is expected to be dark with respect to the ground state and the single excited manifold.

In the case of two HF molecules, in order to describe the bare molecular states, we use the basis $|v_{s},v_{a}\rangle$ formed by the two quantum numbers $v_{s}$ and $v_{a}$. The quantum numbers $v_{s}$ and $v_{a}$ correspond to the vibrational excitation of the symmetric and antisymmetric stretching normal modes respectively. For more details, see Section S1 of the supporting information. The two HF molecules are oriented parallel with respect to the cavity polarization axis and separated by $800\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ to ensure that the non-cavity-induced intermolecular couplings are negligible. Figure 2 c) shows the corresponding energy level diagram. The two single excited states $e$ and $d_{1}$ are energetically degenerate. Both correspond to the first excitation of the symmetric stretching mode ($|10\rangle$) and the antisymmetric stretching mode ($|01\rangle$) of the two HF molecules, respectively. This explains why the $e$ state (symmetric mode) can be excited from the ground state, while the other state is dark because of symmetry. Similarly, the two double excited states $f$ and $d_{2}$ are energetically degenerate and correspond to the second excitation of symmetric stretching mode and the antisymmetric stretching mode $|20\rangle$ and $|02\rangle$. Consequently, the $f$ state is a bright state while $d_{2}$ is dark. The remaining state $f_{2}$ is $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}$Δでるた\Deltaroman_Δでるた higher in energy and is formed by a simultaneous single excitation of both the symmetric and antisymmetric stretching mode $|11\rangle$.

By resonantly coupling the two HF molecules to the single cavity mode, again a pair of hybrid polaritonic states $UP^{(1)}$ and $LP^{(1)}$ with a Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT of $60\text{\,}{\mathrm{cm}}^{-1}$ is formed (see Fig. 2 d)). Since the $d_{1}$ state (antisymmetric stretching mode $|01\rangle$) is dark, it is not affected by the coupling to the cavity. The double excitation manifold is more complex, with a total of six states. Both $f$ ($|20\rangle$) and $d_{2}$ ($|02\rangle$) are only slightly or not affected by the cavity interaction and do not change their predominantly molecular character. The four remaining states, $UP^{(2)}$, $MP^{(2)}$, $d_{3}$, and $LP^{(2)}$, all have a photonic contribution, see TABLE S2 in Section S1 of the Supporting Information. However, only $UP^{(2)}$, $MP^{(2)}$, and $LP^{(2)}$ are true hybrid polaritonic states, while the state $d_{3}$ ($|01,1\rangle$) is formed by adding a cavity excitation to the dark state $d_{1}$ ($|01,0\rangle$) without further hybridization. The splitting between the three polaritonic states is $53\text{\,}{\mathrm{cm}}^{-1}$, corresponding to a total Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(2)}$Ωおめが𝑅2\Omega_{R}^{(2)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT of $106\text{\,}{\mathrm{cm}}^{-1}$. The two states $UP^{(2)}$ and $LP^{(2)}$ are mainly characterized by a hybridization of the two bare states $|10,1\rangle$ and $|00,2\rangle$. The middle polariton $MP^{(2)}$ is mostly a linear combination of $|00,2\rangle$ and $|11,0\rangle$, which corresponds to $f_{2}$ in the case of two uncoupled molecules Fig. 2 c).

III.2 DQC spectra for a single HF molecule

Figs. 3 a) and  3 b) show the absolute value of the DQC spectra for a single HF molecule without cavity and resonantly coupled with a single cavity mode, respectively. The normalized absolute values are used in the following to simplify the comparison between spectra obtained with different cPES and number of molecules.

Figure 3: Absolute value of the normalized DQC spectra of a single HF molecule a) without a cavity and b) coupled to a single cavity mode with $\omega_{c}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG. The coupling strength $\lambda_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is $0.03\text{\,}\mathrm{\text{au}}$ and the dephasing $\gamma$γがんま is $10\text{\,}{\mathrm{cm}}^{-1}$. The black horizontal dashed lines mark the energy of the final states, and all signals are labeled $e$ and $g$, indicating that the initial state is the ground state or an intermediate state. The red lines with arrows highlight relevant energy differences.

As illustrated in Fig. 2 a), the HF molecule without cavity is a simple three-level system and the resulting DQC spectra exhibit two peaks (see Fig. 3 a)). Both peaks are resonant at $8393\text{\,}{\mathrm{cm}}^{-1}$ on the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT axis and separated by the molecular anharmonicity $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}$Δでるた\Deltaroman_Δでるた of $173\text{\,}{\mathrm{cm}}^{-1}$ along the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT axis. The one at ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3=4281{\mathrm{cm}}^{-1}$Ωおめが3times4281centimeter1\Omega_{3}=$4281\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG is due to the $g\rightarrow e$ transition, and the signal at $4108\text{\,}{\mathrm{cm}}^{-1}$ corresponds to the $e\rightarrow f$ transition. For the coupled single-molecule single-cavity mode system, we observe three distinct resonances (horizontal dashed lines) on the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT axis in the DQC spectra (see Fig. 3 b)). At ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8369{\mathrm{cm}}^{-1}$Ωおめが2times8369centimeter1\Omega_{2}=$8369\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8369 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, corresponding to the final state $f$, we observe four peaks, which are also the most intense of the whole DQC spectra. These peaks form a pair of doublets, where the two doublets are separated by the molecular anharmonicity $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}$Δでるた\Deltaroman_Δでるた of $173\text{\,}{\mathrm{cm}}^{-1}$. The two peaks of each doublet are separated by the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT of $60\text{\,}{\mathrm{cm}}^{-1}$. The doublet around ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3=4290{\mathrm{cm}}^{-1}$Ωおめが3times4290centimeter1\Omega_{3}=$4290\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = start_ARG 4290 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG is due to the transition from the ground state to the pair of hybrid polaritonic states $UP^{(1)}$ and $LP^{(1)}$, and the other is due to the transition from these hybrid states to the $f$ state. Comparing the DQC spectra shown in Fig. 3 a) and b) the stabilization of the state $f$ due to the cavity interaction is clearly visible. The other two resonances on the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT axis associated with the two hybrid states $UP^{(2)}$ and $LP^{(2)}$ are weaker and separated by a Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(2)}$Ωおめが𝑅2\Omega_{R}^{(2)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT of $90\text{\,}{\mathrm{cm}}^{-1}$. At ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8507{\mathrm{cm}}^{-1}$Ωおめが2times8507centimeter1\Omega_{2}=$8507\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8507 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, final state $LP^{(2)}$, four peaks are distinguishable at $4206\text{\,}{\mathrm{cm}}^{-1}$ ($LP^{(1)}\rightarrow LP^{(2)}$), at $4245\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow LP^{(1)}$), at $4260\text{\,}{\mathrm{cm}}^{-1}$ ($UP^{(1)}\rightarrow LP^{(2)}$) and at $4303\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow UP^{(1)}$) for the chosen dephasing. By careful analysis of the underlying transitions, these signals can be grouped into two doublets, with each peak approximately separated by the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT of $60\text{\,}{\mathrm{cm}}^{-1}$. For the final state $UP^{(2)}$, ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8595{\mathrm{cm}}^{-1}$Ωおめが2times8595centimeter1\Omega_{2}=$8595\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8595 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, only three very weak signals can be observed at $4242\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow LP^{(1)}$), at $4305\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow UP^{(1)}$) and at $4348\text{\,}{\mathrm{cm}}^{-1}$ ($LP^{(1)}\rightarrow UP^{(2)}$). The intensity of the fourth possible transition ($UP^{(1)}\rightarrow UP^{(2)}$) is too weak to be seen in the spectra. The first two peaks are again approximately split by the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT.

Comparisons between DQC signals based on the full CBO-HF surface, the linear CBO-HF surface, and the ETC surfaces are shown as difference spectra $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮$Δでるた𝒮\Delta\mathcal{S}roman_Δでるた caligraphic_S in Fig. 4. $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮$Δでるた𝒮\Delta\mathcal{S}roman_Δでるた caligraphic_S is calculated according to Eq. 9

$$\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮=\left|{𝒮}_{CBO}\right|-\left|{𝒮}_{CBO}^{lin}\right|\text{ and}\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮=\left|{𝒮}_{CBO}\right|-\left|{𝒮}_{ETC}\right|$$Δでるた𝒮subscript𝒮𝐶𝐵𝑂superscriptsubscript𝒮𝐶𝐵𝑂𝑙𝑖𝑛 andΔでるた𝒮subscript𝒮𝐶𝐵𝑂subscript𝒮𝐸𝑇𝐶\Delta\mathcal{S}=\left|\mathcal{S}_{CBO}\right|-\left|\mathcal{S}_{CBO}^{lin}% \right|\quad\text{ and}\quad\Delta\mathcal{S}=\left|\mathcal{S}_{CBO}\right|-% \left|\mathcal{S}_{ETC}\right|roman_Δでるた caligraphic_S = | caligraphic_S start_POSTSUBSCRIPT italic_C italic_B italic_O end_POSTSUBSCRIPT | - | caligraphic_S start_POSTSUBSCRIPT italic_C italic_B italic_O end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_l italic_i italic_n end_POSTSUPERSCRIPT | and roman_Δでるた caligraphic_S = | caligraphic_S start_POSTSUBSCRIPT italic_C italic_B italic_O end_POSTSUBSCRIPT | - | caligraphic_S start_POSTSUBSCRIPT italic_E italic_T italic_C end_POSTSUBSCRIPT | (9)

Figure 4: Difference of the DQC signal $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮$Δでるた𝒮\Delta\mathcal{S}roman_Δでるた caligraphic_S of a single HF molecule coupled to a photon mode with $\omega_{c}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG between a) full CBO-HF and linear CBO-HF and b) full CBO-HF and ETC. Individual DQC spectra are normalized, and the absolute value is used to calculate the difference. The coupling strength $\lambda_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is $0.03\text{\,}\mathrm{\text{au}}$ for both frequencies and the dephasing $\gamma$γがんま is $10\text{\,}{\mathrm{cm}}^{-1}$. The energies of the final states are marked with black dashed lines for the full CBO-HF case and with dashed dotted lines for the linear CBO-HF and ETC cases, respectively.

In both difference spectra the peaks corresponding to the final state $f$ are the most affected. All peaks associated with the $f$ state are red-shifted by approximately 2030 in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and approximately 1020 in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT when the DSE contribution is included, see Fig. 4 a), or the SCF procedure is performed, see Fig. 4 b)). As shown in Fig. 4 a) the $UP^{(2)}$ resonance and the $LP^{(2)}$ resonance and the corresponding peaks are almost unaffected when the DSE contribution is not included. In contrast, when the ETC surface is used to determine the DQC spectra, the resonances of $UP^{(2)}$ and the resonances of $LP^{(2)}$ are affected, see Fig. 4 b). Both are blue shifted by approximately $20\text{\,}{\mathrm{cm}}^{-1}$ in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT for the case without SCF. Also, the intensity distribution changes, $UP^{(2)}$ loses intensity while $LP^{(2)}$ gains it when going from the ETC spectra to the full CBO-HF results.

For completeness, the real ($Re$) and imaginary ($Im$) parts of the DQC spectra of the coupled single-molecule single-cavity mode system are shown in the supporting information section S3, Fig. S6.

III.3 DQC spectra for a pair of HF molecules

The absolute value of the normalized DQC spectra for two identical HF molecules are shown in Fig. 5 for the case without a cavity in a) and resonantly coupled with a single cavity mode in b).

Figure 5: Absolute value of the normalized DQC spectra of two parallel oriented HF molecules a) without a cavity and b) coupled to a single cavity mode with $\omega_{c}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG. The coupling strength $\lambda_{0}$λらむだ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is $0.03\text{\,}\mathrm{\text{au}}$ and the dephasing $\gamma$γがんま is $10\text{\,}{\mathrm{cm}}^{-1}$. The black horizontal dashed lines mark the energy of the final states, and all signals are labeled $e$ and $g$, indicating that the initial state is the ground state or an intermediate state. The red lines with arrows highlight relevant energy differences.

The DQC spectra of two parallel oriented HF molecules not interacting with a cavity field is shown in Fig. 5 a) and is nearly identical to the single-molecule situation, see Fig. 3 a) for comparison. The bright state $f_{2}$, which is described by a simultaneous single excitation of both the symmetric and antisymmetric stretching mode, is not visible because of an exact cancellation of the two Liouville paths. The corresponding energy splittings between the $g\rightarrow e$ transition and the $e\rightarrow f_{2}$ transition are identical, leading to a harmonic ladder that is not visible in a DQC spectrum. A black dotted line in Fig. 5 a) indicates the energetic position of the corresponding hypothetical $f_{2}$ resonances. The two additional vibrational states ($d_{1}$ and $d_{2}$) present in the case of two molecules are dark, as illustrated in Fig. 2 c), since they are associated with excitations of the antisymmetric linear combination of the two stretching modes. As a consequence, the resulting DQC spectrum has only two peaks and, thus representing an effective three-level system. However, once the two HF molecules are resonantly coupled to the cavity mode, the resulting DQC spectrum, shown in Fig. 5 b), is clearly distinguishable from the single molecule case. We observe four distinct resonances (horizontal dashed lines) on the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT axis in the DQC spectra, since the other two possible double excited states $d_{2}$ and $d_{3}$ are dark for symmetry reasons. Similarly to the case of a single HF molecule, we observe four peaks at ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8374{\mathrm{cm}}^{-1}$Ωおめが2times8374centimeter1\Omega_{2}=$8374\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8374 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, corresponding to the final state $f$, which are also the most intense of the whole DQC spectrum. Interestingly, the resonance of the $f$ state is nearly identical in both the one-molecule and two-molecule cases. The four associated peaks can again be grouped as a pair of doublets separated by the molecular anharmonicity $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}$Δでるた\Deltaroman_Δでるた of about $173\text{\,}{\mathrm{cm}}^{-1}$ and the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT of $65\text{\,}{\mathrm{cm}}^{-1}$. Due to the rescaling of the coupling strength according to Eq. 8 ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT is similar to the one-molecule case. The doublet around ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3=4274{\mathrm{cm}}^{-1}$Ωおめが3times4274centimeter1\Omega_{3}=$4274\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = start_ARG 4274 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG is due to the transition from the ground state to the pair of hybrid polaritonic states $UP^{(1)}$ and $LP^{(1)}$, and the other one is due to the transition from these hybrid states to the $f$ state. The three hybrid polaritonic states $UP^{(2)}$, $MP^{(2)}$, and $LP^{(2)}$ have visible resonances on the ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT axis. These three resonances are separated by about $53\text{\,}{\mathrm{cm}}^{-1}$, corresponding to a total Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(2)}$Ωおめが𝑅2\Omega_{R}^{(2)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT of $106\text{\,}{\mathrm{cm}}^{-1}$, consistent with the observed increase of the Rabi splitting in the single molecule case. At ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8496{\mathrm{cm}}^{-1}$Ωおめが2times8496centimeter1\Omega_{2}=$8496\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8496 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, final state $LP^{(2)}$, we see an intense peak at $4248\text{\,}{\mathrm{cm}}^{-1}$ comparable to the peaks corresponding to the final state $f$ and two weaker signals at $4197\text{\,}{\mathrm{cm}}^{-1}$ ($LP^{(1)}\rightarrow LP^{(2)}$) and at $4303\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow UP^{(1)}$). These three signals are approximately separated by the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT. The high intensity of the central peak can be explained by the overlap of the $g\rightarrow LP^{(1)}$ and $UP^{(1)}\rightarrow LP^{(2)}$ transitions. For the final state $MP^{(2)}$, ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8551{\mathrm{cm}}^{-1}$Ωおめが2times8551centimeter1\Omega_{2}=$8551\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8551 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, only two signals can be observed at $4244\text{\,}{\mathrm{cm}}^{-1}$ and at $4306\text{\,}{\mathrm{cm}}^{-1}$ separated again by the Rabi splitting ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_R^{(1)}$Ωおめが𝑅1\Omega_{R}^{(1)}roman_Ωおめが start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT. These two signals are formed by the transitions $g\rightarrow LP^{(1)}$ and $LP^{(1)}\rightarrow MP^{(2)}$ and the transitions $UP^{(1)}\rightarrow MP^{(2)}$ and $g\rightarrow UP^{(1)}$ respectively. As mentioned in section III.1 the final state $MP^{(2)}$ is formed by hybridization of the two bare states $|00,2\rangle$ and $|11,0\rangle$. Without coupling to a cavity, this $|11,0\rangle$ state, labeled $f_{2}$ in Fig. 5 a), is not visible because of a cancellation of the Liouville diagrams. In the cavity, however, this state becomes visible due to the hybridization, which induces anharmonicity in the corresponding excitation ladder. The identical process is observed when the cavity is resonant with the first hot transition and discussed in Section S2 of the supporting information. Similarly to the single-molecule case, only two very weak signals corresponding to the final state $UP^{(2)}$, ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2=8610{\mathrm{cm}}^{-1}$Ωおめが2times8610centimeter1\Omega_{2}=$8610\text{\,}{\mathrm{cm}}^{-1}$roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = start_ARG 8610 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, can be observed at $4242\text{\,}{\mathrm{cm}}^{-1}$ ($g\rightarrow LP^{(1)}$) and at $4359\text{\,}{\mathrm{cm}}^{-1}$ ($UP^{(1)}\rightarrow UP^{(2)}$). All other possible transitions are too weak to be observed in the spectra. There is no clear indication that any of the three dark states $d_{1}$ and $d_{2}$ takes part in the signal and thus these states are not needed to explain the DQC signals. As for the single molecule case, for completeness we show the real ($Re$) and imaginary ($Im$) parts of the DQC spectra of the two-molecule case coupled to a single-cavity mode in the supporting information section S3 Fig. S7.

In Fig. 6 the comparison between the DQC signals based on the full CBO-HF surface, the linear CBO-HF surface, and the ETC surfaces for the two-molecule case is shown as difference spectra $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮$Δでるた𝒮\Delta\mathcal{S}roman_Δでるた caligraphic_S calculated according to Eq. 9.

Figure 6: Difference $\mathrm{<ruby><rb>\Delta </rb><rt>でるた</rt></ruby>}𝒮$Δでるた𝒮\Delta\mathcal{S}roman_Δでるた caligraphic_S of the DQC signal of two parallel HF molecule coupled to a photon mode with $\omega_{c}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG between a) full CBO-HF and linear CBO-HF and b) full CBO-HF and ETC. Individual DQC spectra are normalized, and the absolute value is used to calculate the difference. The coupling strength $\lambda_{c}$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is $0.03\text{\,}\mathrm{\text{au}}$ for both frequencies and the dephasing $\gamma$γがんま is $10\text{\,}{\mathrm{cm}}^{-1}$. The energies of the final states are marked with black dashed lines for the full CBO-HF case and with dashed dotted lines for the linear CBO-HF and ETC cases, respectively.

In agreement with the single-molecule case, the peaks corresponding to the state $f$ are most affected in both difference spectra. For the comparison between the full CBO-HF and linear CBO-HF results shown in Fig. 6 a), the $f$ peaks are only red-shifted by about $20\text{\,}{\mathrm{cm}}^{-1}$ in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, while no general trend for shifts in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT is apparent. When the DSE term is not included (see Fig. 6 a)), the resonances of the two hybrid polaritonic states $MP^{(2)}$ and $LP^{(2)}$ are affected. In particular, for the $LP^{(2)}$ state, a red shift is observed in both ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT for the strong signal around $4248\text{\,}{\mathrm{cm}}^{-1}$ without DSE included. Interestingly, the weak $UP^{(2)}$ remains almost unchanged, similar to the case of a single molecule. In contrast, when the ETC surface is used to determine the DQC spectra, the observed changes in the difference spectra in the single-molecule case and the two-molecule case are quite similar. For the $f$ signals, almost the same red shift is observed in both ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT, compare Fig. 4 b) and Fig. 6 b). The resonances of three hybrid polaritonic states $UP^{(2)}$, $MP^{(2)}$, and $LP^{(2)}$ are all blue shifted by roughly $20\text{\,}{\mathrm{cm}}^{-1}$ in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT for the ETC case. Consistent with the single-molecule results, the intensity distribution changes, $UP^{(2)}$ loses intensity while $LP^{(2)}$ gains it when going from the ETC spectra to the full CBO-HF results. In general, the observed effect on the DQC spectra of the SCF treatment is the same for one and two HF molecules, where the DSE included in the CBO-HF approach shows additional changes in the spectra when going from one to two molecules.

As a final step, we want to further analyze the differences in the DQC spectra when going from one to two HF molecules for the full CBO-HF approach and using the ETC model. The corresponding difference spectra are shown in Fig. 7

Figure 7: Difference of the absolute values of the DQC spectra of two parallel oriented HF molecules and a single HF molecule coupled to a photon mode with $\omega_{0}=$4281\text{\,}{\mathrm{cm}}^{-1}$$ωおめが start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = start_ARG 4281 end_ARG start_ARG times end_ARG start_ARG power start_ARG roman_cm end_ARG start_ARG - 1 end_ARG end_ARG, $\lambda_{c}=$0.03\text{\,}\mathrm{\text{au}}$$λらむだ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = start_ARG 0.03 end_ARG start_ARG times end_ARG start_ARG au end_ARG and the dephasing $\gamma$γがんま is $10\text{\,}{\mathrm{cm}}^{-1}$. a) The full CBO-HF cPESs are used to construct the DQC spectra, and b) the ETC model cPESs are used. The energies of the final states are marked with black dashed lines for the two molecule case.

Analyzing the difference spectrum in Fig. 7 a) for the full CBO-HF case, two main changes appear when going from one to two HF molecules. The most striking change is the increase in intensity for the signal around $4248\text{\,}{\mathrm{cm}}^{-1}$, which corresponds to the final state $LP^{(2)}$. As discussed for the two-molecule DQC spectrum shown in Fig. 5 b), this intense peak is due to two overlapping signals that are still separated in the single-molecule DQC spectrum (see Fig. 3 b)). The other change is a blue shift in ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT of all signals corresponding to the $f$ resonance. For the ETC case shown in Fig. 7 b) all peaks associated with the final state $f$ are shifted to lower frequencies in both ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_2$Ωおめが2\Omega_{2}roman_Ωおめが start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and ${\mathrm{<ruby><rb>\Omega </rb><rt>おめが</rt></ruby>}}_3$Ωおめが3\Omega_{3}roman_Ωおめが start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT. With respect to polaritonic sates, only smaller changes are observed when going from one to two HF molecules. The distinct differences in the comparison between the full CBO-HF results and those obtained using the ETC model are due solely to the SCF process, which allows the electronic structure to respond to the cavity field mode. This response seems to be particularly relevant when going beyond a single-molecule situation because of cavity-induced molecular interactions.

IV Summary and Conclusion

Based on the recently formulated cavity Born-Oppenheimer Hartree-Fock ansatz [28, 32] and inspired by the work of Saurabh and Mukamel [27] we simulated two-dimensional DQC spectra for the rather anharmonic case of one and two diatomic hydrogen fluoride HF molecules coupled to an infrared cavity. In both cases, the molecular system is coupled to a single-photon mode that is resonant with the first vibrational transition in HF. Using ab initio cPES at the CBO-HF level of theory, we could demonstrate how ground-state vibrational excitation manifolds are modified upon coupling to a cavity mode. As a consequence, the molecular anharmonicities, the molecule-cavity interaction, and the self-polarization of the molecule-cavity system are naturally included. Even in the rather simple diatomic case of HF, coupling to an optical cavity results in a rather complex and ”asymmetric” DQC signal compared to model systems reported in the literature [27]. This complexity comes from the non-idealized second excitation manifold, where due to the anharmonicity of HF both nearly pure molecular states and the expected hybrid polaritonic states contribute, as shown in Fig. 2. The significantly different character of these states leads to a mixture of rather weak and rather strong transitions. We were able to show that neglecting the DSE contribution or not performing the SCF procedure has a significant effect on the DQC signals. In both cases, the resonances of the $f$ state are affected, which is mainly an excited molecular final state. Therefore, it allows us to monitor the effect of DSE and SCF on the LP and UP states in the single excitation manifold. However, for the single-molecule case without SCF, the $UP^{(2)}$ state and the $LP^{(2)}$ state also change. For the case of two HF molecules, the DSE contribution not only influences the $f$ resonances but also leads to a strong signal corresponding to the final state $LP^{(2)}$, which is clearly weaker without DSE or SCF. This $LP^{(2)}$ signal is also the striking difference when going from a single molecule to a pair of molecules. In contrast, the resonances of the $f$ state are nearly identical as the number of molecules increases. An important result of our study is that we show that the DQC spectrum is highly sensitive to the effect of the DSE. In particular, when considering the case of two molecules without a cavity, the molecular state $|11\rangle$ corresponding to two simultaneously excited molecules do not contribute to the DQC spectrum due to the fact that the two molecules do not interact. However, when considering two molecules inside the cavity, the cavity induce an effective interaction between the two molecules resulting from the DSE which can be directly probed by the DQC spectrum through the MP⁽²⁾ state Therefore, we show that the DQC technique offer a unique insight into the cavity-mediated intermolecular interactions. DQC spectroscopy can be used to have a better understanding of many-body effects in molecular systems under vibrational strong coupling which in turn may thus provide a deeper mechanistic insight into chemical reactions under vibrational strong coupling.

Supplementary Material

See the supplementary material for an analysis of the vibrational polaritonic eigenstates in terms of bare states, DQC spectra for a cavity resonant with the first hot transition of the HF molecule and the real and imaginary parts of the DQC spectra.

Author Declaration Section

Conflict of Interest Statement

The authors have no conflicts to disclose

Author Contributions

Thomas Schnappinger: Data curation (lead); Formal analysis (lead); Conceptualization (equal); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Cyril Falvo: Conceptualization (equal); Investigation (equal); Writing – review & editing (equal). Markus Kowalewski: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.