# Equivariant graph neural network interatomic potential for
Green-Kubo thermal conductivity in phase change materials

Sung-Ho Lee  Univ. Grenoble Alpes, CEA, Leti, F-38000, Grenoble, France  Jing
Li  Univ. Grenoble Alpes, CEA, Leti, F-38000, Grenoble, France  Valerio
Olevano  Institut Néel, CNRS & Univ. Grenoble Alpes, F-38042, Grenoble, France
Benoit Sklénard<EMAIL_ADDRESS>Univ. Grenoble Alpes, CEA, Leti,
F-38000, Grenoble, France

###### Abstract

Thermal conductivity is a fundamental material property that plays an
essential role in technology, but its accurate evaluation presents a challenge
for theory. In this work, we demonstrate the application of $E(3)$-equivariant
neutral network interatomic potentials within Green-Kubo formalism to
determine the lattice thermal conductivity in amorphous and crystalline
materials. We apply this method to study the thermal conductivity of germanium
telluride (GeTe) as a prototypical phase change material. A single deep
learning interatomic potential is able to describe the phase transitions
between the amorphous, rhombohedral and cubic phases, with critical
temperatures in good agreement with experiments. Furthermore, this approach
accurately captures the pronounced anharmonicity that is present in GeTe,
enabling precise calculations of the thermal conductivity. In contrast, the
Boltzmann transport equation including only three-phonon processes tends to
overestimate the thermal conductivity by approximately a factor of 2 in the
crystalline phases.

##  I Introduction

Thermal conductivity, as an intrinsic material property, holds significant
implications in technology as it plays a crucial role in determining the
thermal management of electronic devices [1, 2] and serves as a key parameter
in thermoelectric device performance [3, 4]. Heat transport in semiconductors
and insulators is primarily governed by lattice vibrations, i.e., phonons.
Substantial efforts have been directed towards precise calculations of lattice
thermal conductivities from a microscopic perspective. The main methods for
calculating lattice thermal conductivity include the Boltzmann transport
equation (BTE) [5, 6, 7], nonequilibrium Green’s function (NEGF) theory [8, 9,
10], nonequilibrium molecular dynamics (NEMD) [11] and the Green-Kubo (GK)
formula [12, 13, 14]. The BTE assesses the response of phonon occupation to a
temperature gradient, typically including three-phonon scattering processes,
which limits its applicability to weakly anharmonic crystalline materials. In
contrast, NEGF treats phonons quantum mechanically, considering interface
scatterings and phonon anharmonicity through self-energies, but it comes with
a computational cost [8]. The NEMD method is based on a direct simulation of
stationary heat fluxes in a nonequilibrium state and naturally include
anharmonic effects. Finally, GK provides the lattice thermal conductivity from
the autocorrelation of heat fluxes usually computed in an equilibrium
molecular dynamics (MD) simulation, accounting for anharmonic effects to all
orders [15]. Furthermore, recent developments extend GK to low temperatures
[13, 14], which makes it a robust approach for a wide range of temperatures
and materials. GK theory provides a unified approach to compute the lattice
thermal conductivity in ordered and disordered solids. For harmonic amorphous
systems, thermal transport can be described by the Allen and Feldman (AF)
theory [16]. However, it has been shown that the AF theory may be inadequate
when anharmonic effects become important [17, 18].

The MD simulation in the GK approach requires a relatively long simulation
time (up to a few nanoseconds) for sufficient statistical sampling and an
accurate description of interactions among atoms. Such long simulation times
are affordable for MD with empirical force fields, but at the price of reduced
accuracy and universality. Ab initio MD has better accuracy but is too
computationally expensive for large systems or long MD simulations.
Extrapolation schemes have been proposed [19] to reduce the computational
cost, but they are unsuitable for disordered solids.

Over the past few years, machine learning has emerged as a practical
alternative for tasks where ab initio methods encountered difficulties.
Notably, machine learning interatomic potentials (MLIP) have proven successful
in rapidly predicting energies, forces, and stress tensors with an accuracy
comparable to first-principle methods. Particularly in thermal transport GK
calculations, MLIPs have been employed using descriptor-based approaches, such
as Behler-Parrinello neural networks or kernel-based methods [20, 21, 22].
Recently, neural network interatomic potentials utilizing message passing
architectures [23, 24, 25, 26] have been proposed as an alternative to hand-
crafted descriptors. In this approach, structures are encoded as graphs, where
atoms are depicted as nodes connected by edges. In initial models, the
information at the nodes and edges was made invariant with respect to the
Euclidean group $E(3)$ (i.e., the group of translations, rotations, and
inversions in Euclidean space), and the atomic representations were limited to
scalar interatomic distances [23]. Such models have since been largely
superseded by architectures built on convolution operations that are
equivariant with respect to the $E(3)$ group. In equivariant approaches,
isometric transformations on the relative atomic displacement vector inputs
are propagated through the network to correspondingly transform the outputs.
Equivariant MLIPs have been shown to achieve substantially improved data
efficiency and unprecedented accuracy compared to their invariant counterparts
[24, 25, 26]. However, in such message passing architectures, many-body
interactions are captured by iteratively propagating information along the
graph at each layer in the network. This has the effect of extending the local
receptive field of an atom to significantly beyond the cutoff radius, which
renders parallelization impractical [27]. Recently, a strictly local
equivariant neural network potential has been proposed to address this
drawback [27]. In this architecture, information is stored as a per-pair
quantity, and instead of nodes exchanging information with its neighbours via
edges, a convolution operation acts on the cutoff sphere in the form of a set
of invariant (scalar) latent features and a set of equivariant (tensor) latent
features that interact at each layer.

In this work, we employ the strictly local $E(3)$-equivariant neural network
potential from Ref. [27] to compute the temperature-dependent thermal
conductivity of germanium telluride (GeTe) in various phases using GK theory.
GeTe is a chalcogenide material employed in many technological applications,
such as phase change nonvolatile memory storage [28, 29], thermoelectricity
[3, 30, 31] and spintronics [32, 33, 34]. It undergoes a ferroelectric phase
transition from the low temperature rhombohedral $\alpha$-GeTe (spacegroup
$R3m$) to the cubic $\beta$-GeTe (spacegroup $Fm\bar{3}m$) at a Curie
temperature of $T_{c}\approx 650-700$ K [35, 36, 37]. Amorphous GeTe also
plays an important role in technological applications such as phase change
memories. Hence, GeTe serves as a prototype material to demonstrate the
universality of the GK method combined with equivariant neural network
interatomic potentials. This approach effectively captures strong
anharmonicity and materials in diverse crystalline phases or amorphous states
for the calculation of lattice thermal conductivity.

##  II Methodology

The thermal conductivity tensor within GK theory is defined as

$\kappa_{\alpha\beta}(T)=\dfrac{1}{k_{\mathrm{B}}T^{2}V}\lim_{\tau\to\infty}\int_{0}^{\tau}dt\,\langle
j_{\alpha}(t)\cdot j_{\beta}(0)\rangle_{T},$ (1)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ the temperature, $V$ the
volume, $j_{\alpha}(t)$ the $\alpha$-th Cartesian component of the macroscopic
heat flux, and $\langle j_{\alpha}(t)\cdot j_{\beta}(0)\rangle_{T}$ the heat
flux autocorrelation function, with the symbol $\langle\cdot\rangle_{T}$
denoting the ensemble average over time and over independent MD trajectories.

The total heat flux of a system of $N$ atoms is given by

$\bm{j}(t)=\sum_{i=1}^{N}\dfrac{d}{dt}\left(\bm{r_{i}}E_{i}\right),$ (2)

where $E_{i}=m_{i}\bm{v}_{i}^{2}/2+U_{i}$ is the total energy (i.e. kinetic
and potential energy) of atom $i$ with mass $m_{i}$, velocity $\bm{v}_{i}$ and
atomic position $\bm{r_{i}}$. In MLIPs, the partitioning $U=\sum_{i}U_{i}$ of
the total energy of the system into atomic contributions $U_{i}$ allows the
total heat flux of a periodic system to be expressed as [12]

$\bm{j}(t)=\sum_{i=1}^{N}\bm{v}_{i}E_{i}-\sum_{i=1}^{N}\sum_{j\neq
i}\bm{r}_{ij}\left(\dfrac{\partial
U_{i}}{\partial\bm{r}_{ij}}\cdot\bm{v}_{j}\right)$ (3)

where the sum over $j$ runs over the atoms that are within the cutoff radius
$r_{c}$ of atom $i$. We implemented the calculation of Eq. (3) in the LAMMPS
code [38]. The term $\partial U_{i}/\partial\bm{r}_{ij}$ is obtained by
automatic differentiation and was also used for the calculation of the virial
tensor [39, 12], which is required to perform simulations in the isothermal-
isobaric (NpT) ensemble.

##  III Machine Learning Interatomic Potential

###  III.1 Training dataset

To create the reference dataset for training the MLIP, ab initio MD
simulations based on density functional theory (DFT) were carried out using
the VASP code [40, 41], with temperatures ranging from 100 K to 2500 K. The
generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) [42]
was used for the exchange-correlation energy and Grimme’s D3 dispersion
correction [43] was applied. The supercells contained 192 and 216 atoms for
the initial rhombohedral and cubic structures, respectively. Then, a total of
6000 structures were extracted from the MD trajectories and recomputed to
obtain more accurate energies, forces, and stress tensors. We used an energy
cutoff of 400 eV and a $2\times 2\times 2$ $k$ mesh to sample the Brillouin
zone. The equivariant NN model was trained using the Allegro package [27].
Details of the training procedure and dataset partitioning are provided in
Appendix A. The training and test GeTe datasets can be found on Zenodo [44].
The root mean squared errors (RMSE) and mean absolute errors (MAE) on the
predicted energies, forces and stress tensors on the test dataset are $0.90$
meV/atom, $29.87$ meV/Å, $0.28$ meV/$\text{\AA}^{3}$ and $1.07$ meV/atom,
$42.97$ meV/Å, $0.37$ meV/$\text{\AA}^{3}$, respectively.

((a))

((b))

Figure 1: Comparison of phonon dispersions computed with DFT and with the MLIP
of LABEL:sub@fig:phonon-rhomb $\alpha$-GeTe and LABEL:sub@fig:phonon-cub
$\beta$-GeTe

###  III.2 Model validation

To validate the MLIP further, the equilibrium geometries of crystalline GeTe
were optimized using the MLIP. For $\alpha$-GeTe, the lattice parameter was
$a=4.42$ Å and the angle $\alpha=57.13^{\circ}$, closely matching the DFT
results ($a=4.41$ Å and $\alpha=57.42^{\circ}$). Similarly, for $\beta$-GeTe,
the MLIP yields $a=4.24$ Å, in excellent agreement with the lattice parameter
from DFT of $a=4.23$ Å.

Moreover, the phonon dispersion obtained from the MLIP is in excellent
agreement with DFT for both $\alpha$ and $\beta$-GeTe, as shown in Fig. 1. In
particular, our model describes optical phonons well, which is usually
challenging for MLIPs [45, 20]. Imaginary soft phonon modes in cubic GeTe are
also well described by the MLIP, which is essential in capturing the phase
transition [46, 47]. Phonon dispersions were computed using the finite
displacement method implemented in Phonopy [48] with $3\times 3\times 3$ and
$5\times 5\times 2$ supercells of the conventional unit cells for cubic and
rhombohedral phases, respectively. For the DFT calculations, we used the same
settings as those used to generate the reference dataset. LO-TO splitting was
not included in our calculations as long-range Coulomb interactions tend to be
screened by free carriers in real samples [49].

((a))

((b))

Figure 2: Evolution of LABEL:sub@fig:lattice-temp-a the lattice parameter $a$
and LABEL:sub@fig:lattice-temp-angle the angle $\alpha$ as a function of
temperature in the NpT MD simulations of crystalline GeTe, compared against
experimental data from Ref. [36, 35, 37]. Simulated lattice parameters in
LABEL:sub@fig:lattice-temp-a were shifted by $-0.1~{}\text{\AA}$.

##  IV Results and discussion

###  IV.1 Finite temperature lattice dynamics

We investigated the lattice dynamics of GeTe through MD simulations across the
$\alpha\rightarrow\beta$ phase transition with our MLIP. For each temperature,
rhombohedral and cubic GeTe supercells were first equilibrated for at least
200 ps in the NpT ensemble at ambient pressure with a 2 fs timestep in order
to obtain the averaged temperature-dependent structural parameters shown in
Fig. 2. The rhombohedral lattice parameter $a$ and angle $\alpha$ reach cubic
values for T>650 K, in good agreement with experimental data. Interestingly,
below the Curie temperature, the cubic phase remains confined within a local
energy minimum and does not undergo transformation into the rhombohedral phase
throughout our MD simulations.

By employing the temperature-dependent effective-potential (TDEP) method [50,
51, 52], the temperature-dependent interatomic force constants (IFCs) were
extracted from configurations sampled from MD trajectories. We used 512-atom
rhombohedral or cubic GeTe supercells using the temperature-dependent
structural parameters from Fig. 2. Each structure was equilibrated in the NVT
ensemble for 80 ps (thermalization step). Then, snapshots were taken from a
600 ps MD simulation in the microcanonical ensemble for each temperature to
extract the IFCs. Details of the TDEP extraction are provided in Appendix B.
Temperature dependent phonon dispersions and spectral functions obtained using
these IFCs are shown in Fig. 3 and 4 for the $\alpha$ and $\beta$ phases,
respectively.

((a))

((b))

Figure 3: Temperature dependent phonon dispersions (white lines) and spectral
functions of $\alpha$-GeTe computed using the TDEP method at
LABEL:sub@fig:tdep-rhomb-300K 300 K and LABEL:sub@fig:tdep-rhomb-600K 600 K
(just before the $R3m\rightarrow Fm\bar{3}m$ phase transition).

((a))

((b))

Figure 4: Temperature dependent phonon dispersions (white lines) and spectral
functions of $\beta$-GeTe computed using the TDEP method at
LABEL:sub@fig:tdep-cub-300K 300 K and LABEL:sub@fig:tdep-cub-700K 700 K.

In Fig. 5, the evolution of the longitudinal and transverse optical phonon
modes ($\Gamma_{6}$ and $\Gamma_{4}$, respectively) is depicted as a function
of temperature. The softening of these two modes up to the Curie temperature
is consistent with previous theoretical studies [47, 46] and is comparable to
experiments [53, 54, 49]. Beyond 650 K, the optical phonons merge, indicating
the transition to the cubic phase where optical phonons exhibit three-fold
degeneracy.

Figure 5: Temperature evolution of A1 and E optical phonon modes computed with
the TDEP method and compared against experimental data from Ref. [53, 54, 49].

###  IV.2 Green-Kubo thermal conductivity

To compute the GK thermal conductivity of both crystalline and amorphous GeTe,
MD simulations with the MLIP were performed at different temperatures. The
amorphous GeTe structure was generated using a melt-quench process (see
Appendix C). Heat flux was computed during MD simulations in the
microcanonical ensemble and the ensemble average was carried out over
independent trajectories of at least 1 ns after equilibration in the NpT
ensemble.

Figure 6(a) shows the averaged heat flux autocorrelation function (HFACF) at
300 K as a function of the correlation time $\tau$ for the initial GeTe
structures in the rhombohedral, cubic, and amorphous phases. The average HFACF
features large oscillations before decaying to zero, mainly due to the kinetic
contribution of the heat flux (first term in Eq. 3). The cumulative time
integration of the HFACF yields the thermal conductivities
$\kappa_{avg}=1/3\sum_{\alpha=1}^{3}\kappa_{\alpha\alpha}$ that are plotted in
Fig. 6(b).

((a))

((b))

Figure 6: LABEL:sub@fig:hfacf-300K Averaged heat flux autocorrelation function
(HFACF), normalized by its value at $\tau=0$ and LABEL:sub@fig:kappa-300K
lattice thermal conductivity computed using GK theory as a function of
correlation time at $T=300$ K for rhombohedral, cubic and amorphous GeTe.
Figure 7: Calculated thermal conductivities using Green-Kubo theory (solid
lines) of amorphous and crystalline GeTe. The thermal conductivity computed
using direct solution of linearized Boltzmann transport equation are also
reported (dotted line for three phonon, dashed line for three and four
phonon). Experimental data are shown for crystalline GeTe from Ref. [55, 30,
56]) (empty symbols) and for amorphous GeTe from Refs. [56] (green symbols).

After testing the convergence with respect to system size (see Appendix D), we
employed supercells containing 360 and 512 atoms for crystalline and amorphous
structures, respectively. Figure 7 shows the computed lattice thermal
conductivities as a function of temperature for both amorphous and crystalline
GeTe together with existing experimental data from Refs. [55, 30, 56].
Comparing against experiments is challenging because experimental values of
lattice thermal conductivities of crystalline GeTe show a large dispersion.
This variability can be attributed to two main reasons. First, thermal
conductivity comprises a lattice contribution and an electronic contribution.
Therefore, experimental lattice thermal conductivity is an indirect
measurement, which is obtained by removing the electronic contribution from
the measured thermal conductivity, typically evaluated using the Wiedmann-
Franz law that introduces an additional approximation from the Lorenz number.
Second, the quality of the samples varies, leading to variations in thermal
conductivity measurements. Despite the significant experimental variations
mentioned above, the calculated GK thermal conductivity values are found to
fall within the range of experimental values.

The GK lattice thermal conductivity for the amorphous phase (solid green line)
is in excellent agreement with the experimental data of Ref. [56] (green
squares). This can be regarded as a direct comparison with the experiment
since the electronic contribution to the thermal conductivity was found to be
negligible in amorphous GeTe [57]. A previous study obtained a similar value
of $0.27\pm 0.05$ W$\cdot$m${}^{-1}\cdot$K-1 at 300 K from GK simulations with
a Behler-Parrinello-type MLIP [22]. The predicted thermal conductivity for
amorphous GeTe is constant until $\sim 450$ K. It then starts to increase,
indicating a transition to a crystalline phase, as evidenced by the evolution
of the radial distribution function (see Fig. 10.) and consistent with the
amorphous-crystalline phase transition temperature observed experimentally
[56].

###  IV.3 BTE thermal conductivity

To compute the thermal conductivity using the BTE, we employed the FourPhonon
code [58] with the TDEP IFCs up to fourth order and an $11\times 11\times 11$
$q$-mesh. This allows a direct comparison between GK and BTE methods, as both
calculations were conducted under the same conditions, utilizing identical
interatomic potentials and temperatures; the only difference lies in the
thermal transport formalism employed. BTE including only three-phonon
processes overestimates the thermal conductivity by about 1.8
W$\cdot$m${}^{-1}\cdot$K-1, which is about twice the GK result at 300 K, and
about three times that at 900 K. Including four-phonon processes significantly
improves BTE results. The BTE thermal conductivity approaches that of GK,
especially at high temperatures in the cubic phase. However, the discrepancy
at 300 K remains about $50\%$. Such overestimation indicates that higher-order
IFCs are mandatory in BTE to capture the strong anharmonicity. Furthermore,
the computational expenses escalate considerably when incorporating four-
phonon processes in BTE. Consequently, the GK method emerges as a viable
approach for assessing the lattice thermal conductivity of materials
exhibiting strong anharmonicity. We observed differences when comparing our
results to previous BTE calculations reported in the literature [47, 56, 59].
These discrepancies are discussed in Appendix E.

##  V Conclusion

In conclusion, we developed an equivariant graph neural network interatomic
potential to study thermal transport in both amorphous and crystalline GeTe.
The potential describes GeTe at a near-ab initio level of accuracy for the
rhombohedral, cubic and amorphous phases with a single model. Furthermore, it
successfully captures phase transitions with transition temperatures in good
agreement with experimental data. When coupled with Green-Kubo theory, it can
determine the lattice thermal conductivity not only for strongly anharmonic
crystals, but also for the amorphous phase.

###### Acknowledgements.

We thank F. Bottin and J. Bouchet for discussions about TDEP calculations.
This work was performed using HPC/AI resources from GENCI–IDRIS (Grant No.
2022-A0110911995) and was partially funded by the European commission through
ECSEL-IA 101007321 project StorAIge and the French IPCEI program.

##  Appendix A Machine Learning Interatomic Potential generation

###  A.1 Training

The equivariant graph neural network (GNN) model was trained using the Allegro
code [27] on an NVIDIA A100 GPU with float32 precision. The dataset contained
6000 reference configurations, randomly split into 5000, 500 and 500
configurations for the training, validation and test sets, respectively [44].
The training was done using a mean squared error (MSE) loss function based on
a weighted sum of energy, forces and stress loss terms:

$\mathcal{L}=\dfrac{\lambda_{E}}{B}\sum_{b}^{B}\left(\dfrac{\hat{E}_{b}-E_{b}}{N}\right)^{2}+\dfrac{\lambda_{F}}{3BN}\sum_{i=1}^{BN}\sum_{\alpha=1}^{3}\left(\hat{F}_{i,\alpha}-F_{i,\alpha}\right)^{2}+\dfrac{\lambda_{\sigma}}{9B}\sum_{b}^{B}\sum_{\alpha=1}^{3}\sum_{\beta=1}^{3}\left(\hat{\sigma}_{b,\alpha\beta}-\sigma_{b,\alpha\beta}\right)^{2}$
(4)

where $B$ is the batch size, $N$ is the number of atoms, $E_{b}$ is the DFT
energy, $\hat{E}_{b}$ is the predicted energy, $F_{i,\alpha}$ is the DFT force
of atom $i$ along direction $\alpha$, $\hat{F}_{i,\alpha}$ is the predicted
force of atom $i$ along direction $\alpha$, $\sigma_{b,\alpha\beta}$ is the
DFT stress tensor $\alpha\beta$ component and $\hat{\sigma}_{b,\alpha\beta}$
is the predicted stress tensor $\alpha\beta$ component. We used the scaling
factors $\lambda_{E}=1$, $\lambda_{F}=1$ and $\lambda_{\sigma}=10$ (with
energy, force and stress units in eV, $\text{eV}\cdot\text{\AA}^{-1}$ and
$\text{eV}\cdot\text{\AA}^{-3}$, respectively).

###  A.2 Hyperparameters

The scalar latent features were generated via a projection onto eight
trainable Bessel functions and the associated two-body multi layer perceptron
(MLP) consisted of four hidden layers of dimensions [128, 256, 512, 1024] with
SiLU nonlinearities [60]. A cutoff function with a cutoff radius of 6.5 Å and
a polynomial envelope function with $p=6$ was used. For the equivariant
features, a maximum rotation order $l_{max}=2$ with 16 features of even parity
was used. The model contained one interaction layer, where the latent MLP
contained three hidden layers of dimensions [1024, 1024, 1024] with SiLU
nonlinearities [60] and the embedding MLP for the tensor product operation
contained a single layer with no nonlinearity. The final output MLP was a
single layer with 128 nodes and no nonlinearity.

Training was performed with a batch size of 5 and the Adam optimizer was used
with an initial learning rate of 0.01, decaying by a factor of 0.5 whenever
the validation loss had not seen an improvement for 20 epochs.

##  Appendix B TDEP

((a))

((b))

Figure 8: Thermal conductivity of GeTe as a function of the number of
configurations used to extract second and third order force constants using
the TDEP method.

For TDEP calculations, we used cutoff radii of 12, 8 and 5 Å to extract
second, third and fourth-order IFCs, respectively. We verified that the
employed cutoff radii lead to converged results. We also tested the
convergence of the thermal conductivity values with respect to the number of
configurations for different temperature, as shown in Fig. 3. We used 300
configurations to extract the IFCs used to compute phonon dispersions and
thermal conductivities presented in Sec. IV.1.

In addition, for the rhombohedral phase, we used equilibrium atomic positions
computed as an average over all MD snapshots to avoid the spurious soft LO
phonon mode at $\Gamma$. For the cubic phase we used ideal atomic positions.

##  Appendix C Amorphous structure

The amorphous structure was generated by a melt-quench (MQ) MD simulation
using our GeTe MLIP. We used a 512-atom cubic GeTe supercell at the
experimental density (0.03327 atoms/$\text{\AA}^{3}$ [61]) as the starting
structure. The structure was first melted at 1800 K for 200 ps, then quenched
to 300 K over 200 ps in the NVT ensemble. Finally, the structure was
equilibrated in the NpT ensemble at 300 K for 1 ns. The density of the final
structure was 0.03328 atoms/$\text{\AA}^{3}$, very close to the experimental
value. The amorphous structure is shown in Fig. 9.

Figure 9: Amorphous GeTe structure generated by melt-quench MD simulation.

For the lattice thermal conductivity calculations, the amorphous structure was
equilibrated for 0.8 ns in the NpT ensemble at different temperatures before
computing the heat flux in the microcanonical ensemble. We computed the time-
averaged radial distribution function (RDF) $g(r)$ at the end of the NpT
simulation for 10 ps (see Fig. 10). Between 300 K and 450 K, the RDFs are very
similar and show no long range order beyond $6~{}\text{\AA}$. In contrast, at
550 K the RDF displays some peaks beyond $6~{}\text{\AA}$, which suggests a
crystallization of the structure.

Figure 10: Radial distribution function of the amorphous structure after
equilibration at different temperatures.

##  Appendix D Convergence of Green-Kubo simulations

Green-Kubo calculations require large system sizes and long simulation times
in order to obtain a converged thermal conductivity. We computed the lattice
thermal conductivity of the rhombohedral phase using supercells with 192, 360,
and 2880 atoms. For each system size, we performed 10 independent MD
simulations with different initial velocities. The averaged thermal
conductivity over the different MD trajectories is shown in Fig. 11, along
with the standard error (shaded region). The thermal conductivity computed for
the 2880-atom supercell exhibits larger standard errors compared to 360- and
192-atom supercells due to the shorter simulation times used. The lattice
thermal conductivities are extracted for correlation times where the averaged
thermal conductivity reaches a plateau (i.e., typically for $\tau>25$ ps).

Figure 11: Lattice thermal conductivity of rhombohedral GeTe as a function of
correlation time computed at 300 K using GK method for different system sizes.

##  Appendix E Comparison with previous theoretical works

Figure 12 shows the comparison of calculated lattice thermal conductivity of
crystalline GeTe using the BTE (circles) or GK (squares) method, compared to
previous theoretical works from Refs. [47, 56, 59]. It should be noted that
for $T\lesssim 650$ K, GeTe is in the rhombohedral ($R3m$) phase, while for
$T\gtrsim 650$ K, GeTe is in the cubic ($Fm\bar{3}m$) phase. In our
calculations, we used the structural parameters obtained from the NpT MD
simulations (which captures the phase transition at $T\approx 650$ K). The BTE
calculations shown in Fig. 12 include three-phonon processes using IFCs
extracted with the TDEP method. BTE thermal conductivities including three-
and four-phonon processes are shown in Fig. 7.

Figure 12: Lattice thermal conductivity of crystalline GeTe (rhombohedral
$R3m$ at the left, and cubic $Fm\bar{3}m$ phase at the right) as a function of
temperature calculated by the Green-Kubo (GK, squares) or the Boltzmann
transport equation (BTE, circles) methods and compared with previous
theoretical works of Ref. [47, 56, 59], all solving the BTE.

((a))

((b))

Figure 13: Evolution of LABEL:sub@fig:lattice-comp-a the lattice parameter $a$
and LABEL:sub@fig:lattice-comp-angle the angle $\alpha$ as a function of
temperature in the NpT MD simulations of crystalline GeTe (starting with a
$R3m$ supercell), compared against simulations from Ref. [62] and experimental
data from Ref. [36, 35, 37]. Simulated lattice parameters in
LABEL:sub@fig:lattice-comp-a were shifted by $-0.1~{}\text{\AA}$.

A similar approach was employed in Ref. [47] (light blue curve in Fig. 12) but
the reported values are significantly different compared to our results for
the rhombohedral phase. The primary source of discrepancy seems to be related
to the convergence with respect to the number of configurations used to
extract second and third order TDEP IFCs. In Ref. [47], the authors used only
24 configurations. According to our convergence test (see Fig. 8), such a
small number of configurations results in strongly underestimated thermal
conductivity values. Another source of discrepancy lies in the temperature
evolution of the structural parameters (see Fig. 13) which may also lead to
some differences in the computed lattice thermal conductivities. For the cubic
phase ($T\geq 700$ K), our results show a good agreement with those of Ref.
[47], presumably because less configurations are required to achieve
convergence compared to the rhombohedral symmetry.

In Ref. [56] (green curve in Fig. 12), the authors used 0 K IFCs and
structural parameters of $\alpha$-GeTe to compute the BTE thermal
conductivity. Therefore, their results exhibit a different temperature
evolution.

Finally, in Ref. [59], temperature-dependent IFCs were computed for the
$\beta$ phase of GeTe up to fourth order (i.e. three- and four-phonon
scattering processes) and BTE was solved under the relaxation time
approximation at 800 K. Their thermal conductivity value considering only
three-phonon scattering is significantly larger compared to the one from our
work or Ref. [47]. Including four-phonon scattering reduces the thermal
conductivity by the same factor as in our work (reduction by a factor 2.23).

## References

  * Krinner  _et al._ [2019] S. Krinner, S. Storz, P. Kurpiers, P. Magnard, J. Heinsoo, R. Keller, J. Lütolf, C. Eichler, and A. Wallraff, Engineering cryogenic setups for 100-qubit scale superconducting circuit systems, EPJ Quantum Technol. 6, 1–29 (2019).
  * Zhang  _et al._ [2022a] Q. Zhang, K. Deng, L. Wilkens, H. Reith, and K. Nielsch, Micro-thermoelectric devices, Nat. Electron. 5, 333 (2022a).
  * Zhang  _et al._ [2020] X. Zhang, Z. Bu, S. Lin, Z. Chen, W. Li, and Y. Pei, GeTe Thermoelectrics, Joule 4, 986 (2020).
  * Gutiérrez Moreno  _et al._ [2020] J. J. Gutiérrez Moreno, J. Cao, M. Fronzi, and M. H. N. Assadi, A review of recent progress in thermoelectric materials through computational methods, Mater. Renew 9, 16 (2020).
  * Broido  _et al._ [2007] D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A. Stewart, Intrinsic lattice thermal conductivity of semiconductors from first principles, Appl. Phys. Lett. 91, 231922 (2007).
  * Sklénard  _et al._ [2021] B. Sklénard, F. Triozon, C. Sabbione, L. Nistor, M. Frei, G. Navarro, and J. Li, Electronic and thermal properties of GeTe/Sb2Te3 superlattices by ab initio approach: Impact of van der waals gaps on vertical lattice thermal conductivity, Appl. Phys. Lett. 119, 201911 (2021).
  * Li  _et al._ [2022] R. Li, J.-X. Wang, E. Lee, and T. Luo, Physics-informed deep learning for solving phonon Boltzmann transport equation with large temperature non-equilibrium, Npj Comput. Mater. 8, 1–10 (2022).
  * Mingo [2006] N. Mingo, Anharmonic phonon flow through molecular-sized junctions, Phys. Rev. B 74, 125402 (2006).
  * Li  _et al._ [2009] J. Li, T. C. Au Yeung, C. H. Kam, Y. Peng, Q.-h. Chen, X. Zhao, and C. Q. Sun, Anharmonic phonon transport in atomic wire coupled by thermal contacts with surface bond reconstruction, J. Appl. Phys. 106, 014308 (2009).
  * Guo  _et al._ [2020] Y. Guo, M. Bescond, Z. Zhang, M. Luisier, M. Nomura, and S. Volz, Quantum mechanical modeling of anharmonic phonon-phonon scattering in nanostructures, Phys. Rev. B 102, 195412 (2020).
  * Müller-Plathe [1997] F. Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, J. Chem. Phys. 106, 6082 (1997).
  * Fan  _et al._ [2015] Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio, and A. Harju, Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations, Phys. Rev. B 92, 094301 (2015).
  * Isaeva  _et al._ [2019] L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach, Nat. Commun. 10, 3853 (2019).
  * Zhang  _et al._ [2022b] Z. Zhang, Y. Guo, M. Bescond, J. Chen, M. Nomura, and S. Volz, Heat conduction theory including phonon coherence, Phys. Rev. Lett. 128, 015901 (2022b).
  * Tristant  _et al._ [2019] D. Tristant, A. Cupo, X. Ling, and V. Meunier, Phonon anharmonicity in few-layer black phosphorus, ACS Nano 13, 10456–10468 (2019).
  * Allen and Feldman [1993] P. B. Allen and J. L. Feldman, Thermal conductivity of disordered harmonic solids, Phys. Rev. B 48, 12581 (1993).
  * Shenogin  _et al._ [2009] S. Shenogin, A. Bodapati, P. Keblinski, and A. J. H. McGaughey, Predicting the thermal conductivity of inorganic and polymeric glasses: The role of anharmonicity, J. Appl. Phys. 105, 034906 (2009).
  * Zhu and Shao [2022] X. Zhu and C. Shao, Effect of anharmonicity on the thermal conductivity of amorphous silica, Phys. Rev. B 106, 014305 (2022).
  * Carbogno  _et al._ [2017] C. Carbogno, R. Ramprasad, and M. Scheffler, Ab initio green-kubo approach for the thermal conductivity of solids, Phys. Rev. Lett. 118, 175901 (2017).
  * Verdi  _et al._ [2021] C. Verdi, F. Karsai, P. Liu, R. Jinnouchi, and G. Kresse, Thermal transport and phase transitions of zirconia by on-the-fly machine-learned interatomic potentials, npj Comput Mater 7, 1 (2021).
  * Korotaev  _et al._ [2019] P. Korotaev, I. Novoselov, A. Yanilkin, and A. Shapeev, Accessing thermal conductivity of complex compounds by machine learning interatomic potentials, Phys. Rev. B 100, 144308 (2019).
  * Sosso  _et al._ [2012] G. C. Sosso, D. Donadio, S. Caravati, J. Behler, and M. Bernasconi, Thermal transport in phase-change materials from atomistic simulations, Phys. Rev. B 86, 104301 (2012).
  * Schütt  _et al._ [2017] K. Schütt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko, and K.-R. Müller, SchNet: A continuous-filter convolutional neural network for modeling quantum interactions, in  _Advances in Neural Information Processing Systems_, Vol. 30 (Curran Associates, Inc., 2017).
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun. 13, 2453 (2022).
  * Schütt  _et al._ [2021] K. Schütt, O. Unke, and M. Gastegger, Equivariant message passing for the prediction of tensorial properties and molecular spectra, in  _Proceedings of the 38th International Conference on Machine Learning_ (PMLR, 2021) pp. 9377–9388.
  * Brandstetter  _et al._ [2022] J. Brandstetter, R. Hesselink, E. van der Pol, E. J. Bekkers, and M. Welling, Geometric and Physical Quantities Improve E(3) Equivariant Message Passing (2022), arXiv:2110.02905 [cs, stat].
  * Musaelian  _et al._ [2023] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nat Commun 14, 579 (2023).
  * Lankhorst  _et al._ [2005] M. H. R. Lankhorst, B. W. S. M. M. Ketelaars, and R. a. M. Wolters, Low-cost and nanoscale non-volatile memory concept for future silicon chips, Nature Mater 4, 347 (2005).
  * Jeong  _et al._ [2021] K. Jeong, H. Lee, C. Lee, L. H. Wook, H. Kim, E. Lee, and M.-H. Cho, Ferroelectric switching in gete through rotation of lone-pair electrons by electric field-driven phase transition, Appl. Mater. Today 24, 101122 (2021).
  * Zhang  _et al._ [2021] Q. Zhang, Z. Ti, Y. Zhu, Y. Zhang, Y. Cao, S. Li, M. Wang, D. Li, B. Zou, Y. Hou, P. Wang, and G. Tang, Achieving Ultralow Lattice Thermal Conductivity and High Thermoelectric Performance in GeTe Alloys via Introducing Cu2Te Nanocrystals and Resonant Level Doping, ACS Nano 15, 19345 (2021).
  * Jiang  _et al._ [2022] Y. Jiang, J. Dong, H.-L. Zhuang, J. Yu, B. Su, H. Li, J. Pei, F.-H. Sun, M. Zhou, H. Hu, J.-W. Li, Z. Han, B.-P. Zhang, T. Mori, and J.-F. Li, Evolution of defect structures leading to high ZT in GeTe-based thermoelectric materials, Nat. Commun. 13, 6087 (2022).
  * Picozzi [2014] S. Picozzi, Ferroelectric Rashba semiconductors as a novel class of multifunctional materials, Front. Phys. 2, 10.3389/fphy.2014.00010 (2014).
  * Rinaldi  _et al._ [2016] C. Rinaldi, J. C. Rojas-Sánchez, R. N. Wang, Y. Fu, S. Oyarzun, L. Vila, S. Bertoli, M. Asa, L. Baldrati, M. Cantoni, J.-M. George, R. Calarco, A. Fert, and R. Bertacco, Evidence for spin to charge conversion in GeTe(111), APL Mater. 4, 032501 (2016).
  * Varotto  _et al._ [2021] S. Varotto, L. Nessi, S. Cecchi, J. Sławińska, P. Noël, S. Petrò, F. Fagiani, A. Novati, M. Cantoni, D. Petti, E. Albisetti, M. Costa, R. Calarco, M. Buongiorno Nardelli, M. Bibes, S. Picozzi, J.-P. Attané, L. Vila, R. Bertacco, and C. Rinaldi, Room-temperature ferroelectric switching of spin-to-charge conversion in germanium telluride, Nat. Electron. 4, 740–747 (2021).
  * Chatterji  _et al._ [2015] T. Chatterji, C. M. N. Kumar, and U. D. Wdowik, Anomalous temperature-induced volume contraction in GeTe, Phys. Rev. B 91, 054110 (2015).
  * Sist  _et al._ [2018] M. Sist, H. Kasai, E. M. J. Hedegaard, and B. B. Iversen, Role of vacancies in the high-temperature pseudodisplacive phase transition in GeTe, Phys. Rev. B 97, 094116 (2018).
  * Chattopadhyay  _et al._ [1987] T. Chattopadhyay, J. X. Boucherle, and H. G. vonSchnering, Neutron diffraction study on the structural phase transition in GeTe, J. Phys. C: Solid State Phys. 20, 1431 (1987).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. i. t. Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271, 108171 (2022).
  * Thompson  _et al._ [2009] A. P. Thompson, S. J. Plimpton, and W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, J. Chem. Phys. 131, 154107 (2009).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010).
  * Lee  _et al._ [2024] S.-H. Lee, J. Li, V. Olevano, and B. Sklénard, Dataset related to article: Equivariant graph neural network interatomic potential for Green-Kubo thermal conductivity in phase change materials, 10.5281/zenodo.10688120 (2024).
  * Liu  _et al._ [2021] P. Liu, C. Verdi, F. Karsai, and G. Kresse, $\alpha\text{$-$}\beta$ phase transition of zirconium predicted by on-the-fly machine-learned force field, Phys. Rev. Mater. 5, 053804 (2021).
  * Wang  _et al._ [2021] C. Wang, J. Wu, Z. Zeng, J. Embs, Y. Pei, J. Ma, and Y. Chen, Soft-mode dynamics in the ferroelectric phase transition of GeTe, npj Comput Mater 7, 1 (2021).
  * Dangić  _et al._ [2021] D. Dangić, O. Hellman, S. Fahy, and I. Savić, The origin of the lattice thermal conductivity enhancement at the ferroelectric phase transition in GeTe, Npj Comput. Mater. 7, 1 (2021).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Steigmeier and Harbeke [1970] E. F. Steigmeier and G. Harbeke, Soft phonon mode and ferroelectricity in GeTe, Solid State Commun. 8, 1275 (1970).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Hellman and Abrikosov [2013] O. Hellman and I. A. Abrikosov, Temperature-dependent effective third-order interatomic force constants from first principles, Phys. Rev. B 88, 144301 (2013).
  * Bottin  _et al._ [2020] F. Bottin, J. Bieder, and J. Bouchet, a-TDEP: Temperature Dependent Effective Potential for Abinit – Lattice dynamic properties including anharmonicity, Comp. Phys. Comm. 254, 107301 (2020).
  * Fons  _et al._ [2010] P. Fons, A. V. Kolobov, M. Krbal, J. Tominaga, K. S. Andrikopoulos, S. N. Yannopoulos, G. A. Voyiatzis, and T. Uruga, Phase transition in crystalline GeTe: Pitfalls of averaging effects, Phys. Rev. B 82, 155209 (2010).
  * Kadlec  _et al._ [2011] F. Kadlec, C. Kadlec, P. Kužel, and J. Petzelt, Study of the ferroelectric phase transition in germanium telluride using time-domain terahertz spectroscopy, Phys. Rev. B 84, 205209 (2011).
  * Acharyya  _et al._ [2020] P. Acharyya, S. Roychowdhury, M. Samanta, and K. Biswas, Ultralow Thermal Conductivity, Enhanced Mechanical Stability, and High Thermoelectric Performance in (GeTe)1-2x(SnSe)x(SnS)x, J. Am. Chem. Soc. 142, 20502 (2020).
  * Ghosh  _et al._ [2020] K. Ghosh, A. Kusiak, P. Noé, M.-C. Cyrille, and J.-L. Battaglia, Thermal conductivity of amorphous and crystalline GeTe thin film at high temperature: Experimental and theoretical study, Phys. Rev. B 101, 214305 (2020).
  * Nath and Chopra [1974] P. Nath and K. L. Chopra, Thermal conductivity of amorphous and crystalline Ge and GeTe films, Phys. Rev. B 10, 3412 (1974).
  * Han  _et al._ [2022] Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, Fourphonon: An extension module to shengbte for computing four-phonon scattering rates and thermal conductivity, Comput. Phys. Commun. 270, 108179 (2022).
  * Xia and Chan [2018] Y. Xia and M. K. Y. Chan, Anharmonic stabilization and lattice heat transport in rocksalt $\beta$-GeTe, Appl. Phys. Lett. 113, 193902 (2018).
  * Hendrycks and Gimpel [2023] D. Hendrycks and K. Gimpel, Gaussian Error Linear Units (GELUs) (2023), arXiv:1606.08415 [cs].
  * Ghezzi  _et al._ [2011] G. E. Ghezzi, J. Y. Raty, S. Maitrejean, A. Roule, E. Elkaim, and F. Hippert, Effect of carbon doping on the structure of amorphous GeTe phase change material, Appl. Phys. Lett. 99, 151906 (2011).
  * Dangić  _et al._ [2018] D. Dangić, A. R. Murphy, E. D. Murray, S. Fahy, and I. Savić, Coupling between acoustic and soft transverse optical phonons leads to negative thermal expansion of GeTe near the ferroelectric phase transition, Phys. Rev. B 97, 224106 (2018).


# A quantum chemistry approach to linear vibro-polaritonic IR spectra with
perturbative electron-photon correlation

Eric W. Fischer<EMAIL_ADDRESS>Institut für Chemie, Humboldt-
Universität zu Berlin, Brook-Taylor-Straße 2, D-12489 Berlin, Germany
Institut für Chemie, Universität Potsdam, Karl-Liebknecht-Straße 24-25,
D-14476 Potsdam-Golm, Germany  Jan A. Syska  Institut für Chemie, Universität
Potsdam, Karl-Liebknecht-Straße 24-25, D-14476 Potsdam-Golm, Germany  Peter
Saalfrank<EMAIL_ADDRESS>Institut für Chemie, Universität
Potsdam, Karl-Liebknecht-Straße 24-25, D-14476 Potsdam-Golm, Germany  Institut
für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Straße 24-25,
D-14476 Potsdam-Golm, Germany

###### Abstract

In the vibrational strong coupling (VSC) regime, molecular vibrations and
resonant low-frequency cavity modes form light-matter hybrid states, named
vibrational polaritons, with characteristic IR spectroscopic signatures. Here,
we introduce a quantum chemistry based computational scheme for linear IR
spectra of vibrational polaritons in polyatomic molecules, which
perturbatively accounts for nonresonant electron-photon interactions under
VSC. Specifically, we formulate a cavity Born-Oppenheimer perturbation theory
(CBO-PT) linear response approach, which provides an approximate but
systematic description of such electron-photon correlation effects in VSC
scenarios, while relying on molecular ab initio quantum chemistry methods. We
identify relevant electron-photon correlation effects at second-order of CBO-
PT, which manifest as static polarizability-dependent Hessian corrections and
an emerging polarizability-dependent cavity intensity component providing
access to transmission spectra commonly measured in vibro-polaritonic
chemistry. Illustratively, we address electron-photon correlation effects
perturbatively in IR spectra of CO2 and Fe(CO)5 vibro-polaritonic models
qualitatively in sound agreement with non-perturbative CBO linear response
theory.

Introduction.— Vibrational polaritons are light-matter hybrid states formed
when molecular vibrational modes strongly interact with quantized modes of
optical low-frequency cavities[1] lying at the heart of the emerging field of
vibro-polaritonic chemistry[2, 3]. Experimentally, vibrational polaritons
exhibit characteristic spectroscopic signatures, specifically transitions to
upper and lower vibro-polaritonic states, which have been probed by both
linear[4, 5, 6, 7, 8, 9] and non-linear[10] infrared (IR) spectroscopic
techniques.

Computationally, vibro-polaritonic IR spectra are commonly obtained from
linear response theory, where two distinct routes exist: First, linear
response approaches based on effective ground state Pauli-Fierz Hamiltonians,
which are fully characterized by purely molecular properties such as the
molecular ground state potential energy surface (PES) and dipole moment,
easily accessible by means of standard quantum chemistry methods.[11, 12, 13,
14] Second, a linear response approach formulated in the cavity Born-
Oppenheimer (CBO) approximation[15, 16, 17, 18], which relies on an extended
electronic structure problem taking into account nonresonant interactions
between electrons and low-frequency cavity modes inducing electron-photon
correlation.[18] The CBO formulation is more complete but also more involved,
since both generalized cavity PES and dipole moments depend on both molecular
and cavity coordinates. From a linear response perspective, both approaches
are formulated in double-harmonic approximation, i.e., harmonic approximation
for molecular and cavity modes in combination with a linear approximation for
the ground state dipole moment. This is particularly beneficial for treating
vibrations and spectra of large molecules or molecular ensembles quantum
mechanically, since only Hessians and linearized dipole moments at equilibrium
structures are required.

Figure 1: Sketch of CO2 and Fe(CO)5 vibro-polaritonic model systems studied
with respect to their linear IR response in this work.

Recently, we showed that effective ground state and CBO formulations differ in
the way they account for electron-photon correlation.[19] Specifically, the
effective ground state formulation can be understood as additional
approximation to the CBO formulation, which we named crude CBO approximation,
that fully neglects electron-photon correlation in the VSC regime.[19] The
nonresonant nature of interactions between electrons and low-frequency cavity
modes characterized by distinct excitation energy scales of both subsystems,
subsequently motivated us to propose a perturbative connection between crude
CBO and correlated CBO formulations, denoted as CBO perturbation theory (CBO-
PT).[19] CBO-PT solves the extended CBO electronic structure problem
perturbatively and allows for a systematic correction of electron-photon
correlation in the crude CBO approach.[19] Further, CBO-PT is conceptually
similar to other perturbative approaches[20, 21], which however constitute a
further approximation[20, 19] or do not rely on the CBO formulation[21].

Here, we introduce an approximation to the non-perturbative CBO linear
response framework of Bonini and Flick[15] by combining CBO-PT with linear
response theory. This approach allows us to systematically correct vibro-
polaritonic IR spectra for electron-photon correlation effects while fully
relying on molecular properties, which are accessible by ab initio quantum
chemistry methods. We derive explicit expressions for electron-photon
correlation corrected vibro-polaritonic Hessian matrix elements and IR
intensities, which are both shown to be determined by the molecule’s static
polarizability, whose possible relevance has been noted only recently.[21, 22,
23] We illustratively apply the CBO-PT linear response approach to vibro-
polaritonic model systems of CO2 and Fe(CO)5 (cf. Fig.1) in line with
Ref.[15], and discuss electron-photon correlation effects in connection to
experimentally relevant transmission spectra[4, 7].

CBO Linear Response Theory.— We first recapitulate the main aspects of CBO
linear-response theory[15], which fully accounts for electron-photon
correlation in vibro-polaritonic IR spectra. We consider a molecular system
with $N_{n}$ nuclei under VSC with $2N_{c}$ quantized transverse field modes
of an IR cavity. We furthermore assume the CBO approximation to be valid,
i.e., non-adiabatic coupling to the excited state manifold is negligible.[18,
19] The CBO electronic ground state problem is described by an electron-photon
time-independent Schrödinger equation (TISE)[17, 18]

$\displaystyle\hat{H}_{ec}\,\ket{\Psi^{(ec)}_{0}(\underline{R},\underline{x})}$
$\displaystyle=E^{(ec)}_{0}(\underline{R},\underline{x})\,\ket{\Psi^{(ec)}_{0}(\underline{R},\underline{x})}\quad,$
(1)

with adiabatic electron-photon ground state, $\ket{\Psi^{(ec)}_{0}}$, and
ground state cavity potential energy surface (cPES), $E^{(ec)}_{0}$, which
parametrically depend on both nuclear, $\underline{R}$, and cavity
displacement coordinates, $\underline{x}$. The electron-photon Hamiltonian
reads[18]

$\displaystyle\hat{H}_{ec}$ $\displaystyle=\hat{H}_{e}+V_{c}+W_{c}\quad,$ (2)

with electronic Hamiltonian, $\hat{H}_{e}=\hat{T}_{e}+V_{\mathrm{coul}}$,
composed of the electronic kinetic energy, $\hat{T}_{e}$, and the molecular
Coulomb potential, $V_{\mathrm{coul}}$. Further,
$V_{c}=\frac{1}{2}\sum^{2N_{c}}_{\lambda,k}\,\omega^{2}_{k}\,x^{2}_{\lambda
k}$, is the harmonic cavity potential characterized by harmonic frequencies,
$\omega_{k}$, and displacement coordinates, $x_{\lambda k}$, with
polarization, $\lambda$, and mode index, $k$, respectively. The third term on
the right-hand side of Eq.(2) constitutes the light-matter interaction
potential

$\displaystyle W_{c}$
$\displaystyle=g_{0}\sum^{2N_{c}}_{\lambda,k}\omega_{k}\left(\underline{e}_{\lambda
k}\cdot\underline{d}_{en}\right)x_{\lambda
k}+\dfrac{g^{2}_{0}}{2}\sum^{2N_{c}}_{\lambda,k}\left(\underline{e}_{\lambda
k}\cdot\underline{d}_{en}\right)^{2},$ (3)

with light-matter interaction linear in the coupling constant, $g_{0}$, and
dipole self-energy (DSE) quadratic in $g_{0}$.[24] The coupling constant,
$g_{0}=\frac{1}{\sqrt{\epsilon_{0}V_{\mathrm{cav}}}}$, is only determined by
the cavity volume, $V_{\mathrm{cav}}$, and the permittivity, $\epsilon_{0}$.
It relates to a mode specific coupling constant,
$g_{k}=\sqrt{2\omega_{k}/\hbar}\,g_{0}$ (cf. Supplementary Information (SI),
Sec.S1). The light-matter interaction, $W_{c}$, in Eq.(3) is mediated by the
polarization-projected molecular dipole operator, $\underline{e}_{\lambda
k}\cdot\underline{d}_{en}$, with cavity polarization vector,
$\underline{e}_{\lambda k}$. Thus, $W_{c}$ couples the electronic subsystem to
low-frequency cavity modes and additionally induces electron-electron
interactions due to the DSE term.[18, 15, 19] In the non-interacting limit
with $g_{0}=0$, the light-matter interaction, $W_{c}$, vanishes such that,
$\hat{H}_{ec}=\hat{H}_{e}+V_{c}$, and Eq.(1) reduces to the electronic TISE
with a constant energy shift induced by the cavity potential, $V_{c}$.

In CBO linear response theory, the ground state cPES in Eq.(1) is harmonically
approximated[15]

$\displaystyle E^{(ec)}_{0}(\underline{R},\underline{x})$
$\displaystyle\approx
E^{(ec)}_{0}(\underline{R}_{0},\underline{x}_{0})+\frac{1}{2}\,\underline{C}^{T}\underline{\underline{H}}^{(ec)}_{0}\,\underline{C}\quad,$
(4)

around a minimum configuration, $(\underline{R}_{0},\underline{x}_{0})$, of
the cPES with coordinate vector,
$\underline{C}=(\underline{\tilde{R}},\underline{x})$, which collects
$N_{\mathrm{vib}}$ mass-weighted cartesian displacement coordinates,
$\tilde{R}_{i}=\sqrt{M_{i}}\,R_{i}$, with atomic mass, $M_{i}$, and $2N_{c}$
cavity displacement coordinates, $\underline{x}$. The vibro-polaritonic
Hessian, $\underline{\underline{H}}^{(ec)}_{0}$, contains second-order
derivatives of the cPES with respect to both mass-weighted molecular cartesian
and cavity displacement coordinates. It specifies a matrix eigenvalue
problem[15]

$\displaystyle\begin{pmatrix}\underline{\underline{H}}^{(ec)}_{QQ}&\underline{\underline{H}}^{(ec)}_{QC}\vspace{0.2cm}\\\
\underline{\underline{H}}^{(ec)}_{QC}&\underline{\underline{H}}^{(ec)}_{CC}\end{pmatrix}\begin{pmatrix}\underline{q}_{m}\vspace{0.3cm}\\\
\underline{c}_{m}\end{pmatrix}=\Omega^{2}_{m}\begin{pmatrix}\underline{q}_{m}\vspace{0.3cm}\\\
\underline{c}_{m}\end{pmatrix}\quad.$ (5)

where $\underline{\underline{H}}^{(ec)}_{0}$ is written as $2\times 2$-block
matrix with molecular ($QQ$), cavity ($CC$) and light-matter interaction
blocks ($QC$), respectively. Eigenvectors,
$(\underline{q}_{m},\underline{c}_{m})^{T}$, resemble vibro-polaritonic normal
modes with molecular, $\underline{q}_{m}$, and cavity components,
$\underline{c}_{m}$, characterized by corresponding harmonic vibro-polaritonic
frequencies, $\Omega_{m}$. The eigensystem of Eq.(5) provides access to linear
vibro-polaritonic IR spectra

$\displaystyle\sigma_{\mathrm{IR}}(\hbar\omega)$
$\displaystyle=\sum^{N_{p}}_{m}I_{m}\,\delta(\hbar\omega-\hbar\Omega_{m})\quad,$
(6)

for $N_{p}=N_{\mathrm{vib}}+2N_{c}$ vibro-polaritonic normal modes. IR
intensities, $I_{m}=\sum_{\kappa}|Z_{\kappa m}|^{2}$, are determined by
cartesian components of mode effective charges

$\displaystyle Z_{\kappa m}$ $\displaystyle=Z^{(Q)}_{\kappa m}+Z^{(C)}_{\kappa
m}\quad,\quad\kappa=x,y,z\quad,$ (7)

with molecular and cavity mode contributions[15]

$\displaystyle Z^{(Q)}_{\kappa m}$
$\displaystyle=\sum^{N_{\mathrm{vib}}}_{i}\dfrac{\partial
D^{(ec)}_{00,\kappa}}{\partial\tilde{R}_{i}}\,q_{mi}\quad,\vspace{0.2cm}$ (8)
$\displaystyle Z^{(C)}_{\kappa m}$
$\displaystyle=\sum^{2N_{c}}_{\lambda,k}\dfrac{\partial
D^{(ec)}_{00,\kappa}}{\partial x_{\lambda k}}\,c_{m\lambda k}\quad,$ (9)

where, $q_{mi}$ and $c_{m\lambda k}$, resemble elements of $\underline{q}_{m}$
and $\underline{c}_{m}$ in Eq.(5), respectively. The IR intensity of the
$m\mathrm{th}$-vibro-polaritonic transition decomposes via Eq.(7) into

$\displaystyle I_{m}$ $\displaystyle=\sum_{\kappa}\left[(Z^{(Q)}_{\kappa
m})^{2}+(Z^{(C)}_{\kappa m})^{2}+2Z^{(Q)}_{\kappa m}Z^{(C)}_{\kappa
m}\right]\quad,$ (10)

such that vibro-polaritonic IR spectra, $\sigma_{\mathrm{IR}}$, exhibit a
molecular, $\sigma^{(Q)}_{\mathrm{IR}}$, a cavity,
$\sigma^{(C)}_{\mathrm{IR}}$, and an additional mixed contribution,
$\sigma^{(X)}_{\mathrm{IR}}$, respectively. Further, $Z_{\kappa m}$ is
determined by cartesian components, $D^{(ec)}_{00,\kappa}$, of a generalized
permanent dipole moment

$\displaystyle\underline{D}^{(ec)}_{00}(\underline{R},\underline{x})$
$\displaystyle=\braket{\Psi^{(ec)}_{0}(\underline{R},\underline{x})}{\underline{d}_{en}}{\Psi^{(ec)}_{0}(\underline{R},\underline{x})}_{\underline{r}}\quad,$
(11)

which is evaluated with respect to the electron-photon ground state,
$\ket{\Psi^{(ec)}_{0}}$, where integration over electronic coordinates is
indicated by $\braket{\dots}_{\underline{r}}$.

Perturbative CBO Linear Response Theory.— We now introduce a systematic
approximation to the CBO linear response approach[15] based on cavity Born-
Oppenheimer perturbation theory (CBO-PT)[19]. CBO-PT approximates nonresonant
interactions between electrons and low-frequency cavity modes and therefore
electron-photon correlation by treating the light-matter interaction
potential, $W_{c}$, in Eq.(3) as perturbation of the electronic subsystem[19]

$\displaystyle\hat{H}_{ec}$ $\displaystyle=\hat{H}_{0}+\lambda\,W_{c}\quad,$
(12)

with zeroth-order Hamiltonian, $\hat{H}_{0}=\hat{H}_{e}+V_{c}$, and formal
perturbation parameter, $\lambda$. In $\hat{H}_{0}$, the cavity potential,
$V_{c}$, is a constant with respect to the bare electronic problem, such that
zeroth-order states are given by bare adiabatic electronic states,
$\ket{\Psi^{(e)}_{\mu}}$. CBO-PT is motivated by distinct excitation energy
scales of electronic and low-frequency cavity mode subsystems, and provides a
systematic approximation of electron-photon correlation in light-matter hybrid
systems under VSC by perturbatively solving the electron-photon TISE (1).[19]

At $n$th-order of CBO-PT, abbreviated as CBO-PT($n$), the approximate cPES
reads, $E^{(ec)}_{0}\approx\sum^{n}_{k=0}\lambda^{k}E^{(k)}_{0}$ (cf. SI,
Sec.S1), which straightforwardly leads to the corresponding Hessian

$\displaystyle\underline{\underline{H}}^{(ec)}_{0}\approx\sum^{n}_{k=0}\begin{pmatrix}\underline{\underline{H}}^{(k)}_{QQ}&\underline{\underline{H}}^{(k)}_{QC}\vspace{0.2cm}\\\
\underline{\underline{H}}^{(k)}_{QC}&\underline{\underline{H}}^{(k)}_{CC}\end{pmatrix}\quad,$
(13)

here given in block-matrix form in analogy to Eq.(5). We restrict our
discussion to perturbation order $n\leq 2$, i.e., up to CBO-PT(2), sufficient
to capture leading-order electron-photon correlation corrections[19]. In
contrast to the non-perturbative CBO approach[15], we work directly in a
representation employing mass-weighted molecular normal modes, which
constitute the CBO-PT(0) reference and relate to coordinates,
$\underline{Q}=(Q_{1},\dots,Q_{\mathrm{vib}})$. Note, molecular normal modes
are not strictly extractable from Eq.(5) due to light-matter interaction
unless $g_{0}=0$. Linear vibro-polaritonic IR spectra in CBO-PT($n$) are given
by

$\displaystyle\sigma^{(n)}_{\mathrm{IR}}(\hbar\omega)$
$\displaystyle=\sum^{N_{p}}_{m}I^{(n)}_{m}\,\delta(\hbar\omega-\hbar\Omega^{(n)}_{m})\quad,$
(14)

with approximate vibro-polaritonic normal-mode frequencies,
$\Omega^{(n)}_{m}$, and IR intensities,
$I^{(n)}_{m}=\sum_{\kappa}|Z^{(n)}_{\kappa m}|^{2}$, with, $Z^{(n)}_{\kappa
m}=Z^{(Q,n)}_{\kappa m}+Z^{(C,n)}_{\kappa m}$. The approximate mode effective
charge, $Z^{(n)}_{\kappa m}$, is determined by components (cf. SI, Sec.S1A for
derivation)

$\displaystyle Z^{(Q,n)}_{\kappa m}$
$\displaystyle=\sum^{N_{\mathrm{vib}}}_{i}\sqrt{\dfrac{\hbar}{2\omega_{i}}}\dfrac{\partial
D^{(n)}_{00,\kappa}}{\partial Q_{i}}\,q_{mi}\quad,\vspace{0.2cm}$ (15)
$\displaystyle Z^{(C,n)}_{\kappa m}$
$\displaystyle=\sum^{2N_{c}}_{\lambda,k}\sqrt{\dfrac{\hbar}{2\omega_{k}}}\dfrac{\partial
D^{(n)}_{00,\kappa}}{\partial x_{\lambda k}}\,c_{m\lambda k}\quad,$ (16)

which result from a linearization of the approximate generalized permanent
dipole moment in the normal mode representation

$\displaystyle\underline{D}^{(n)}_{00}(\underline{R},\underline{x})$
$\displaystyle=\braket{\Phi^{(n-1)}_{0}(\underline{R},\underline{x})}{\underline{d}_{en}}{\Phi^{(n-1)}_{0}(\underline{R},\underline{x})}_{\underline{r}}\,.$
(17)

The latter expression is determined by the perturbatively approximated
adiabatic electron-photon ground state,
$\ket{\Phi^{(n-1)}_{0}}=\sum^{n-1}_{k=0}\lambda^{k}\ket{\Psi^{(k)}_{0}}$ (cf.
SI, Sec.S1). Note, in contrast to Eq.(7), both $Z^{(Q,n)}_{\kappa m}$ and
$Z^{(C,n)}_{\kappa m}$ depend on prefactors with molecular normal mode,
$\omega_{i}$, and cavity frequencies, $\omega_{k}$, which follow from the
normal mode representation (cf. SI, Sec.S1A).

We now provide explicit expressions for Hessian matrix elements and IR
intensities up to CBO-PT(2) as derived in SI, Sec.S1. The sum of CBO-PT(0) and
CBO-PT(1) Hessians, which do not account for electron-photon correlation, is
determined by elements (cf. SI, Sec.S1B and S1C)

$\displaystyle\left(\underline{\underline{H}}^{(0)}_{QQ}+\underline{\underline{H}}^{(1)}_{QQ}\right)_{ij}=$
$\displaystyle\,\omega^{2}_{i}\,\delta_{ij}+g^{2}_{0}N_{c}\sum^{2}_{\lambda}d^{(i)}_{\lambda}d^{(j)}_{\lambda}\,,\vspace{0.2cm}$
(18)
$\displaystyle\left(\underline{\underline{H}}^{(0)}_{QC}+\underline{\underline{H}}^{(1)}_{QC}\right)_{i,\lambda
k}=$ $\displaystyle\,g_{0}\,\omega_{k}\,d^{(i)}_{\lambda}\quad,\vspace{0.2cm}$
(19)
$\displaystyle\left(\underline{\underline{H}}^{(0)}_{CC}+\underline{\underline{H}}^{(1)}_{CC}\right)_{\lambda
k,\lambda^{\prime}k^{\prime}}=$
$\displaystyle\,\omega^{2}_{k}\,\delta_{kk^{\prime}}\delta_{\lambda\lambda^{\prime}}\quad.$
(20)

Here, the DSE term induces inter-normal-mode couplings in Eq.(18). Light-
matter interaction elements in Eq.(19) relate to polarization-projected
permanent dipole moment derivatives,
$d^{(i)}_{\lambda}=\partial_{Q_{i}}\left(\underline{e}_{\lambda}\cdot\underline{d}_{00}\right)|_{Q_{i}=0}$.
We assume here identical polarization vectors for all modes with index, $k$.
The CBO-PT(1) IR spectrum, $\sigma^{(1)}_{\mathrm{IR}}$, is fully determined
by the zeroth-order adiabatic electronic ground state, $\ket{\Psi^{(e)}_{0}}$,
such that, $\underline{D}^{(1)}_{00}=\underline{d}_{00}$ in Eq.(17) with
first-order IR intensity given by

$\displaystyle I^{(Q,1)}_{m}$ $\displaystyle=\sum_{\kappa}|Z^{(Q,1)}_{\kappa
m}|^{2}\quad,\vspace{0.2cm}$ (21) $\displaystyle Z^{(Q,1)}_{\kappa m}$
$\displaystyle=\displaystyle\sum^{N_{\mathrm{vib}}}_{i}\sqrt{\dfrac{\hbar}{2\omega_{i}}}\dfrac{\partial
d_{00,\kappa}}{\partial Q_{i}}\,q^{(1)}_{mi}\quad,$ (22)

which reduces in the non-interacting limit with $g_{0}=0$ to
$q^{(1)}_{mi}=\delta_{mi}$ and the molecular intensity,
$I^{(Q,1)}_{m}=I^{(Q,1)}_{i}$. Notably, in CBO-PT(1) the vibro-polaritonic IR
spectrum is fully determined by molecular intensity components,
$I^{(Q,1)}_{m}$, since the cavity mode effective charge vanishes,
$Z^{(C,1)}_{\kappa m}=0$. Thus, vibro-polaritonic contributions enter only via
cavity-induced linear combinations of molecular normal mode components,
$q^{(1)}_{mi}$, and the CBO-PT(1) IR spectrum, $\sigma^{(1)}_{\mathrm{IR}}$,
exclusively probes the effective matter response of the light-matter hybrid
system under VSC. Further, CBO-PT(2) Hessian matrix element corrections in
Franck-Condon approximation are given by (cf. SI, Sec.S1D and S2)

$\displaystyle\left(\underline{\underline{H}}^{(2)}_{QQ}\right)_{ij}$
$\displaystyle=-\dfrac{g^{4}_{0}}{4}N^{2}_{c}\sum^{2}_{\lambda,\lambda^{\prime}}d^{(i)}_{\lambda}\alpha^{0}_{\lambda\lambda^{\prime}}d^{(j)}_{\lambda^{\prime}}\quad,\vspace{0.2cm}$
(23)
$\displaystyle\left(\underline{\underline{H}}^{(2)}_{QC}\right)_{i,\lambda k}$
$\displaystyle=-\dfrac{g^{3}_{0}}{2}\omega_{k}\,N_{c}\sum^{2}_{\lambda^{\prime}}\alpha^{0}_{\lambda\lambda^{\prime}}d^{(i)}_{\lambda^{\prime}}\quad,\vspace{0.2cm}$
(24) $\displaystyle\left(\underline{\underline{H}}^{(2)}_{CC}\right)_{\lambda
k,\lambda^{\prime}k^{\prime}}$
$\displaystyle=-g^{2}_{0}\,\omega_{k}\,\omega_{k^{\prime}}\,\alpha^{0}_{\lambda\lambda^{\prime}}\quad,$
(25)

which are all determined by polarization-projected ground state (superscript
$0$) static polarizability tensor elements,
$\alpha^{0}_{\lambda\lambda^{\prime}}$. Notably, the leading-order correction
in Eq.(25) enters the cavity block and exhibits a matter-induced perturbative
coupling of distinct cavity modes[19] in agreement with non-perturbative CBO
linear response theory[15].

Figure 2: Linear vibro-polaritonic IR spectra for selected single-molecule
models under VSC with a single cavity mode at coupling strength,
$g_{0}=0.03\,\sqrt{E_{h}}/ea_{0}$. Top-row: z-polarized IR spectra of the
antisymmetric CO2-stretching mode under VSC with a single cavity mode,
$\hbar\omega_{c}=\hbar\omega_{\mathrm{as}}=2400\,\mathrm{cm}^{-1}$, with a)
CBO-PT(1) and CBO-PT(2) IR spectra, $\sigma^{(1)}_{\mathrm{IR}}(\hbar\omega)$
and $\sigma^{(2)}_{\mathrm{IR}}(\hbar\omega)$, besides bare molecular
spectrum, $\sigma^{(0)}_{\mathrm{IR}}(\hbar\omega)$, b) molecular, cavity and
mixed CBO-PT(2) contributions,
$\sigma^{(Q,2)}_{\mathrm{IR}}(\hbar\omega),\sigma^{(C,2)}_{\mathrm{IR}}(\hbar\omega)$
and $\sigma^{(X,2)}_{\mathrm{IR}}(\hbar\omega)$ and c) comparison of
normalized $\sigma^{(2)}_{\mathrm{IR}}(\hbar\omega)$ and
$\sigma^{(C,2)}_{\mathrm{IR}}(\hbar\omega)$. Bottom-row: Polarization-averaged
linear vibro-polaritonic IR spectra for the CO-stretch band of a single
Fe(CO)5 molecule under VSC with a single cavity mode,
$\hbar\omega_{c}=\hbar\omega_{e^{\prime}}=2052\,\mathrm{cm}^{-1}$ with d)-f)
providing same information as a)-c). Stick spectra in d) resemble individual
polarization-dependent contributions to $\bar{\sigma}^{(2)}_{\mathrm{IR}}$
(cf. main text for details). The respective cavity mode frequency is indicated
by a grey vertical line in all panels.

The CBO-PT(2) IR spectrum, $\sigma^{(2)}_{\mathrm{IR}}(\hbar\omega)$, is
determined by second-order vibro-polaritonic frequencies, $\Omega^{(2)}_{m}$,
and intensities derived from second-order mode effective charges as

$\displaystyle I^{(2)}_{m}$ $\displaystyle=\sum_{\kappa}|Z^{(2)}_{\kappa
m}|^{2}\quad,\vspace{0.2cm}$ (26) $\displaystyle Z^{(Q,2)}_{\kappa m}$
$\displaystyle=\sum^{N_{\mathrm{vib}}}_{i}\sqrt{\dfrac{\hbar}{2\omega_{i}}}\dfrac{\partial
D^{(Q,2)}_{00,\kappa}}{\partial Q_{i}}\,q^{(2)}_{mi}\quad,\vspace{0.2cm}$ (27)
$\displaystyle Z^{(C,2)}_{\kappa m}$
$\displaystyle=\sum^{2N_{c}}_{\lambda,k}\sqrt{\dfrac{\hbar}{2\omega_{k}}}\dfrac{\partial
D^{(C,2)}_{00,\kappa}}{\partial x_{\lambda k}}\,c^{(2)}_{m\lambda k}\quad,$
(28)

which contains both molecular and cavity contributions with dipole derivatives
explicitly given by (cf. SI Sec.S3)

$\displaystyle\dfrac{\partial D^{(Q,2)}_{00,\kappa}}{\partial Q_{i}}$
$\displaystyle=d^{(i)}_{00,\kappa}-\dfrac{g^{2}_{0}}{2}N_{c}\sum^{2}_{\lambda}\alpha^{0}_{\kappa\lambda}\,d^{(i)}_{\lambda}\quad,\vspace{0.2cm}$
(29) $\displaystyle\dfrac{\partial D^{(C,2)}_{00,\kappa}}{\partial x_{\lambda
k}}$ $\displaystyle=-g_{0}\,\omega_{k}\,\alpha^{0}_{\kappa\lambda}\quad.$ (30)

with second-order dipole moment, $\underline{D}^{(2)}_{00}$, following from
Eq.(17). Note, in the non-interacting limit, $q^{(2)}_{mi}=\delta_{mi}$ and
$\frac{\partial D^{(C,2)}_{00,\kappa}}{\partial x_{\lambda k}}=0$, such that
the CBO-PT(2) IR spectrum reduces to the molecular spectrum. In Eq.(29) and
(30), we find a static polarization-dependent correction of the molecular
intensity and in particular a non-vanishing cavity component, which is
determined by the light-matter interaction strength, $g_{0}$, and the static
polarizability, $\alpha^{0}_{\kappa\lambda}$, of the molecular subsystem.
Notably, the related cavity dipole component, $D^{(C,2)}_{00,\kappa}$, is
linear in $g_{0}$ and $x_{\lambda k}$ (cf. SI, Sec.S3) in agreement with
numerical results reported for non-perturbative CBO linear response theory
under VSC (cf. supporting information of Ref.[15]) In contrast to CBO-PT(1),
the non-vanishing CBO-PT(2) cavity intensity component allows for addressing
transmission spectra, which exclusively probe the cavity response,
$\sigma^{(C)}_{\mathrm{IR}}$, of the light-matter hybrid system commonly
measured in vibro-polaritonic chemistry.[4, 7, 8] In this context, we like to
emphasize, that Eq.(30) solely emerges due to the fact that we
(perturbatively) account for electron-photon correlation effects in the
description of the vibro-polaritonic system by exclusively relying on
molecular ground state information. The latter is obtainable from standard
quantum chemistry with polarizabilities being practically accessible by means
of response theory.

IR spectra under VSC. — In the remainder of this work, we discuss IR spectra
of CO2 and Fe(CO)5 vibro-polaritonic models[15] obtained from second-order
CBO-PT linear response theory with a focus on electron-photon correlation
effects. We consider polarization-averaged vibro-polaritonic IR spectra

$\displaystyle\bar{\sigma}^{(n)}_{\mathrm{IR}}(\hbar\omega)$
$\displaystyle=\dfrac{1}{K}\sum^{K}_{\alpha}\sigma^{(n)}_{\alpha}(\hbar\omega)\quad,\quad
n=0,1,2\quad,$ (31)

to mimic random molecular orientation, where the sum runs over all unique
polarization components, $\alpha=x,y,z$, with $K=1,2,3$, depending on the
number of contributing terms. Molecular structures, harmonic frequencies,
dipole derivatives and static polarizability tensors elements are obtained
from density functional theory (DFT) using the TPSSh[25, 26] and B3LYP[27]
functionals and a Def2TZVP basis set[28, 29]. Computational details are
provided in SI, Sec.S4A. We furthermore consider spontaneous emission-induced
Lorentzian peak broadening (and neglect other dissipative channels) with FWHM,
$\kappa=41\,\mathrm{cm}^{-1}$ (cf. SI, Sec.S4B) for an experimentally
motivated cavity quality factor, $Q=\frac{\hbar\omega_{c}}{\kappa}$[30, 31],
with $Q=50$ (CO2) and $Q=59$ (Fe(CO)5).[8]

CO2 under VSC.— We first discuss the antisymmetric stretching mode of CO2
strongly coupled to a single resonant cavity mode,
$\hbar\omega_{c}=\hbar\omega_{\mathrm{as}}$, with
$\hbar\omega_{\mathrm{as}}=2400\,\mathrm{cm}^{-1}$ (TPSSh/Def2TZVP) at
coupling strength, $g_{0}=0.03\,\sqrt{E_{h}}/ea_{0}$. Here, the IR spectrum
reduces to,
$\sigma^{(n)}_{\mathrm{IR}}(\hbar\omega)=\sigma^{(n)}_{z}(\hbar\omega)$, since
only the dipole component along the molecular z-axis contributes. In Fig.2a,
CBO-PT(1) and CBO-PT(2) IR spectra are shown for CO2 under VSC besides the
bare molecular spectrum, $\sigma^{(0)}_{\mathrm{IR}}(\hbar\omega)$. We observe
a Rabi-splitting characterized by lower and upper vibro-polaritonic
transitions for both CBO-PT(1) and CBO-PT(2). Uncorrelated CBO-PT(1) spectra
exhibit nearly identical molecular and photonic contributions for both
polariton states, contrasted by a dominantly photonic lower and dominantly
molecular upper polariton state in correlated CBO-PT(2) spectra in agreement
with non-perturbative CBO results[15] (cf. SI, Sec.S4B). Further, the
correlated CBO-PT(2) spectrum exhibits a pronounced peak asymmetry in favour
of the lower polariton transition, a spectral red-shift and a pronounced
asymmetric peak splitting relative to the cavity frequency. Rabi splittings,
$\Omega^{(1)}_{R}=121\,\mathrm{cm}^{-1}$ and
$\Omega^{(2)}_{R}=134\,\mathrm{cm}^{-1}$, indicate a slightly stronger
effective light-matter interaction in presence of electron-photon correlation,
despite the fact that both results have been obtained at identical light-
matter coupling.

Figure 3: Linear vibro-polaritonic IR spectra for CO2 ensemble models under
VSC with a single cavity mode,
$\hbar\omega_{c}=\hbar\omega_{\mathrm{as}}=2400\,\mathrm{cm}^{-1}$. a)
Ensemble and effective CBO-PT(n) IR spectra,
$\sigma^{(n)}_{\mathrm{ens}}(\hbar\omega)$ and
$\sigma^{(n)}_{\mathrm{eff}}(\hbar\omega)$ for $n=,1,2$, b) molecular, cavity
and mixed CBO-PT(2) contributions,
$\sigma^{(Q,2)}_{\mathrm{ens}}(\hbar\omega),\sigma^{(C,2)}_{\mathrm{ens}}(\hbar\omega)$
and $\sigma^{(X,2)}_{\mathrm{ens}}(\hbar\omega)$ for explicit ensemble and c)
comparison of normalized $\sigma^{(2)}_{\mathrm{ens}}(\hbar\omega)$ and
$\sigma^{(C,2)}_{\mathrm{ens}}(\hbar\omega)$ for explicit ensemble. The
respective cavity mode frequency is indicated by a grey vertical line.

We now turn to the details of the CBO-PT(2) IR spectrum, Fig.2b, where we
compare molecular, $\sigma^{(Q,2)}_{\mathrm{IR}}$, cavity,
$\sigma^{(C,2)}_{\mathrm{IR}}$, and mixed contributions,
$\sigma^{(X,2)}_{\mathrm{IR}}$, following from Eq.(10). Although,
$\sigma^{(Q,2)}_{\mathrm{IR}}$ is dominant relative to
$\sigma^{(C,2)}_{\mathrm{IR}}$, $\sigma^{(X,2)}_{\mathrm{IR}}$ determines the
peak asymmetry of $\sigma^{(2)}_{\mathrm{IR}}$ by “shifting” intensity from
the upper to the lower vibro-polaritonic transition. This effect directly
translates to the sign of the photonic contribution in linear combinations
forming upper (+) and lower (-) polariton states (cf. SI, Sec.S4B): Mixed
intensity contributions satisfy,
$I^{(X,2)}_{\mp}\propto\pm|Z^{(Q,2)}_{\kappa\mp}||Z^{(C,2)}_{\kappa\mp}|$,
since the dipole derivative in Eq.(30) determining $Z^{(C,2)}_{\kappa\mp}$ is
negative. Further, we recall that IR transmission spectra of vibro-polaritonic
systems measured in micro-fluidic optical cavities probe the cavity response
of the light-matter hybrid system.[4, 7, 8] The cavity response translates
here into the cavity component, $\sigma^{(C,2)}_{\mathrm{IR}}$, which appears
in CBO-PT(2) but is not accounted for in uncorrelated CBO-PT(1) models[15]. In
Fig.2c, we compare full and cavity CBO-PT(2) IR spectra, where the intensity
has been normalized for reasons of comparison. We find nearly identical
results with only minor intensity deviations for the lower polariton peak in
the spectral region of interest. However, when the full IR spectrum is
considered, then $\sigma^{(2)}_{\mathrm{IR}}$ still contains molecular
contributions related to CO2 bending modes at low frequencies (cf. SI,
Sec.S4B). Those peaks are absent in $\sigma^{(C,2)}_{\mathrm{IR}}$ since
bending modes do not couple to the “high frequency” cavity mode.

Fe(CO)5 under VSC.— As a second example, we discuss the CO-stretch band of
Fe(CO)5 under VSC, which constitutes a both experimentally and theoretically
relevant molecular multi-mode example.[8, 15, 32] The CO-stretch band
decomposes into a doubly-degenerate equatorial normal mode (symmetry
$e^{\prime}$ in molecular point group $D_{3h}$) and a single axial normal mode
($a^{\prime\prime}_{2}$ in $D_{3h}$) with harmonic frequencies,
$\hbar\omega_{e^{\prime}}=2052\,\mathrm{cm}^{-1}$ and
$\hbar\omega_{a^{\prime\prime}_{2}}=2078\,\mathrm{cm}^{-1}$ (TPSSh/Def2TZVP,
cf. SI, Sec.S4A). We consider a single cavity mode tuned resonant to the
equatorial $e^{\prime}$-modes, $\hbar\omega_{c}=\hbar\omega_{e^{\prime}}$, at
light-matter interaction, $g_{0}=0.03\,\sqrt{E_{h}}/ea_{0}$. Here,
$e^{\prime}$-modes are $(x,z)$-polarized, while the
$a^{\prime\prime}_{2}$-mode is $y$-polarized. Therefore, all three molecular
modes can enter VSC depending on the cavity mode polarization, which
constitutes a scenario with both strongly coupled and decoupled molecular
modes. In Fig.2d, polarization-averaged CBO-PT(1) and CBO-PT(2) IR spectra are
shown besides the bare molecular CBO-PT(0) spectrum. Note, the
$a^{\prime\prime}_{2}$-peak in the latter is only visible as blue-shifted
shoulder relative to the dominant transition for the herein chosen peak
broadening. Further, we display CBO-PT(2) stick spectra to highlight the
polarization-dependence of the individual spectral contributions. According to
the averaged character, we only discuss effective polaritonic signatures
resulting from broadening effects. We find the $e^{\prime}$-peak to split
under VSC but in contrast to CO2, we find here additionally a pronounced
molecular contribution related to uncoupled molecular modes. Note, those do
not refer to “dark states” as present in the ensemble scenario below. The
correlated CBO-PT(2) IR spectrum exhibits a strong red-shift and significantly
more intense peaks relative to the uncorrelated CBO-PT(1) result. Both share a
slight peak asymmetry in favour of the upper polariton peak. Further, a
comparison of Rabi splittings, $\Omega^{(1)}_{R}=169\,\mathrm{cm}^{-1}$ and
$\Omega^{(2)}_{R}=228\,\mathrm{cm}^{-1}$, indicates also here an enhanced
effective light-matter interaction in presence of electron-photon correlation.
By analysing, $\bar{\sigma}^{(2)}_{\mathrm{IR}}$, in terms of,
$\bar{\sigma}^{(Q,2)}_{\mathrm{IR}},\bar{\sigma}^{(C,2)}_{\mathrm{IR}}$ and
$\bar{\sigma}^{(X,2)}_{\mathrm{IR}}$, reveals a trend similar to the CO2
example (cf. Fig.2e): $\bar{\sigma}^{(Q,2)}_{\mathrm{IR}}$ is the dominant
contribution relative to $\bar{\sigma}^{(C,2)}_{\mathrm{IR}}$, but the cross
term again “shifts” intensity from the upper to the lower vibro-polaritonic
transition. Importantly, in contrast to the CO2 example,
$\bar{\sigma}^{(Q,2)}_{\mathrm{IR}}$ contains here in addition to the matter
response of the light-matter hybrid system bare molecular contributions, which
relate to uncoupled CO-stretching modes. In order to access the light-matter
response without bare molecular contributions, we compare normalized,
$\bar{\sigma}^{(2)}_{\mathrm{IR}}$ and $\bar{\sigma}^{(C,2)}_{\mathrm{IR}}$,
in Fig.2f: Here, we find the upper polariton peak to significantly differ in
intensity for $\bar{\sigma}^{(C,2)}_{\mathrm{IR}}$, which directly relates to
the absence of additional molecular contributions of uncoupled CO-modes
present in $\bar{\sigma}^{(2)}_{\mathrm{IR}}$. Thus, correlation corrected
CBO-PT(2) IR spectra allow for avoiding artificial molecular intensity
contributions related to uncoupled molecular modes in vibrational multi-mode
systems under VSC, which are not probed in vibro-polaritonic transmission
spectroscopy.[4, 7, 8]

Collective VSC effects in CO2 ensembles.— We finally address collective strong
coupling effects in two different molecular ensemble models of CO2 under VSC
(cf. SI, Sec.5), which have also been studied in Ref.[15]: First, an explicit
ensemble of $M$ molecules in the dilute gas limit[33], which interact only via
DSE contributions at light-matter coupling, $g_{0}=0.01\,\sqrt{E_{h}}/ea_{0}$.
Second, an effective “ensemble” model, which contains a single CO2 molecule
under VSC with a single cavity mode characterized by an effective enhanced
interaction constant, $g_{\mathrm{eff}}=\sqrt{M}\,g_{0}$. In Fig.3a, we
compare ensemble, $\sigma^{(n)}_{\mathrm{ens}}$, and effective IRspectra,
$\sigma^{(n)}_{\mathrm{eff}}$, for CBO-PT(1) and CBO-PT(2) for $M=20$ parallel
aligned CO2 molecules. Intensities have been scaled by a factor of $M^{-1}$
for reasons of comparison. We find an excellent agreement between spectra
obtain at both levels of CBO-PT in agreement with Ref.[15]. Further, for the
ensemble model, we observe qualitatively identical results for the different
spectral contributions to $\sigma^{(2)}_{\mathrm{ens}}$ in Fig.3b as in the
single-molecule model discussed in Fig.(2)b. In addition, both the full IR
spectrum and the bare cavity response are very similar with a minor intensity
mismatch for the upper polariton peak (cf. Fig.3c). We like to note here that
the explicit ensemble hosts $M-1$ “dark states” with vanishing intensity at
the bare molecular frequency, $\hbar\omega_{\mathrm{as}}=\hbar\omega_{c}$,
which resemble purely molecular linear combinations of all CO2 vibrations
participating in VSC and are to be contrasted by uncoupled molecular states
with non-vanishing intensity in the IR spectrum of Fe(CO)5 under VSC. From the
ensemble analysis, we conclude that transmission spectra obtained from CBO-PT
linear response theory can be calculated via effectively scaled single-
molecule approaches subject to an effective light-matter interaction constant,
$g_{\mathrm{eff}}=\sqrt{M}\,g_{0}$, when inter-molecular interactions are
neglected in the ensemble in line with the non-perturbative CBO approach[15].

Conclusions.— We introduced a cavity Born-Oppenheimer perturbation theory
(CBO-PT) linear response approach to calculate electron-photon correlation
corrected vibro-polaritonic IR spectra of polyatomic molecules under VSC. This
approach approximates non-perturbative CBO linear response theory[15] by
perturbatively approximating nonresonant interactions between electrons and
low-frequency cavity modes, while fully relying on molecular ab initio quantum
chemistry methods. In zeroth-order, we obtain the bare molecular spectrum and
in first-order the linear vibro-polaritonic IR spectrum in absence of
electron-photon correlation. Electron-photon correlation is accounted for at
second-order of CBO-PT linear response theory, where it manifests as static
polarizability-dependent cavity intensity component and related corrections of
Hessian matrix elements, most notably a matter-mediated inter-cavity mode
coupling. Those corrections capture the main characteristics of non-
perturbative CBO linear response theory[15] and provide additional physical
insight due to their explicit relation to molecular and cavity mode
properties. A comparison of uncorrelated CBO-PT(1) and correlation corrected
CBO-PT(2) IR spectra for CO2 and Fe(CO)5 vibro-polaritonic models reveals the
impact of electron-photon correlation on vibro-polaritonic intensity ratios,
Rabi splittings and experimentally relevant transmission spectra, which probe
only the cavity response of the light-matter hybrid system. In both model
scenarios, CBO-PT(2) linear response theory exhibits significant qualitative
agreement with non-perturative CBO linear response results reported in
Ref.[15]. In summary, CBO-PT linear response theory constitutes a promising
since easily accessible approach to vibro-polaritonic IR spectra of polyatomic
molecular systems under VSC accounting for non-trivial electron-photon
correlation effects. A systematic quantitative comparison of CBO-PT and non-
perturbative CBO linear response theories based on controllable ab initio wave
function approaches[16] will be subject of future work.

## Acknowledgements

The authors acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy within
“Sonderforschungsbereich 1636” of University Potsdam, project A05,
“Understanding and controlling reactivity under vibrational and electronic
strong coupling”. E.W. Fischer acknowledges the kind hospitality of Michael
Roemelt and his group at Humboldt-Universität zu Berlin, and helpful
discussions with Thomas Schnappinger (Stockholm). J.A. Syska acknowledges
support and funding by the International Max Planck Research School for
Elementary Processes in Physical Chemistry.

## Data Availability Statement

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

## Conflict of Interest

The authors have no conflicts to disclose.

## References

  * [1] Ebbesen, T.W. Hybrid light-matter states in a molecular and material science perspective. Acc. Chem. Res. 49, 2403 (2016).
  * [2] Hirai, K.; Hutchison, J.A.; Uji-i, H. Recent progress in vibropolaritonic chemistry. ChemPlusChem 85, 1981 (2020).
  * [3] Nagarajan, K.; Thomas, A.; Ebbesen, T.W. Chemistry under vibrational strong coupling. J. Am. Chem. Soc. 143, 16877 (2021).
  * [4] Shalabney, A.; George, J.; Hutchison, J.A.; Pupillo, G.; Genet, C.; Ebbesen, T.W. Coherent coupling of molecular resonators with a microcavity mode. Nat. Commun. 6, 5981 (2015).
  * [5] Shalabney, A.; George, J.; Hiura, H., Hutchison, J.A.; Genet, C.; Hellwig, P.; Ebbesen, T.W. Enhanced Raman Scattering from Vibro-Polariton Hybrid States. Angew. Chem., Int. Ed. 54, 7971-7975 (2015).
  * [6] George, J.; Shalabney, A.; Hutchison, J.A.; Genet, C.; Ebbesen, T.W. Liquid-phase vibrational strong coupling. J. Phys. Chem. Lett. 6, 1027 (2015).
  * [7] Long, J.P.; Simpkins, B.S. Coherent coupling between a molecular vibration and Fabry-Perot optical cavity to give hybridized states in the strong coupling limit. ACS Photonics 2, 130 (2015).
  * [8] George, J.; Chervy, T.; Shalabney, A.; Devaux, E.; Hiura, H.; Genet, C.; Ebbesen, T.W. Multiple Rabi Splittings under Ultrastrong Vibrational Coupling. Phys. Rev. Lett. 117, 153601 (2016).
  * [9] Chervy, T.; Thomas, A.; Akiki, E.; Vergauwe, R.M.A.; Shalabney, A.; George, J.; Devaux, E.; Hutchison, J.A.; Genet, C.; Ebbesen, T.W. Vibro-polaritonic IR emission in the strong coupling regime. ACS Photonics 5, 217 (2018).
  * [10] Xiang, B.; Ribeiro, R.F.; Dunkelberger, A.D.; Wang, J.; Li, Y.; Simpkins, B.S.; Owrutsky, J.C.; Yuen-Zhou, J.; Xiong, W. Two-dimensional infrared spectroscopy of vibrational polaritons. Proc. Natl. Acad. Sci. 115, 4845 (2018).
  * [11] Li, T. E.; Subotnik, J. E.; Nitzan, A.; Cavity molecular dynamics simulations of liquid water under vibrational ultrastrong coupling. Proc. Natl. Acad. Sci. 117, 18324 (2020).
  * [12] Fischer, E.W.; Saalfrank, P. Ground state properties and infrared spectra of anharmonic vibrational polaritons of small molecules in cavities. J. Chem. Phys. 154, 104311 (2021).
  * [13] Lieberherr, A. Z.; Furniss, S. T. E.; Lawrence, J. E.; Manolopoulos, D. E.; Vibrational strong coupling in liquid water from cavity molecular dynamics. J. Chem. Phys. 158, 234106 (2023).
  * [14] Gómez, J.A.; Vendrell, O. Vibrational Energy Redistribution and Polaritonic Fermi Resonances in the Strong Coupling Regime. J. Phys. Chem. A 127, 1598 (2023).
  * [15] Bonini, J.; Flick, J. Ab initio linear-response approach to vibro-polaritons in the cavity Born-Oppenheimer approximation. J. Chem. Theory Comput. 18, 2764 (2022).
  * [16] Schnappinger, T.; Kowalewski, M. Ab-Initio Vibro-Polaritonic Spectra in Strongly Coupled Cavity-Molecule Systems. J. Chem. Theory Comput. 19, 9278 (2023).
  * [17] Flick, J.; Ruggenthaler, M.; Appel, H.; Rubio, A. Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry. PNAS 114, 3026 (2017).
  * [18] Flick, J.; Appel, H.; Ruggenthaler, M.; Rubio, A. Cavity Born-Oppenheimer approximation for correlated electron-nuclear-photon systems. J. Chem. Theory Comput. 13, 1616, (2017).
  * [19] Fischer, E.W.; Saalfrank, P. Beyond Cavity Born-Oppenheimer: On Nonadiabatic Coupling and Effective Ground State Hamiltonians in Vibro-Polaritonic Chemistry. J. Chem. Theory Comput. 19, 7215-7229 (2023).
  * [20] Galego, J.; Climent, C; Garcia-Vidal, F.J.; Feist, J. Cavity Casimir-Polder forces and their effects in ground-state chemical reactivity.  Phys. Rev. X 9, 021057 (2019).
  * [21] Szidarovszky, T. An efficient and flexible approach for computing rovibrational polaritons from first principles. arXiv preprint arXiv:2304.02346 (2023).
  * [22] Haugland, T. S.; Philbin, J. P.; Ghosh, T. K.; Chen, M.; Koch, H.; Narang, P. Understanding the polaritonic ground state in cavity quantum electrodynamics. arXiv preprint arXiv:2307.14822 (2023).
  * [23] Schäfer, C.; Fojt, J.; Lindgren, E., Erhart, P. Machine Learning for Polaritonic Chemistry: Accessing chemical kinetics. arXiv preprint arXiv:2311.09739 (2023).
  * [24] The interaction parameter $g_{0}$ is equivalent to $\lambda_{c}$ employed by Bonini and Flick.
  * [25] Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 91 (2003).
  * [26] Staroverov, V. N.; Scuseria, G. E.; Tao J.; Perdew, J. P. Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes. J. Chem. Phys. 119, 12129 (2003).
  * [27] Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 98, 5648, (1993).
  * [28] Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * [29] Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 8 1057 (2006).
  * [30] Ulusoy, I.S.; Vendrell, O. Dynamics and spectroscopy of molecular ensembles in a lossy microcavity. J. Chem. Phys. 153, 044108 (2020).
  * [31] Fischer, E.W.; Saalfrank, P. Cavity-induced non-adiabatic dynamics and spectroscopy of molecular rovibrational polaritons studied by multi-mode quantum models. J. Chem. Phys. 157, 034305 (2022).
  * [32] Chen, T.-T.; Du, M.; Yang, Z.; Yuen-Zhou, J.; Xiong, W. Cavity-enabled enhancement of ultrafast intramolecular vibrational redistribution over pseudorotation. Science 378, 794 (2022).
  * [33] Sidler, D.; Schnappinger, T.; Obzhirov, A.; Ruggenthaler, M.; Kowalewski, M.; Rubio, A. Unraveling a Cavity Induced Molecular Polarization Mechanism from Collective Vibrational Strong Coupling. arXiv preprint arXiv:2306.06004 (2023).