††thanks: Denotes equal contribution

# Molecular van der Waals fluids in cavity quantum electrodynamics

John P. Philbin¶<EMAIL_ADDRESS>Harvard John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA
College of Letters and Science, University of California, Los Angeles, CA
90095, USA  Tor S. Haugland¶ Department of Chemistry, Norwegian University of
Science and Technology, 7491 Trondheim, Norway  Tushar K. Ghosh¶ Department of
Chemistry, Purdue University, West Lafayette, IN 47907, USA  Enrico Ronca
Dipartimento di Chimica, Biologia e Biotecnologie, Università degli Studi di
Perugia, Via Elce di Sotto, 8, 06123, Perugia, Italy  Max Planck Institute for
the Structure and Dynamics of Matter and Center Free-Electron Laser Science,
Luruper Chaussee 149, 22761 Hamburg, Germany  Ming Chen<EMAIL_ADDRESS>Department of Chemistry, Purdue University, West Lafayette, IN 47907, USA
Prineha Narang<EMAIL_ADDRESS>Harvard John A. Paulson School of Engineering
and Applied Sciences, Harvard University, Cambridge, MA 02138, USA  College of
Letters and Science, University of California, Los Angeles, CA 90095, USA
Henrik Koch<EMAIL_ADDRESS>Department of Chemistry, Norwegian University
of Science and Technology, 7491 Trondheim, Norway  Scuola Normale Superiore,
Piazza dei Cavalieri, 7, 56124 Pisa, Italy

###### Abstract

Intermolecular van der Waals interactions are central to chemical and physical
phenomena ranging from biomolecule binding to soft-matter phase transitions.
However, there are currently very limited approaches to manipulate van der
Waals interactions. In this work, we demonstrate that strong light-matter
coupling can be used to tune van der Waals interactions, and, thus, control
the thermodynamic properties of many-molecule systems. Our analyses reveal
orientation dependent single molecule energies and interaction energies for
van der Waals molecules (for example, H2). For example, we find intermolecular
interactions that depend on the distance between the molecules $R$ as $R^{-3}$
and $R^{0}$. Moreover, we employ non-perturbative ab initio cavity quantum
electrodynamics calculations to develop machine learning-based interaction
potentials for molecules inside optical cavities. By simulating systems
ranging from $12$ H2 to $144$ H2 molecules, we demonstrate that strong light-
matter coupling can tune the structural and thermodynamic properties of
molecular fluids. In particular, we observe varying degrees of orientational
order as a consequence of cavity-modified interactions, and we explain how
quantum nuclear effects, light-matter coupling strengths, number of cavity
modes, molecular anisotropies, and system size all impact the extent of
orientational order. These simulations and analyses demonstrate both local and
collective effects induced by strong light-matter coupling and open new paths
for controlling the properties of molecular clusters.

††preprint: APS/123-QED††preprint: APS/123-QED

Van der Waals interactions are ubiquitous in chemistry and physics, playing
important roles in diverse scientific fields ranging from DNA base stacking to
2D material interlayer interactions.Hobza and Šponer (2002); Novoselov  _et
al._ (2016); Sternbach  _et al._ (2021) There has been a long history of
attempting to elucidate the origin of van der Waals interactions;Maitland  _et
al._ (1981); Stone (2013) the first quantum mechanical derivation was
performed by London in the 1930s using second-order perturbation theory.London
(1937) London found that two molecules that do not have permanent dipoles
(e.g. H2), which we refer to as van der Waals molecules, have an attractive
interaction between them that scales with the distance between the molecules
$R$ as $R^{-6}$.London (1937) This $R^{-6}$ attractive force is commonly used
as the long-distance asymptotic form of van der Waals interactions in many
force fields and to correct van der Waals interactions in ab initio
calculations, which have both achieved great successes in modeling
thermodynamic properties in a variety of systems.Halgren (1992); Grimme  _et
al._ (2010) Despite van der Waals interactions being central to many
properties of molecular and condensed matter systems, limited approaches have
been proposed to manipulate intermolecular van der Waals interactions.
However, applied electromagnetic fields have been shown to modify van der
Waals interactions between atoms and molecules,Thirunamachandran (1980);
Milonni and Smith (1996); Sherkunov (2009); Fiscelli  _et al._ (2020) and
Haugland et al.Haugland  _et al._ (2021) recently showed numerically that van
der Waals interactions are significantly altered by strong light-matter
coupling in optical cavities. These studies open the possibility of
controlling the properties and structure of molecular fluids by tuning the
light-matter coupling parameters, the coupling strength and frequency.

Figure 1: (A) Schematic representation of the findings from our simulations of
a fluid of H2 molecules outside and inside a cavity. Specifically,
orientational order can be observed inside a cavity whereas the H2 molecules
can rotate freely outside of a cavity. The dashed lines represent the
different intermolecular interaction length scales outside and inside a
cavity. (B) Diagram describing the computational workflow used in this work.
Ab initio cavity QED energies and corresponding symmetry preserving features
(see Fig. S3, Table S3 and Section SIV.A.1 for details of symmetry preserving
features) of many $2$H2 configurations are used to develop neural network-
based intermolecular pair potentials capable of being utilized in path
integral molecular dynamics simulations of fluids of H2 molecules.

The goal of this work is to understand how the structure of molecular van der
Waals fluids can be modulated using enhanced vacuum electromagnetic quantum
fields, and we focus on the impact that a single strongly coupled photon mode
can have on the properties of a model molecular van der Waals fluids. To this
end, we leverage recent developments in cavity quantum electrodynamics (QED)
simulations and neural network pair potentials to simulate molecular fluids of
H2 molecules strongly coupled to a single photon mode (Fig. 1). By analyzing
how cavity-modified single molecule energies and cavity-mediated
intermolecular interactions depend on the orientation of the H2 molecules both
relative to the cavity polarization vector and relative to one another, we can
explain how cavities impact the structure and orientational order of molecular
van der Waals fluids. The findings reported herein should readily be
transferable to other molecules and light-matter regimes (e.g. vibrational
polaritons) given the generality of the cavity QED Hamiltonian used in this
work.Ribeiro  _et al._ (2018); Rivera  _et al._ (2019); Thomas  _et al._
(2019); Li  _et al._ (2020); Garcia-Vidal  _et al._ (2021); Li  _et al._
(2021a) We also discuss how the light-matter coupling strength, number of
cavity modes, temperature, anisotropic polarizabilities of molecules, quantum
nuclear effects, and molecular concentrations can all impact the extent of
orientational order observed in any particular cavity QED experiment.Vahala
(2003); Cortese  _et al._ (2017); Joseph  _et al._ (2021); Fukushima  _et al._
(2022); Sandeep  _et al._ (2022)

In molecular dynamics (MD) simulations, the nuclei move along electronic
potential energy surfaces. In the cavity case, where the photon contributions
are added, these surfaces have been termed polaritonic potential energy
surfaces.Galego  _et al._ (2015); Lacombe  _et al._ (2019); Fregoni  _et al._
(2022) In both cases, the total potential energy of $N$ H2 molecules can be
calculated as a many-body expansion,

$E_{\text{total}}=\sum_{A}E_{A}+\sum_{\left\langle
A,B\right\rangle}E_{AB}+{\sum_{\left\langle A,B,C\right\rangle}E_{ABC}}+...,$
(1)

where $E_{A}$ represents the single-molecule energies, $E_{AB}$ represents the
intermolecular interaction energies between all unique pairs of molecules, and
so on for higher-body terms. In this work, we focus on contributions to the
total energy in Eq. 1 arising from at most two-body interactions. The three-
body and higher-body terms are significantly smaller than the two-body
interactions per interaction, see the Supplementary Information (SI) for
details. Outside the cavity, the one-body term does not depend on the
orientation of the H2 molecule. On the other hand, inside the cavity, the
molecule-field interaction causes the one-body energies to depend on the
orientation of the H2 molecules with respect to the optical cavity
polarization vector, $\bm{\varepsilon}$. Furthermore, the two-body energies
depends on the orientation between the two molecules as well as their
orientation relative to the field as a consequence of the anisotropic
polarizability of H2 molecules, in contrast to isotropic polarizabilities of
atoms.Thirunamachandran (1980); Milonni and Smith (1996); Sherkunov (2009);
Fiscelli  _et al._ (2020)

Figure 2: (A-B) Snapshots taken at thermal equilibrium from molecular dynamic
(MD) simulations in the case of (A) no cavity (orange) and (B) cavity-modified
one-body and two-body terms (blue). (C) The impact of quantum nuclear effects
are demonstrated by comparing the molecular bond axis to cavity polarization
vector ($\theta_{A\varepsilon}$), angular probability distribution function,
$P\left(\theta_{A\varepsilon}\right)$ for path integral molecular dynamics
(PIMD) simulations of H2, D2, T2, and a classical MD simulation of H2. (D)
Molecular bond axis of molecule $A$ to molecular bond axis of molecule $B$
($\theta_{AB}$) angular probability distribution function,
$P\left(\theta_{AB}\right)$ and (E) $P\left(\theta_{A\varepsilon}\right)$ are
shown for PIMD simulations for no cavity (orange), cavity (blue), and cavity-
modified one-body term but no cavity two-body term (green) cases. (F)
$P\left(\theta_{A\varepsilon}\right)$ is shown for two different PIMD
simulations containing different numbers of H2 molecules within the same
cavity volume (i.e. changing the molecular density). All PIMD simulations
shown in this figure were performed using neural networks trained with CCSD
(no cavity) or QED-CCSD-12-SD1 with $\lambda=0.1$ a.u. (cavity) calculated
energies. All entropic contributions to angle distribution functions are
removed.

We calculate $E_{A}$ and $E_{AB}$ by solving the Schrödinger equation for the
cavity QED Hamiltonian in the dipole approximation with a single photon mode
using accurate coupled cluster (QED-CCSD-12-SD1) and near exact full
configuration interaction (QED-FCI-5).Haugland  _et al._ (2020) Our single
photon mode has a coupling constant of $\lambda=0.1$ a.u. and energy of
$\hbar\omega_{c}=13.6$ eV unless specified otherwise. This coupling constant
is rather large as it corresponds to the coupling of at least $5$ independent
modes where each has an effective volume of $0.9$ nm3. We detail below how the
cavity-modified local interactions and cavity-induced collective effects
depend on $\lambda$. More than $100,000$ H2 dimer configurations are used as
inputs to a fully-connected neural network that serves as our intermolecular
pair potential, which is trained and tested against the calculated energies.
The trained potential energy functions were carefully tested, and, in the SI,
we demonstrate that our machine learning models are fully capable of
reproducing the potential energy surfaces. In Fig. 1B, we show the
computational workflow used in this work schematically. In this study, we
focus on path integral molecular dynamics (PIMD) simulations in order to
account for quantum nuclear effects. Our PIMD simulations of fluids of H2
molecules were performed with a fixed number of molecules ($N$), temperature
($T$), and volume ($V$). All PIMD simulations presented herein were performed
with a molecular density of $13$ molecules per nm3, temperature of $70$ K, and
$N=12$ unless otherwise specified. More details on the simulations, including
comparisons of QED-CCSD-12-SD1 with QED-FCI-5, comparisons of MD with PIMD,
and additional parameter regimes (e.g. smaller $\lambda$ values), are provided
in the SI.

Figure 3: (A) Energy difference, $\Delta E$, between a single H2 molecule
inside a cavity aligned perfectly along the cavity polarization vector,
$\bm{\varepsilon}$, and different angles relative to the cavity polarization
vector. The inset shows the energy of a single molecule within a cavity
increases with $\lambda^{2}$. (B) Intermolecular interaction energies,
$E_{AB}$, and fits to a Lennard-Jones type potential given by Eq. 3 (dashed
lines) and cavity-modified Lennard Jones type potential given by Eq. 4 (solid
line). (C) Intermolecular interaction energies, $E_{AB}$, at $25$ Å for
various high symmetry molecular orientations and cavity polarizations. All
calculations shown in this figure were performed using QED-CCSD-12-SD1 with
$\lambda=0.1$ a.u.

The structural properties of the molecular van der Waals fluids are analyzed
using PIMD simulation trajectories. In Fig. Fig. S13, we summarize the main
findings of our PIMD and classical MD simulations. Fig. Fig. S13A and Fig.
Fig. S13B show representative thermal equilibrium configurations for the no
cavity (orange) and cavity (blue) scenarios, respectively. The impact of the
cavity-modified interactions are observable in the orientational order of the
H2 molecules both relative to the cavity polarization vector
($\theta_{A\varepsilon}$, Figs. Fig. S13C, E and F) and relative to other H2
molecules ($\theta_{AB}$, Fig. Fig. S13D). Specifically, Figs. Fig. S13C-F
show that the cavity-modified energies enhance the probability of finding two
molecules oriented parallel to one another (i.e. $\theta_{AB}=0,\pi$) and
perpendicular to the cavity polarization vector (i.e.
$\theta_{A\varepsilon}=\frac{\pi}{2}$). However, the extent of this
orientational order depends on many factors, including the magnitude of
quantum nuclear effects, the light-matter coupling strengths, molecular
anisotropies, and number of molecules. To elucidate the importance of quantum
nuclear effects, we compare the orientational order observed in PIMD
simulations of H2, D2, and T2 with a classical MD simulation of H2 in Fig.
Fig. S13C; the degree of orientational order monotonically increases upon
increasing the molecular masses from H2 to D2 to T2 (which reduces quantum
nuclear effects) and is further enhanced when quantum nuclear effects are
completely removed as in the classical MD simulation. Next, in Figs. Fig.
S13D-F, we show how cavity-modified one-body energies and two-body
intermolecular energies each impact the orientational order. Fig. Fig. S13D
and Fig. Fig. S13E demonstrate that the cavity-modified one-body energies are
the dominant driver of the orientational order for the case of $12$ H2
molecules. The orange lines in Figs. Fig. S13D,E show that the H2 molecules
have no preferred orientation axis outside the cavity, consistent with the
global rotational symmetry of the electronic and nuclear Hamiltonian in
absence of the cavity. However, the presence of the bilinear coupling and
dipole self-energy terms break this symmetry such that H2 molecules prefer to
orient their bond axis in specific orientations relative to the cavity
polarization vector and relative to one another. In particular, the dipole
self-energy term outcompetes the bilinear coupling term and is responsible for
the $12$ molecule simulations preferentially aligning perpendicular to the
cavity polarization vector (Fig. 3A). However, Figs. Fig. S13E,F demonstrate
that the cavity-modified one-body energies lead to this perpendicular
alignment whereas the cavity-modified two-body intermolecular interactions
attempt to align the molecules parallel to the cavity polarization vector.
Specifically, the green line in Fig. Fig. S13E shows that the cavity-modified
one-body term causes H2 molecules to preferentially align perpendicular to the
cavity polarization vector (i.e. $\theta_{A\varepsilon}=\frac{\pi}{2}$), and
the inclusion of cavity-modified two-body interactions begins to counteract
this effect as seen in the blue line in Fig. Fig. S13E reducing the
orientational alignment. This effect of the two-body interactions causing the
H2 molecules to preferentially align parallel to the cavity polarization
vector (i.e. $\theta_{A\varepsilon}=0,\pi$) and the collective nature of the
cavity-modified intermolecular interactions are highlighted in Fig. Fig. S13F
and Fig. S13. We find that for a small number of molecules (e.g. $N=12$) the
one-body term dominates and the molecules preferentially align perpendicular
to the cavity polarization vector, but as $N$ increases to $144$ H2 molecules
with a fixed coupling and cavity volume the orientational order is lost due
the cavity-modified one-body and two-body effects perfectly canceling one
another. Additionally, the extent of orientational order induced by the cavity
decreases as the light-matter coupling strength decreases as shown in Fig. S8
and explained analytically below.

Although we performed non-perturbative ab initio cavity QED calculations,
perturbation theory can be used to further analyze and explain the major
findings of our PIMD and MD simulations. We summarize our key findings here
and in Fig. 3, and the complete analysis is provided in the SI. The cavity
modifications to the one-body energies, $E_{A}$, results in the H2 molecules
aligning their bonds orthogonal to the cavity polarization. This occurs
because H2 is most polarizable along its bond axis, and, from perturbation
theory, we can obtain an expression for the cavity-modified one-body energy as

$E_{A}^{\text{cavity}}\approx E_{A}^{\text{no
cavity}}+c\,(\alpha_{\parallel}\cos^{2}{\theta_{A\varepsilon}}+\alpha_{\perp}\sin^{2}{\theta_{A\varepsilon}}),$
(2)

where $\alpha_{\parallel}$ and $\alpha_{\perp}$ are the polarizabilities of
molecular hydrogen along its bond axis and perpendicular axes, respectively,
and $c$ is a positive scalar constant proportional to the molecule-cavity
coupling squared (i.e. $c\propto\lambda^{2}$). Eq. 2 is in agreement with the
ab initio calculations shown in Fig. 3A. Interestingly, the dipole self-energy
term increases the energy of a single molecule in a cavity more than the
bilinear coupling term decreases the energy (Eq. S12); thus, the lowest energy
orientation of a single molecule in a cavity is such that its most polarizable
axis is perpendicular to the cavity polarization vector (or vectors in terms
of multimode cavities).

In terms of the cavity modifications to the two-body energies, Fig. 3B shows
the intermolecular interaction between two H2 molecules as a function of the
center-to-center distance ($R$). The impact of the cavity on this dissociation
curve at first glance appears modest, even for the rather large light-matter
coupling of $\lambda=0.1$ a.u., but these modifications can impact the
structural and thermodynamic properties of molecular van der Waals systems for
a few reasons. First, a standard intermolecular van der Waals interaction
potential given by

$E_{AB}^{\text{no cavity}}=\frac{c_{6}}{R^{6}}+E_{\text{short-range}},$ (3)

where $E_{\text{short-range}}$ accounts for the short-range repulsion between
van der Waals molecules and the $R^{-6}$ term is the usual attractive London
dispersion interaction, is not applicable inside an optical cavity (Fig.
3B).Thirunamachandran (1980); Milonni and Smith (1996); Sherkunov (2009);
Fiscelli  _et al._ (2020) A modified interaction potential that includes
angle-dependent terms that scale as $R^{-3}$ and $R^{0}$ is necessary inside
an optical cavity such that the interaction between two van der Waals
molecules is given by

$E_{AB}^{\text{cavity}}=\frac{c_{0}}{R^{0}}+\frac{c_{3}}{R^{3}}+\frac{c_{6}}{R^{6}}+E_{\text{short-
range}}.$ (4)

These interactions arise as early as second-order perturbation theory (see SI
Eq. S9).Thirunamachandran (1980) The $R^{0}$ interaction between a single pair
of molecules is rather weak ($c_{0}\propto\lambda^{4}$) as shown in Fig. 3C.
However, due to its long-range nature, a single molecule interacts with all
other molecules, and, thus, the collective effect of this interaction can
become large in many-molecule simulations. Importantly, this interaction
strength depends on the orientations of both molecular bonds relative to the
cavity polarization (Fig. 3C). Specifically, the interaction energy is
minimized when the molecular bonds of both molecules are parallel to the
cavity polarization vector, because the interaction strength of this term is
approximately related to the product of the polarizability of each molecule
along $\bm{\varepsilon}$
($c_{0}\propto\alpha_{A\varepsilon}\alpha_{B\varepsilon}$). And because
$c_{0}$ is always negative, this $R^{0}$ intermolecular interaction increases
the probability of finding H2 molecules parallel to the cavity polarization
vector and decreases the probability to find the molecules perpendicular to
the polarization vector (Fig. Fig. S13E,F). The collective nature of this
interaction is demonstrated in Fig. Fig. S13F and Fig. S13 where the
orientational order depends on the number of H2 molecules for simulations with
the same simulation volume but different molecular densities. At $N=144$, the
orientational order due to the two-body interactions have become so large that
they entirely cancel out the orientational effects from the cavity modified
one-body energies that are dominated by dipole self-energy effects for $N=12$
molecules. As $N$ increases further, we expect that the system will completely
flip, and instead align parallel to the polarization vector. This is
demonstrated in the SI (Fig. S13), but the number of molecules required
($N\geq 1000$) is too large to justify in a realistic system with the coupling
we are using currently. Both the cavity-modified $R^{-6}$ and cavity-induced
$R^{-3}$ interactions scale with $\lambda^{2}$ at lowest order. Importantly,
the $R^{-3}$ interaction is not a result of the cavity inducing a dipole
moment in the H2 molecules but rather an interaction taking place via the
cavity mode. As discussed in the SI in more detail, the intermolecular angle
and molecule-cavity angle dependencies of the perturbation potential combine
to create the orientational order shown throughout Fig. Fig. S13.

In summary, we have demonstrated that strong light-matter coupling to a single
photon mode can have profound impacts on the properties of molecular van der
Waals fluids by combining ab initio cavity QED calculations with path integral
molecular dynamics simulations of many H2 molecules. We found that cavity-
modified single molecule and intermolecular interaction energies result in
significantly changed molecular orientational order, even in the fluid phase.
We look forward to seeing future experimental and theoretical studies that aim
to elucidate how processes such as ion and molecular diffusion, intermolecular
energy transfer,Zhong  _et al._ (2016); Du  _et al._ (2018); Xiang  _et al._
(2020) and chemical reactivityHerrera and Spano (2016); Thomas  _et al._
(2019); Yang and Cao (2021); Li  _et al._ (2021b); Simpkins  _et al._ (2021);
Philbin  _et al._ (2022) are impacted by the unique properties of molecular
fluids in cavity QED reported here.

###### Acknowledgements.

We thank Jonathan Curtis, Davis Welakuh, Wenjie Dou, and Rosario R. Riso for
helpful discussions. This work was primarily supported by the Department of
Energy, Photonics at Thermodynamic Limits Energy Frontier Research Center,
under Grant No. DE-SC0019140 and European Research Council under the European
Union’s Horizon 2020 Research and Innovation Programme grant agreement No.
101020016. An award of computer time was provided by the INCITE program. This
research also used resources of the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725 J.P.P. also acknowledges support from the Harvard University
Center for the Environment. T.K.G. and M.C. acknowledge support from Purdue
startup funding. T.S.H. and H.K. also acknowledges funding from the Research
Council of Norway through FRINATEK project 275506. P.N. acknowledges support
as a Moore Inventor Fellow through Grant No. GBMF8048 and gratefully
acknowledges support from the Gordon and Betty Moore Foundation as well as
support from a NSF CAREER Award under Grant No. NSF-ECCS-1944085. E.R
acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon Europe Research and Innovation Programme (Grant n.
ERC-StG-2021-101040197 - QED-SPIN).

## References

  * Hobza and Šponer (2002) P. Hobza and J. Šponer, J. Am. Chem. Soc. 124, 11802 (2002).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, aac9439 (2016).
  * Sternbach  _et al._ (2021) A. J. Sternbach, S. H. Chae, S. Latini, A. A. Rikhter, Y. Shao, B. Li, D. Rhodes, B. Kim, P. J. Schuck, X. Xu, X. Y. Zhu, R. D. Averitt, J. Hone, M. M. Fogler, A. Rubio, and D. N. Basov, Science 371, 617 (2021).
  * Maitland  _et al._ (1981) G. C. Maitland, G. D. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham,  _Intermolecular Forces: Their Origin and Determination_ (Oxford University Press, USA, 1981).
  * Stone (2013) A. Stone,  _The Theory of Intermolecular Forces_, 2nd ed. (Oxford University Press, Oxford, 2013) p. 352.
  * London (1937) F. London, Trans. Faraday Soc. 33, 8b (1937).
  * Halgren (1992) T. A. Halgren, J. Am. Chem. Soc. 114, 7827 (1992).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Thirunamachandran (1980) T. Thirunamachandran, Mol. Phys. 40, 393 (1980).
  * Milonni and Smith (1996) P. W. Milonni and A. Smith, Phys. Rev. A 53, 3484 (1996).
  * Sherkunov (2009) Y. Sherkunov, J. Phys. Conf. Ser. 161, 012041 (2009).
  * Fiscelli  _et al._ (2020) G. Fiscelli, L. Rizzuto, and R. Passante, Phys. Rev. Lett. 124, 013604 (2020).
  * Haugland  _et al._ (2021) T. S. Haugland, C. Schäfer, E. Ronca, A. Rubio, and H. Koch, J. Chem. Phys. 154, 094113 (2021).
  * Ribeiro  _et al._ (2018) R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem. Sci. 9, 6325 (2018).
  * Rivera  _et al._ (2019) N. Rivera, J. Flick, and P. Narang, Phys. Rev. Lett. 122, 193603 (2019).
  * Thomas  _et al._ (2019) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Science 363, 615 (2019).
  * Li  _et al._ (2020) T. E. Li, J. E. Subotnik, and A. Nitzan, Proc. Natl. Acad. Sci. U. S. A. 117, 18324 (2020).
  * Garcia-Vidal  _et al._ (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, eabd0336 (2021).
  * Li  _et al._ (2021a) T. E. Li, A. Nitzan, and J. E. Subotnik, Angew. Chemie 133, 15661 (2021a).
  * Vahala (2003) K. J. Vahala, Nature 424, 839 (2003).
  * Cortese  _et al._ (2017) E. Cortese, P. G. Lagoudakis, and S. De Liberato, Phys. Rev. Lett. 119, 043604 (2017).
  * Joseph  _et al._ (2021) K. Joseph, S. Kushida, E. Smarsly, D. Ihiawakrim, A. Thomas, G. L. Paravicini-Bagliani, K. Nagarajan, R. Vergauwe, E. Devaux, O. Ersen, U. H. F. Bunz, and T. W. Ebbesen, Angew. Chem. Int. Ed. 60, 19665 (2021).
  * Fukushima  _et al._ (2022) T. Fukushima, S. Yoshimitsu, and K. Murakoshi, J. Am. Chem. Soc. 144, 12177 (2022).
  * Sandeep  _et al._ (2022) K. Sandeep, K. Joseph, J. Gautier, K. Nagarajan, M. Sujith, K. G. Thomas, and T. W. Ebbesen, J. Phys. Chem. Lett. 13, 1209 (2022).
  * Galego  _et al._ (2015) J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 5, 41022 (2015).
  * Lacombe  _et al._ (2019) L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys. Rev. Lett. 123, 083201 (2019).
  * Fregoni  _et al._ (2022) J. Fregoni, F. J. Garcia-Vidal, and J. Feist, ACS Photonics 9, 1096 (2022).
  * Haugland  _et al._ (2020) T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio, and H. Koch, Phys. Rev. X 10, 041043 (2020).
  * Zhong  _et al._ (2016) X. Zhong, T. Chervy, S. Wang, J. George, A. Thomas, J. A. Hutchison, E. Devaux, C. Genet, and T. W. Ebbesen, Angew. Chem. Int. Ed. 55, 6202 (2016).
  * Du  _et al._ (2018) M. Du, L. A. Martínez-Martínez, R. F. Ribeiro, Z. Hu, V. M. Menon, and J. Yuen-Zhou, Chem. Sci. 9, 6659 (2018).
  * Xiang  _et al._ (2020) B. Xiang, R. F. Ribeiro, M. Du, L. Chen, Z. Yang, J. Wang, J. Yuen-Zhou, and W. Xiong, Science 368, 665 (2020).
  * Herrera and Spano (2016) F. Herrera and F. C. Spano, Phys. Rev. Lett. 116, 238301 (2016).
  * Yang and Cao (2021) P. Y. Yang and J. Cao, J. Phys. Chem. Lett. 12, 9531 (2021).
  * Li  _et al._ (2021b) X. Li, A. Mandal, and P. Huo, Nat. Commun. 12, 1315 (2021b).
  * Simpkins  _et al._ (2021) B. S. Simpkins, A. D. Dunkelberger, and J. C. Owrutsky, J. Phys. Chem. C 125, 19081 (2021).
  * Philbin  _et al._ (2022) J. P. Philbin, Y. Wang, P. Narang, and W. Dou, J. Phys. Chem. C 126, 14908 (2022).
  * White  _et al._ (2020) A. F. White, Y. Gao, A. J. Minnich, and G. K. L. Chan, J. Chem. Phys. 153, 224112 (2020).
  * Eisenschitz and London (1930) R. Eisenschitz and F. London, Zeitschrift für Phys. 60, 491 (1930).
  * London (1930) F. London, Zeitschrift für Phys. 63, 245 (1930).
  * Dahlke and Truhlar (2007) E. E. Dahlke and D. G. Truhlar, J. Chem. Theory Comput. 3, 46 (2007).
  * Barron (2017) J. T. Barron, Continuously differentiable exponential linear units (2017), arXiv:1704.07483 .
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014), arXiv:1412.6980 .
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Bussi and Parrinello (2007) G. Bussi and M. Parrinello, Phys. Rev. E 75, 056707 (2007).
  * Ceriotti  _et al._ (2009) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 103, 030603 (2009).
  * Ceriotti  _et al._ (2010a) M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010a).
  * Ceriotti  _et al._ (2011) M. Ceriotti, D. E. Manolopoulos, and M. Parrinello, J. Chem. Phys. 134, 084104 (2011).
  * Ceriotti  _et al._ (2010b) M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 133, 124104 (2010b).
  * Ceriotti  _et al._ (2014) M. Ceriotti, J. More, and D. E. Manolopoulos, Comput. Phys. Commun. 185, 1019 (2014).

Supplementary Information:
Molecular van der Waals fluids in cavity quantum electrodynamics

###### Contents

  1. I Ab Initio Calculations
  2. II Perturbation Theory
  3. III Many-body Interactions
  4. IV Molecular Dynamics
    1. IV.1 Training Potential Energy Functions for Simulating Fluids of H2
      1. IV.1.1 Neural Network-based Pairwise Interactions
      2. IV.1.2 Single Molecule Potential Energies
    2. IV.2 Molecular Dynamics
      1. IV.2.1 Classical Molecular Dynamics
      2. IV.2.2 Path Integral Molecular Dynamics
    3. IV.3 Radial Distribution Functions
    4. IV.4 Angular Distribution Functions
  5. V Additional Results
    1. V.1 Comparison of Radial Distribution Functions
    2. V.2 Comparison of Classical MD and PIMD
    3. V.3 Comparison of QED-FCI-5 and QED-CCSD-12-SD1
    4. V.4 $\lambda$ Dependent Molecular Alignment

##  I Ab Initio Calculations

The Hamiltonian used in the ab initio calculations is the single mode Pauli-
Fierz Hamiltonian in the length gauge

$\displaystyle H$
$\displaystyle=H_{e}+\lambda\sqrt{\frac{\omega_{c}}{2}}((\bm{d}-\expectationvalue{\bm{d}})\cdot\bm{\varepsilon})(b+b^{\dagger})$
(S1)
$\displaystyle+\frac{\lambda^{2}}{2}((\bm{d}-\expectationvalue{\bm{d}})\cdot\bm{\varepsilon})^{2}+\omega_{c}b^{\dagger}b,$

where $H_{e}$ is the electronic Hamiltonian, $\lambda$ is the bilinear
coupling, $\omega_{c}$ is the cavity frequency, $\bm{d}$ is the molecular
dipole, $\bm{\varepsilon}$ is the cavity polarization vector, and $b$ and
$b^{\dagger}$ are the photon annihilation and creation operators,
respectively.

All electronic structure calculations are run using an aug-cc-pVDZ basis set.
The optical cavity is described by a single linearly polarized mode coupling
parameter $\lambda$ is set to $0.1$ a.u. and the cavity energy
$\hbar\omega_{c}$ is $13.6$ eV, unless otherwise specified.

The large value for the coupling is partially justified by the single mode
approximation. For cavity-induced changes in the ground state, each cavity
mode will to second order in perturbation theory (see Eq. S12) enter the
energy independently. For larger frequencies, the bilinear contribution from
each mode cancels part of the dipole self-energy. For smaller frequencies
compared to electronic excitation energies, we find that only contributions
from the dipole self-energy are significant. Therefore, in the low-frequency
regime, the coupling from $N_{\rm modes}$ modes is given by an effective
coupling $\lambda^{2}_{\rm eff}\approx N_{\rm modes}\lambda^{2}$.

As shown and discussed in Ref. Haugland  _et al._ (2021), cavity quantum
electrodynamics Hartree-Fock (QED-HF) and current QED density functional
theory (QEDFT) implementations do not describe intermolecular forces properly,
especially van der Waals interactions in which they fail to predict an
attractive interaction between van der Waals molecules. Therefore, we
performed the ab initio simulations with QED coupled cluster (QED-CCSD-12-SD1)
and QED full configuration interaction (QED-FCI).White  _et al._ (2020) QED-
CCSD-12-SD1 is an extension of QED-CCSD-1, as described in Ref. Haugland  _et
al._ (2020), with two-photon excitations. The QED-CCSD-12-SD1 cluster operator
is

$T=T_{1}+T_{2}+S_{1}b^{\dagger}+S_{2}b^{\dagger}+\gamma_{1}b^{\dagger}+\gamma_{2}(b^{\dagger})^{2},$
(S2)

where $T_{1}$ and $T_{2}$ are singles and doubles electron excitations,
$S_{1}b^{\dagger}$ and $S_{2}b^{\dagger}$ are singles and doubles coupled
electron-photon excitations, and $\gamma_{1}b^{\dagger}$ and
$\gamma_{2}(b^{\dagger})^{2}$ are singles and doubles photon excitations. The
reference state is QED-HF as described in Ref. Haugland  _et al._ (2020). QED-
FCI calculations are run with up to five photons (QED-FCI-5) to ensure that
the energy with respect to photon number is converged.

We use QED-CCSD-12-SD1 instead of QED-CCSD-1 (equivalent to QED-CCSD-1-SD1)
because the two-photon excitations are important for properly modeling the
two-body interactions, as tested against QED-FCI-5 calculations. Without two-
photon excitations, the two-body interactions have the wrong sign in the case
of molecules separated by large distances (e.g. molecules separated by more
than $1$ nm). This is visualized in Fig. S1.

Figure Fig. S1:  Calculated intermolecular interaction energies for a C2v
configuration of $2$H2 with the cavity polarization vector parallel to the
center-to-center intermolecular distance vector. All calculations shown in
this figure were performed with $\lambda=0.1$ a.u

In all of our calculations, we use a linearly polarized optical cavity with a
single photon frequency and single polarization vector. In most experiments as
of today, the optical cavity is not limited to just one polarization, but
rather it hosts two degenerate cavity modes with orthogonal polarizations
(both cavity mode polarization vectors are perpendicular to the cavity
wavevector). Since the molecular orientations aligns with the transversal
polarization, we expect that a standard optical cavity, which has both
polarizations, will interact with the system differently. In particular, we
expect that for few molecules, the molecules will orient along the wavevector
$\bm{k}$, perpendicular to both cavity polarization vectors. For many
molecules, we expect that the molecules will align perpendicular to $\bm{k}$,
in the plane defined from the two transversal polarization vectors.

##  II Perturbation Theory

As we demonstrate throughout this work, strong coupling to a single photon
mode fundamentally changes the length scales and orientational dependence in
which van der Waals molecules interact with one another. In this section, we
explain these observations by performing perturbation theory in a similar
spirit as Fritz London did in 1930Eisenschitz and London (1930); London (1930,
1937) but with additional perturbative potentials associated with coupling to
the cavity. This analysis shows cavity-mediated intermolecular interactions
between van der Waals molecules that scale with $R^{-3}$ and distance
independent, $R^{0}$, interactions in addition to modifications to London
dispersion forces that have an $R^{-6}$ dependence.Thirunamachandran (1980);
Milonni and Smith (1996); Sherkunov (2009); Fiscelli  _et al._ (2020)

The total Hamiltonian is given by $H=H^{0}+H^{1}$ with

$H^{0}=H_{e,A}+H_{e,B}+\omega_{c}b^{\dagger}b$ (S3)

where $b^{\dagger}$ and $b$ are photon creation and annihilation operators for
the cavity mode of frequency $\omega_{c}$ and $H_{e,A}$ and $H_{e,B}$ refer to
the electronic Hamiltonians of molecules $A$ and $B$, respectively. The
perturbative Hamiltonian ($H^{1}$) includes the dipolar coupling between
molecules $A$ and $B$, in the spirit of London’s first derivation of van der
Waals interactions, and the light-matter coupling to a single cavity mode

$H^{1}=-\frac{\bm{d}_{A}\cdot\bm{d}_{B}}{R^{3}}+\frac{3(\bm{d}_{A}\cdot\bm{R})(\bm{d}_{B}\cdot\bm{R})}{R^{5}}+\lambda\sqrt{\frac{\omega_{c}}{2}}(\bm{\varepsilon}\cdot\Delta\bm{d}_{A}+\bm{\varepsilon}\cdot\Delta\bm{d}_{B})(b+b^{\dagger})+\frac{\lambda^{2}}{2}(\bm{\varepsilon}\cdot\Delta\bm{d}_{A}+\bm{\varepsilon}\cdot\Delta\bm{d}_{B})^{2}$
(S4)

where $\Delta\bm{d}_{A}=\bm{d}_{A}-\langle\bm{d}_{A}\rangle$ and
$\Delta\bm{d}_{B}=\bm{d}_{B}-\langle\bm{d}_{B}\rangle$ are the fluctuations of
molecule $A$ and molecule $B$’s dipoles, respectively and $\bm{d}_{A}$ and
$\bm{d}_{B}$ are the dipole operators for molecule $A$ and molecule $B$,
respectively. Recall that in this work we are working with van der Waals
molecules such that both molecules do not have permanent dipoles (i.e.
$\langle\bm{d}_{A}\rangle=\langle\bm{d}_{B}\rangle=0$).

The first-order correction to the energy is given by

$E^{1}=\matrixelement{g}{H^{1}}{g}$ (S5)

where $\ket{g}$ denotes the ground state of the total system,
$\ket{g}=\ket{g_{A}}\ket{g_{B}}\ket{g_{c}}$ where molecule $A$, molecule $B$,
and the cavity are in their ground states. In this illustrative perturbation
theory, we are interested in the asymptotic behavior for when molecule $A$ and
molecule $B$ are far away from one another; thus, the antisymmetry of the
total electronic wavefunctions is ignored. Substituting in Eq. S4 into Eq. S5,
we obtain

$\displaystyle E^{1}$
$\displaystyle=\frac{\lambda^{2}}{2}(\matrixelement{g_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)^{2}}{g_{A}}+\matrixelement{g_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)^{2}}{g_{B}})$
$\displaystyle=\frac{\lambda^{2}}{2}(E^{1}_{A}+E^{1}_{B})$ (S6)

where
$E^{1}_{A}=\matrixelement{g_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)^{2}}{g_{A}}$
and
$E^{1}_{B}=\matrixelement{g_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)^{2}}{g_{B}}$
are the dipole self-energies of molecule $A$ and molecule $B$, respectively.
In Eq. S6 we have used the facts that there are no photons in the ground state
of the cavity ($\matrixelement{g_{c}}{b^{\dagger}b}{g_{c}}=0$) and that for
van der Waals molecules, by definition, there is no permanent dipole
($\matrixelement{g_{A}}{\bm{d}_{A}}{g_{A}}=\langle\bm{d}_{A}\rangle=0$ and
$\matrixelement{g_{B}}{\bm{d}_{B}}{g_{B}}=\langle\bm{d}_{B}\rangle=0$). The
fact that molecules $A$ and $B$ do not have permanent dipoles allows us to
express $E^{1}_{A}$ and $E^{1}_{B}$ with a different formula, i.e.

$\displaystyle E^{1}_{A}$
$\displaystyle=\matrixelement{g_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)^{2}}{g_{A}}$
(S7)
$\displaystyle=\matrixelement{g_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)\hat{I}\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)}{g_{A}}$
$\displaystyle=\sum_{e_{A}}|\matrixelement{e_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)}{g_{A}}|^{2}\;\;,$

where $\ket{e_{A}}$ is an excited state of molecule $A$. An important
observation here is that both $E^{1}_{A}$ and $E^{1}_{B}$ are single molecule
terms and are always positive; we will return to these facts after deriving
the second-order energy correction.

The second-order correction to the energy is given by

$E^{2}=-\sum_{e}\frac{\left|\matrixelement{e}{H^{1}}{g}\right|^{2}}{E_{e}-E_{g}}$
(S8)

where $\ket{g}$ is the ground state of the bi-molecule system with energy
$E_{g}$ and $\ket{e}$ indicates an excited state of the bi-molecule system
with energy $E_{e}$. Substituting Eq. S4 into Eq. S8 along with some
simplifications we obtain the second-order correction to the energy to be

$\displaystyle E^{2}=$
$\displaystyle-\sum_{e_{A}e_{B}}\frac{\left|\matrixelement{e_{A}e_{B}}{V_{AB}}{g_{A}g_{B}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}-\lambda^{2}\sum_{e_{A}e_{B}}\frac{\matrixelement{e_{A}e_{B}}{V_{AB}}{g_{A}g_{B}}\matrixelement{e_{A}}{\bm{d_{A}}\cdot\bm{\varepsilon}}{g_{A}}\matrixelement{e_{B}}{\bm{d_{B}}\cdot\bm{\varepsilon}}{g_{B}}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}$
$\displaystyle-\frac{\lambda^{2}\omega_{c}}{2}\left[\sum_{e_{A}}\frac{\left|\matrixelement{e_{A}}{\bm{d_{A}}\cdot\bm{\varepsilon}}{g_{A}}\right|^{2}}{\omega_{c}+E_{e_{A}}-E_{g_{A}}}+\sum_{e_{B}}\frac{\left|\matrixelement{e_{B}}{\bm{d_{B}}\cdot\bm{\varepsilon}}{g_{B}}\right|^{2}}{\omega_{c}+E_{e_{B}}-E_{g_{B}}}\right]$
$\displaystyle-\frac{\lambda^{4}}{4}\left[\sum_{e_{A}}\frac{\left|\matrixelement{e_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)^{2}}{g_{A}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}}+\sum_{e_{B}}\frac{\left|\matrixelement{e_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)^{2}}{g_{B}}\right|^{2}}{E_{e_{B}}-E_{g_{B}}}+4\sum_{e_{A}e_{B}}\frac{\left|\matrixelement{e_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)}{g_{A}}\right|^{2}\left|\matrixelement{e_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)}{g_{B}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}\right]$
$\displaystyle=E^{2}_{AB,d^{0}}+\lambda^{2}E^{2}_{AB,d^{1}}+\frac{\lambda^{2}}{2}(E^{2}_{A,d^{1}}+E^{2}_{B,d^{1}})+\frac{\lambda^{4}}{4}(E^{2}_{A,d^{2}}+E^{2}_{B,d^{2}}+E^{2}_{AB,d^{2}})$
(S9)

where we defined

$V_{AB}=-\frac{\bm{d}_{A}\cdot\bm{d}_{B}}{R^{3}}+\frac{3(\bm{d}_{A}\cdot\bm{R})(\bm{d}_{B}\cdot\bm{R})}{R^{5}}\;\;\ldotp$
(S10)

$E^{2}_{AB,d^{0}}$, $E^{2}_{AB,d^{1}}$, $E^{2}_{A,d^{1}}$, $E^{2}_{B,d^{1}})$,
$E^{2}_{A,d^{2}}$, $E^{2}_{B,d^{2}}$, and $E^{2}_{AB,d^{2}}$ are defined as

$\displaystyle E^{2}_{AB,d^{0}}$
$\displaystyle=-\sum_{e_{A}e_{B}}\frac{\left|\matrixelement{e_{A}e_{B}}{V_{AB}}{g_{A}g_{B}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}$
(S11a) $\displaystyle E^{2}_{AB,d^{1}}$
$\displaystyle=-\sum_{e_{A}e_{B}}\frac{\matrixelement{e_{A}e_{B}}{V_{AB}}{g_{A}g_{B}}\matrixelement{e_{A}}{\bm{d_{A}}\cdot\bm{\varepsilon}}{g_{A}}\matrixelement{e_{B}}{\bm{d_{B}}\cdot\bm{\varepsilon}}{g_{B}}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}$
(S11b) $\displaystyle E^{2}_{A,d^{1}}$
$\displaystyle=-\omega_{c}\sum_{e_{A}}\frac{\left|\matrixelement{e_{A}}{\bm{d_{A}}\cdot\bm{\varepsilon}}{g_{A}}\right|^{2}}{\omega_{c}+E_{e_{A}}-E_{g_{A}}}$
(S11c) $\displaystyle E^{2}_{B,d^{1}}$
$\displaystyle=-\omega_{c}\sum_{e_{B}}\frac{\left|\matrixelement{e_{B}}{\bm{d_{B}}\cdot\bm{\varepsilon}}{g_{B}}\right|^{2}}{\omega_{c}+E_{e_{B}}-E_{g_{B}}}$
(S11d) $\displaystyle E^{2}_{A,d^{2}}$
$\displaystyle=-\sum_{e_{A}}\frac{\left|\matrixelement{e_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)^{2}}{g_{A}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}}$
(S11e) $\displaystyle E^{2}_{B,d^{2}}$
$\displaystyle=-\sum_{e_{B}}\frac{\left|\matrixelement{e_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)^{2}}{g_{B}}\right|^{2}}{E_{e_{B}}-E_{g_{B}}}$
(S11f) $\displaystyle E^{2}_{AB,d^{2}}$
$\displaystyle=-4\sum_{e_{A}e_{B}}\frac{\left|\matrixelement{e_{A}}{\left(\bm{d_{A}}\cdot\bm{\varepsilon}\right)}{g_{A}}\right|^{2}\left|\matrixelement{e_{B}}{\left(\bm{d_{B}}\cdot\bm{\varepsilon}\right)}{g_{B}}\right|^{2}}{E_{e_{A}}-E_{g_{A}}+E_{e_{B}}-E_{g_{B}}}\;\;,$
(S11g)

where $\ket{g_{A}}$ ($\ket{g_{B}}$) is the ground state of molecule $A$ ($B$)
with energy $E_{g_{A}}$ ($E_{g_{B}}$), $\ket{e_{A}}$ ($\ket{e_{B}}$) indicates
an excited state of molecule $A$ ($B$) with energy $E_{e_{A}}$ ($E_{e_{B}}$),
and $\matrixelement{e_{A}}{\bm{d}_{A}}{g_{A}}$
($\matrixelement{e_{B}}{\bm{d}_{B}}{g_{B}}$) is the transition dipole moment
of molecule $A$ ($B$) associated with the excited state. Eq. II is an
important result in this work, and the physical interpretation, origin, and
implications of each term are worth exploring in detail. $E^{2}_{AB,d^{0}}$ in
Eq. II is the typical attractive London dispersion interaction with its
prototypical $R^{-6}$ dependence (as each $V_{AB}$ scales with $R^{-3}$). The
remaining terms all arise from interactions through the cavity mode.
$E^{2}_{AB,d^{1}}$ contains a single $V_{AB}$ matrix element giving an
$R^{-3}$ of this term. Interestingly, this term also contains dot products of
transition dipole moments ($\matrixelement{e_{A}}{\bm{d}_{A}}{g_{A}}$) with
the cavity polarization vector ($\bm{\varepsilon}$). This $R^{-3}$ term is
central to this work as it says that van der Waals molecules inside a cavity
have this interesting interaction length scale that also has unique, coupled
molecule-molecule and molecular-cavity angle dependencies. $E^{2}_{A,d^{1}}$
and $E^{2}_{B,d^{1}}$are very similar to $E^{1}_{A}$ and $E^{1}_{B}$ except
that $E^{2}_{A,d^{1}}$ and $E^{2}_{B,d^{1}}$ arise from the bilinear coupling
term and have the opposite sign as $E^{1}_{A}$ and $E^{1}_{B}$. Specifically,
to second-order in the coupling $\lambda$, the one-body energy (e.g. molecule
$A$) is given by

$\displaystyle E_{A}^{\text{cavity}}$ $\displaystyle=E_{A}^{\text{no
cavity}}+\frac{\lambda^{2}}{2}(E^{1}_{A}+E^{2}_{A,d^{1}})$ (S12)
$\displaystyle=E_{A}^{\text{no
cavity}}+\frac{\lambda^{2}}{2}\sum_{e_{A}}\left|\matrixelement{e_{A}}{\bm{d}_{A}\cdot\bm{\varepsilon}}{g_{A}}\right|^{2}$
$\displaystyle\quad-\frac{\lambda^{2}\omega_{c}}{2}\sum_{e_{A}}\frac{\left|\matrixelement{e_{A}}{\bm{d}_{A}\cdot\bm{\varepsilon}}{g_{A}}\right|^{2}}{w_{c}+E_{e_{A}}-E_{g_{A}}}\;\;\ldotp$

A similar energy term can be derived for molecule $B$ as well. We want to
emphasize that $E^{1}_{A}$ arises from the dipole self-energy term in first-
order perturbation theory (Eq. S6) and $E^{2}_{A,d^{1}}$ arises from the
bilinear coupling term in second-order perturbation theory (Eq. II).
Interestingly, $E^{1}_{A}$ and $E^{2}_{A,d^{1}}$ only exactly cancel if the
cavity frequency is much larger than the electronic transition energies
($\omega_{c}\gg E_{e_{A}}-E_{g_{A}}$). Thus, for H2 molecules with a cavity in
the electronic regime ($\omega_{c}=13.6$ eV here) the total energy of a single
molecule ends up increasing with $\lambda^{2}$ (main text Fig. 3A). For H2
molecules, the one-body energy reaches a minimum when the molecular bond is
perpendicular to the cavity polarization vector
($\theta_{A\varepsilon}=\frac{\pi}{2}$). Intuitively, this occurs because H2
is most polarizable along its bond axis which leads to
$\sum_{e_{A}}\left|\matrixelement{e_{A}}{\bm{d}_{A}\cdot\bm{\varepsilon}}{g_{A}}\right|^{2}/(E_{e_{A}}-E_{g_{A}})=\bm{\varepsilon}^{T}\bm{\alpha}\bm{\varepsilon}$
being largest when $\theta_{A\varepsilon}=0,\pi$.

$E^{2}_{A,d^{2}}$, $E^{2}_{B,d^{2}}$, and $E^{2}_{AB,d^{2}}$ arise from two
factors of the dipole self-energy part of Eq. S4 and, thus, scale with
$\lambda^{4}$. While $E^{2}_{A,d^{2}}$ and $E^{2}_{B,d^{2}}$ are corrections
to the one-body energies, $E^{2}_{AB,d^{2}}$ impacts the two-body energies
(i.e. intermolecular interaction energy). Furthermore, this term has no $R$
dependence, and, thus, $E^{2}_{AB,d^{2}}$ is the first term that we have
discussed that gives rise to the collective orientational order reported in
the main text. The magnitude of this term is greatest when both molecules have
their bonds oriented along the cavity polarization vector
($\bm{\varepsilon}$), because
$\bm{\varepsilon}^{T}\bm{\alpha}_{A}\bm{\varepsilon}$ and
$\bm{\varepsilon}^{T}\bm{\alpha}_{B}\bm{\varepsilon}$ are both largest in the
case which both of their bonds are oriented parallel to $\bm{\varepsilon}$.
And because of the negative sign in front of this infinite range interaction
term, it contributes to lowering the energy of molecular configurations in
which the molecular bonds of the hydrogen molecules are oriented parallel to
the cavity polarization vector, as shown in Fig. 3C of the main text.

##  III Many-body Interactions

The many-body expansion,

$E=\sum_{A}E_{A}+\sum_{AB}E_{AB}+\sum_{ABC}E_{ABC}+\dots$ (S13)

is a routinely used expansion for modeling and gaining insight into
intermolecular forces.Dahlke and Truhlar (2007) For van der Waals type
intermolecular forces, the higher-order interactions such as $E_{ABC}$ quickly
become negligible with distance and they can be assumed to be much smaller
than the lower-order terms at large distances. QED electronic structure
calculations allow us to test if the three-body and higher-order terms can be
ignored for the strong light-matter coupling cavity QED Hamiltonian with
similar parameters used in the calculations of the main text. In Table S1 and
Fig. S2, we show the intermolecular interactions for molecules separated far
apart, $25$ Å. As expected, QED-HF does not capture the dynamic correlation
and cannot describe the intermolecular forces arising from neither the cavity
nor the van der Waals forces. QED-CCSD-1 captures the dynamic correlation, but
the sign of the two-body interaction is not consistent with QED-FCI. Adding
just one more term to the cluster operator of QED-CCSD-1, the two-photon
$(b^{\dagger})^{2}$ term in QED-CCSD-12-SD1, yields a sufficient description
of the two-body interactions. For QED-CCSD-12-SD1, we find that the higher-
order terms quickly approach zero even for the very strong coupling
$\lambda=0.1$ a.u. From perturbation theory, we find that the $N$-body
interactions are sensitive to the light-matter coupling strength and scale as
$\lambda^{2N}$ (see Fig. S2).

A few additional key points about the many-body expansion of van der Waals
interactions in the context of the nonrelativistic cavity QED Hamiltonian
given in Eq. S1 are worth mentioning here. Because the three-body interactions
have opposite sign to the two-body interactions (Table S1), we expect that the
collective orientational order induced by the infinite range cavity-induced
interactions would be reduced by including the three-body terms in the
molecular dynamics simulations. While the three-body terms are insignificant
on a per interaction basis, the lack of distance ($R$) dependence in the
cavity-induced interactions, see Eq. II, results in all molecules in the
simulation interacting with all other molecules independent of how far away
they are from each other. In a simulation with $n$ molecules, there are
$n(n-1)/2\sim n^{2}$ two-body interactions, $n(n-1)(n-2)/6\sim n^{3}$ three-
body interactions, and similarly for higher-order terms (Table S2). Therefore,
there must exist a number of molecules where the total three-body energy is
larger than the total two-body energy. This makes it very challenging to
extrapolate our results to truly macroscopic systems. Extending these
microscopic equations and calculations to truly macroscopic systems remains an
open question.

Method | 1-body | 2-body | 3-body | 4-body
---|---|---|---|---
QED-HF | 204.9 | 0.0000 | 0.0000 | 0.0000
QED-CCSD-1 | 107.5 | 0.3238 | -0.0571 | 0.0042
QED-CCSD-12-SD1 | 107.1 | -0.5600 | 0.0104 | -0.0004
QED-FCI-5 | 106.7 | -0.6601 | $\ldots$ | $\ldots$
Table Table S1: Cavity-induced $N$-body effects for different QED electronic structure methodologies with $\lambda=0.1$ a.u. The cavity energy is $\hbar\omega_{c}=13.6$ eV and polarization perpendicular to all molecules. The molecules are placed on the edges of a line ($E_{AB}$), equilateral triangle ($E_{ABC}$) and square ($E_{ABCD}$), all with side lengths of $25$ Å. All numbers in the table are meV. QED-FCI-5 is too computationally expensive for more than two H2 molecules in the aug-cc-pVDZ basis set. Figure Fig. S2: $N$-body effects for different coupling strengths $\lambda$. All calculations are performed on $N$ H2 molecules with QED-CCSD-12-SD1. The cavity energy is $\hbar\omega_{c}=13.6$ eV and polarization perpendicular to all molecules. The molecules are placed on the edges of a line ($E_{AB}$), equilateral triangle ($E_{ABC}$) and square ($E_{ABCD}$), all with side lengths of $25$ Å. | 1-body | 2-body | 3-body | 4-body
---|---|---|---|---
Scaling with coupling | $\lambda^{2}$ | $\lambda^{4}$ | $\lambda^{6}$ | $\lambda^{8}$
Number of terms | $n\choose 1$ | $n\choose 2$ | $n\choose 3$ | $n\choose 4$
Table Table S2: The number of interactions and scaling of the cavity-induced
interaction energy in the $N$th body of the $N$-body expansion for a system
with $n$ molecules.

##  IV Molecular Dynamics

###  IV.1 Training Potential Energy Functions for Simulating Fluids of H2

####  IV.1.1 Neural Network-based Pairwise Interactions

We developed neural network-based potential energy functions (NNPs) for the
pairwise interaction of a pair of hydrogen molecules using ${\it ab~{}initio}$
energy data with CCSD, FCI, QED-CCSD-12-SD1, and QED-FCI levels of theory. The
potential energy functions have the forms,

$E_{\rm AB}^{\rm no~{}cavity}=c_{\rm
exp}\exp(-aR)-\frac{c_{6}\\{\theta\\}}{R^{6}}$ (S14) $E_{\rm AB}^{\rm
cavity}=E_{\rm 2b}^{\rm
no~{}cavity}-\frac{c_{3}\\{\theta\\}}{R^{3}}+\frac{c_{0}\\{\theta\\}}{R^{0}}$
(S15)

where $c_{\rm exp},{a},c_{6},c_{3},c_{0}$ are represented by neural networks
(NNs). Each NN takes symmetry preserved features of a pair of molecules as
input. Symmetry preserved features that have been selected as the input for
the machine learning (ML) model to get the pairwise interaction energy are
shown pictorially in Fig. S3 and are listed in Table S3. In the case without
the cavity field, the interaction energies are obtained using the input
features $\theta_{{\bf R}A},\theta_{{\bf R}B},\theta_{AB},{\left\|\bf
R\right\|}$. For the cavity case, additional terms that depend on the cavity
polarization vector are added. In particular,
$\theta_{A\varepsilon},\theta_{B\varepsilon},\text{ and }\theta_{{\bf
R}\varepsilon}$ are added and $\left\|\bf R\right\|$ is replaced by $R_{\rm
cap}$ and $R_{\rm cap}={\rm C}\tanh(\left\|\bf R\right\|/{\rm C})$, where C is
a cutoff distance. In order to account for molecular and exchange symmetries,
$\cos 2\theta$ and $\sin 2\theta$ are used for any
$\theta\in\Theta\equiv\\{\theta_{{\bf R}A},\theta_{{\bf R}B},\theta_{\rm
AB},\theta_{A\varepsilon},\theta_{B\varepsilon},\theta_{{\bf
R}\varepsilon}\\}$. For each of $c_{\rm exp}$, $a$, $c_{6}$, $c_{3}$, we are
using F($\Theta,R_{\rm cap}$)+F($\tilde{\Theta},\tilde{R}_{\rm cap}$) where
$\tilde{\Theta}$ and $\tilde{R}_{\rm cap}$ are calculated by switching the
index of the two molecules. For $c_{0}$, only Type 1 features as tabulated in
Table S3 were used.

The neural network model has four fully-connected layers including a linear
output layer. The other three linear layers have CELU activation
functions.Barron (2017) The number of neurons per layer is 64 in our model. To
train the model, we used energy data points of pair configurations that are
generated using a classical MD simulation of liquid H2. $10^{5}$ pair
configurations generated by MD simulation were used to compute energies with
CCSD level of theory for training model when no cavity is present. While the
pair configurations generated by MD simulation were good enough to train a
model without a cavity, long range pair configurations are extremely important
to train the model with a cavity. Similarly, short range pair configurations
are very crucial to accurately reproduce the corrected short range repulsion
energies in the potential energy functions in the presence of a cavity. While
MD of liquid H2 produces good random configurations with various possible
orientations, the probability of finding short range pair configurations is
low in an MD simulation. In order to include sufficient number of
configurations at short range, we randomly select $10\%$ of the total
configurations obtained from MD simulation of liquid H2 molecules and scale
the intermolecular distance to be within $2-5$ Å. A similar strategy was
followed to generate very long range configurations between $18-90$ Å for
$10\%$ of the total configurations. A total of $121,000$ data points,
including both the additional short range and long range configurations, were
used to the train the NN model to the QED-CCSD-12-SD1 calculated energies in
the cavity case. For training using the QED-FCI calculated data, we use a
smaller data set of $30,000$ calculated energies. In order to train the model
on this smaller data set, we initialize each NN with the parameters obtained
from our QED-CCSD-12-SD1 fits, which was trained using a larger data set of
$121,000$ calculated energies. We use the Adam optimizer Kingma and Ba (2014)
with $\beta_{1}=0.90$ and $\beta_{2}=0.99$. And we utilize a constant learning
rate of $10^{-5}$ and a batch size of $32$. $90\%$ of the total data points
were used in the training data set and the remaining $10\%$ were used as a
test data set. All training and testing protocols were implemented with
PyTorch.Paszke  _et al._ (2019)

Type of feature | Features
---|---
Type 1 |  $\cos 2\theta_{A\varepsilon}$, $\sin 2\theta_{A\varepsilon}$,$\cos 2\theta_{B\varepsilon}$, $\sin 2\theta_{B\varepsilon}$
Type 2 |  $\cos 2\theta_{\bf R\varepsilon}$, $\sin 2\theta_{\bf R\varepsilon}$, $\cos 2\theta_{{\bf R}A}$, $\sin 2\theta_{{\bf R}A}$
|  $\cos 2\theta_{{\bf R}B}$, $\sin 2\theta_{{\bf R}B}$, $\cos
2\theta_{{AB}}$, $\sin 2\theta_{{AB}}$
Type 3 | $\rm C\tanh(\left\|\bf R\right\|/C)$
Table Table S3: Input features involved in the energy contributions. Figure
Fig. S3: Symmetry preserved features that are considered while generating the
pair interaction potential using a neural network based machine learning model
are shown here. Various angles between a pair of molecule which are considered
as input features are shown. $\bf R$ is the distance vector of the center of
mass (COM) of molecule $A$ and molecule $B$. $\varepsilon$ represents the
cavity polarization. Orientation of the molecules are completely specified by
various angles $\\{\theta\\}$. Figure Fig. S4: Energy of a single H2 molecule
inside a cavity with respect to cavity polarization vector, ${\varepsilon}$
using it ab initio QED-CCSD-12-SD1 and ML. Single molecular energy at
${\varepsilon}=0.0$ was set to zero while plotting energies of both QED-
CCSD-12-SD1 and ML.

The energies of the ab initio (CCSD) calculations and the ML predicted
energies of the pairs of molecules without a cavity field are shown in the
Fig. S9A. A linearity plot shows the accuracy of the predicted energy using
our ML model. Apart from the linearity plot, we scanned potential energy
curves for a few selected orientations of pairs of molecules. These results
show that the ML predicted potential energy curves for pairs of hydrogen
molecules are in good agreement with the potential energy curves obtained from
ab initio calculations. These plots are shown in Fig. S9B. A linearity plot
comparing the ab initio (QED-CCSD-12-SD1) calculations and the ML predicted
energies with the cavity field turned on are shown in Fig. S10A. Potential
energy curves (Fig. S10B) were scanned for D2h configuration of a pair of
molecules along three different cavity polarization directions with respect to
the molecular bond axis. These plots shows that our ML model accurately
reproduces the ab initio potential energy curves.

####  IV.1.2 Single Molecule Potential Energies

Single molecule potential energies involve intra-molecular chemical bonds and
the cavity-modified single molecule contributions. Intra-molecular chemical
bonds were modeled within the harmonic approximation. We like to emphasize
that the intra-molecular interaction energy does not play a significant role
in determining the properties that we focused on in this study.

Single molecule energies in the presence of a cavity field is important.
Training of the cavity-modified single molecule energies has been done with a
linear regression method. The following form of energy function is trained for
the single molecule energies,

$E_{\rm A}=\sum_{n=1}^{3}C_{n}\sin 2n{\theta}+\sum_{n=0}^{2}D_{n}\cos
2n{\theta}$ (S16)

where $\theta$ is the angle between the molecular bond axis and the cavity
polarization vector. $C_{n}$ and $D_{n}$ are the trainable parameters. Fig. S4
shows the accuracy of fitting single molecular energies with respect to the ab
initio, QED-CCSD-12-SD1 calculations.

###  IV.2 Molecular Dynamics

Molecular dynamics (MD) simulations were used to compute the statistical
properties of fluids of $H_{2}$ molecules at $70$ K by employing the potential
energy functions, generated by our machine learning models. For computing the
statistical behaviour of the system both classical MD and path integral MD
(PIMD) were used.

####  IV.2.1 Classical Molecular Dynamics

NVT ensemble MD simulations were carried out using Langevin dynamics with a
time step of $1.0$ femtosecond (fs) and the friction coefficient for the
Langevin dynamics was chosen $0.0005$ a.u (20.7 ps-1). Random initial atomic
velocities and random initial positions were provided to run MD. In order to
use ML potentials generated with PyTorch, we also implement the MD engine with
PyTorch. The integrator used here is described in Ref. Bussi and Parrinello
(2007). Forces were computed using the PyTorch autograd module and the PyTorch
MD simulations were performed using GPUs.

Since we are simulating a fluid system, the system was confined within a
spherical volume, similar to a cluster of molecules. In practice, a stiff
harmonic potential was used to confine the center of the mass of each molecule
within a spherical volume with radius $R_{\rm c}$ (see Fig. S5). Adopting such
a boundary condition was necessary in order to account for non-decaying nature
of the pair interaction potential inside of an optical cavity. In order to
simulate various different system sizes, $R_{\rm c}$ is scaled appropriately
to preserve the overall molecular density.

####  IV.2.2 Path Integral Molecular Dynamics

In the previous section, we discussed the MD simulations in which the nuclei
were considered as classical particles. However, for light nuclei such as
hydrogen atoms, this assumption could lead to serious problems in predicting
the statistical properties because of strong quantum nuclei effects,
especially at low temperatures. In order to account for quantum nuclei effects
in our MD simulations, we performed path integral molecular dynamics (PIMD)
simulations.

Usually PIMD simulations require a large number of beads to converge
thermodynamics properties at low temperatures. Herein, we used the generalized
Langevin equation (GLE) in PIMD, which can significantly reduce the number of
beads.Ceriotti  _et al._ (2009, 2010a, 2011) In the GLE formulation,Ceriotti
_et al._ (2010b) each bead of the simulated system is coupled to several
extended degrees of freedom with an appropriate drift matrix and a diffusion
matrix to approximate a friction kernel function. We used $8$ extra degrees of
freedom in GLE and the drift matrix and diffusion matrix used in GLE were
generated by an online tool called GLE4MD (http://gle4md.org/) with the
maximum physical frequency set to $\omega_{\rm max}=9608$ cm-1. With the GLE
formulation, we observed that using $32$ beads are able to converge the
simulations whereas more than $128$ beads are needed to converge the results
without the GLE formulation. We have developed an interface to i-PI Ceriotti
_et al._ (2014) to run the PIMD simulations using our ML potentials.

###  IV.3 Radial Distribution Functions

Figure Fig. S5: Schematic diagram of the radius cutoff that are used in
computing radial distribution functions. $R_{\rm c}$ is the distance at which
a high energy potential barrier has been applied. $R_{\rm 1}$ is the radius of
core region where surface effects due to the spherical boundary are minimal
and molecules found within the radius of $R_{\rm 2}$ are used to compute the
histogram of pairwise distance for the calculations of the radial distribution
functions.

The radial distribution functions (g(r)) of fluid of H2 molecules are computed
from the PIMD trajectories of $1,000$ molecules. As the system we simulated
has a spherical volume without any periodic boundary, computing a bulk-like
g(r) (i.e. a g(r) that converges to $1$ in the long distance limit) is not
straightforward. In order to compute g(r) from such a spherical system, the
following steps are taken. First, a bulk-like core region is chosen within a
certain cutoff distance $R_{1}$.

$\bar{h}(\left|\bf r\right|)=\frac{1}{N_{1}}\sum_{i,R_{i}<R_{1}}h(\left|{\bf
r}-{\bf r}_{i}\right|)$ (S17)

For the $i^{\rm th}$ molecule located at ${\bf r}_{i}$ with $R_{i}=\left|{\bf
r}_{i}\right|<R_{1}$, $h(\left|{\bf r}-{\bf r}_{i}\right|)$ is the histogram
of all distance between any other molecules and the $i^{th}$ molecule with
$(\left|{\bf r}-{\bf r}_{i}\right|)<R_{2}$, $R_{1}+R_{2}<R_{\rm c}$ and
$N_{1}$ is the number of molecules inside $R_{1}$. Second, the average over
each frame of MD or PIMD as well as the average over number of beads was
computed in the calculations of the radial distribution functions. Lastly, the
averaged $\bar{h}(\left|\bf r\right|)$ was normalized by the average density
and $4\pi r^{2}$. In this study, $R_{1}=6.0$ Å and $R_{2}=12$ Å was used.

###  IV.4 Angular Distribution Functions

We also computed angular distribution functions for the angle between the
molecular bond axis of molecule $A$ and the molecular bond axis of molecule
$B$ ($\theta_{\rm AB}$) and angular distribution functions for the angle
between the molecular bond axis of molecule $A$ and the cavity polarization
vector ($\theta_{\rm A\varepsilon}$). The probability distributions of
$\theta_{\rm AB}$ and $\theta_{\rm A\varepsilon}$ are proportional to
sin($\theta_{\rm AB}$) and sin($\theta_{\rm A\varepsilon}$), respectively, if
molecules A and B can rotate freely without any interactions. In order to
emphasize the energy contribution, we computed the potentials of mean force by
scaling the probability distributions of $\theta_{\rm AB}$ and $\theta_{\rm
A\varepsilon}$ with their corresponding sine functions. In the case of PIMD,
the average over each frame and the average over the number of beads are
considered when computing the histograms.

##  V Additional Results

###  V.1 Comparison of Radial Distribution Functions

We compute the radial distribution function at three different situations when
(1) cavity polarization is not active, (2) cavity-modified one-body term is
active but cavity modified two-body term is not active, and (3) both cavity
modified one-body and two-body terms are active. We have observed
differentiable changes in radial distribution function for three different
situations. This indicates the difference in equilibrium structure when cavity
polarization is on. The results are shown in Fig. S12.

###  V.2 Comparison of Classical MD and PIMD

In this section we compare the results of our classical MD and the PIMD
simulations with $\lambda=0.1$ a.u. Based on Fig. S6, it is evident that
classical MD and PIMD qualitatively follow the same trend when angular
distribution function of $\theta_{\rm A\varepsilon}$ and $\theta_{\rm AB}$ are
compared. In particular, one observes a strong orientational alignment of the
molecules along direction of the cavity polarization vector occurring inside
of an optical cavity. Inclusion of nuclear quantum effects does not change the
overall conclusion. However, the extent of alignment of the molecules inside
the cavity in our PIMD simulations is considerably reduced compared to our
classical MD simulations.

Figure Fig. S6: Angular distribution functions of molecular bond axis of
molecule $A$ to the molecular bond axis of molecule $B$ ($\theta_{\rm AB}$)
and angular distribution functions of molecular bond axis of molecule $A$ to
the cavity polarization vector ($\theta_{\rm A\varepsilon}$) for $1,000$ H2
molecules of a (A) classical MD simulation and (B) PIMD simulation are shown.
Pair interaction potentials used for the MD simulation were obtained by
training an ML model with the calculated energies from QED-CCSD-12-SD1 level
of theory.

###  V.3 Comparison of QED-FCI-5 and QED-CCSD-12-SD1

Here we compare our results of classical MD simulations using the ML
potentials obtained from QED-FCI-5 and QED-CCSD-12-SD1 calculations. As
summarized in Fig. S7, we see that classical MD with ML potentials that are
obtained from the two different levels of $\it ab~{}initio$ calculations
qualitatively match each other. However, the intensities in the angular
distribution functions of $\theta_{\rm A\varepsilon}$ and $\theta_{\rm AB}$
for the two cases are different. These differences are due to the quantitative
differences in predicting the interaction energies using these two methods
(see Fig. S1).

Figure Fig. S7: Angular distribution functions of molecular bond axis of
molecule $A$ to the molecular bond axis of molecule $B$ ($\theta_{\rm AB}$)
and angular distribution functions of molecular bond axis of molecule $A$ to
the cavity polarization vector ($\theta_{\rm A\varepsilon}$) for $1,000$ H2
molecules of a classical MD trajectory with the NN potentials obtained from
training the ML model on (A) QED-CCSD-12-SD1 and (B) QED-FCI-5 data sets.

###  V.4 $\lambda$ Dependent Molecular Alignment

Two different $\lambda$ values were considered in our study. In the main text,
we focused our discussion on the results with $\lambda=0.1$ a.u. In this
section, we study the properties of a system with $\lambda=0.02$ a.u. and
compare these results with the results obtained using $\lambda=0.1$ a.u.

In order to train a model with $\lambda=0.02$ a.u. important NN parameters for
$c_{0}$ and $c_{3}$ were transferred and scaled from our training model with
$\lambda=0.1$ a.u. together with the perturbation theory analysis. The
accuracy of the model has been tested by plotting the energies obtained from
the NNPs against the ab initio energies. A linearity plot is obtained as shown
in Fig. S11A. Additionally, scanned potential energy curves of several
selected pair configurations are in good agreement with ab initio potential
energy curves. Some of these plots are shown in Fig. S11B. The accuracy of our
ML model is further justified with in Fig. S11C, where we show that our ML
model correctly predicts the long range interaction energy with different
directions of the cavity polarization vector.

A significant difference in the angular distribution functions of $\theta_{\rm
A\varepsilon}$ is observed when the results of two different $\lambda$ values
are compared for $1,000$ H2 molecules. The distribution function of
$\theta_{\rm A\varepsilon}$ for $1,000$ H2 molecules with $\lambda=0.02$ a.u.
(Fig. S8A) shows molecular alignment perpendicular to the cavity polarization
($\theta_{\rm A\varepsilon}=\frac{\pi}{2}$). On the other hand, we observe in
Fig. S6A that the angular distribution function of $\theta_{\rm A\varepsilon}$
is maximized in the direction of cavity polarization vector ($\theta_{\rm
A\varepsilon}=0,\pi$) when $\lambda=0.1$ a.u. This can be explained from our
perturbation theory analysis where we showed that the cavity-modifications to
the single molecule energies scale with $\lambda^{2}$ and the extremely long
range pairwise interaction scales with $\lambda^{4}$. Thus, the importance of
the pairwise interaction decreases much faster than the single molecule energy
contribution as $\lambda$ decreases. In this particular example of $1,000$ H2
molecules with $\lambda=0.02$ a.u., the single molecule energy dominates
whereas, with $\lambda=0.1$ a.u., the pairwise interaction energy dominates.
$\theta_{\rm AB}$ qualitatively follow the same trend as we observed for
$1,000$ H2 molecules with $\lambda=0.1$ a.u.; however, the intensity of the
peak is reduced which suggests a weaker synchronization of molecular
orientations. This is shown in the inset of Fig. S8A.

From the above discussion, we understand that the energy contributions from a
single molecule can be altered by (1) changing the number of molecules with a
fixed $\lambda$, and (2) changing the value of $\lambda$ for a fix number of
molecules. We ran simulations considering these two possibilities. For the
first possibility, we reduced the number of molecules from $1,000$ to $108$
while keeping $\lambda$ equal to $0.1$ a.u., and we compute the angular
distribution function for $\theta_{\rm A\varepsilon}$. We find that in the
$108$ molecule simulation the preferential alignment of the molecules is
perpendicular to the cavity polarization vector, which is opposite to the
alignment of $1,000$ molecules with $\lambda=0.1$ a.u. (aligned parallel to
the cavity polarization vector). These results are shown in Fig. S6A and Fig.
S8B. For the second possibility, we simulate $1,000$ molecules with a reduced
value of $\lambda=0.02$ a.u. The angular distribution function of $\theta_{\rm
A\varepsilon}$ in this simulation is qualitatively similar to the results
obtained in the first possibility with the molecular alignment perpendicular
to the cavity polarization vector (see Fig. S6A and Fig. S8B). All of our
numerical simulation results reported in this section further confirm the
conceptual validity of our perturbation theory analysis.

Figure Fig. S8: Angular distribution functions of molecular bond axis of
molecule $A$ to the molecular bond axis of molecule $B$ ($\theta_{\rm AB}$)
and angular distribution functions of molecular bond axis of molecule $A$ to
the cavity polarization vector ($\theta_{\rm A\varepsilon}$) for $1,000$ H2
molecules of a classical MD trajectory with the NNPs obtained from the
training ML model on (A) QED-CCSD-12-SD1 and $\lambda=0.02$ a.u. coupling
constant are shown. A zoom-in figure of $\theta_{\rm A\varepsilon}$ is shown
in the inset. (B) Angular distribution functions of molecular bond axis of
molecule $A$ to the cavity polarization vector ($\theta_{\rm A\varepsilon}$)
of $108$ molecules with $\lambda=0.1$ a.u. (dashed line) and $1,000$ molecules
with $\lambda=0.02$ a.u. (solid line) are shown. Figure Fig. S9: (A) Pairwise
interaction energies obtained from ab initio, CCSD calculation (without
cavity) and ML predicted energies are plotted. (B) Scanned potential energy
curve for D2h, C2v and D∞h configuration of a pair of molecules using NNPs and
from ab initio calculation are shown.  Figure Fig. S10: (A) Pairwise
interaction energies obtained from ab initio, QED-CCSD-12-SD1 calculation
(with cavity) and ML predicted energies are plotted. (B) Scanned potential
energy curve for D2h configuration with three different direction of cavity
polarization using NNPs and from ab initio calculation are shown. Figure Fig.
S11: (A) Pairwise interaction energies obtained from ab initio, QED-
CCSD-12-SD1 calculation and ML predicted energies with $\lambda=0.02$ a.u. are
plotted. (B) Scanned potential energy curves for D2h configuration of a pair
of molecules using NNPs and from ab initio calculation are shown. Distance
($R$) between molecule $A$ and molecule $B$ over which potential energy is
scanned is shown in the inset of the figure. (C) Scanned potential energy
curves for D2h configuration at the long range are shown. ML model can
accurately distinguish different configurations at long distance. Figure Fig.
S12: Radial distribution function generated using PIMD trajectory with 1000 H2
molecules using pair potential obtained through a ML training on ab initio
calculation with QED-CCSD-12-SD1 and $\lambda=0.1$ a.u. Figure Fig. S13: (A-C)
Snapshots taken at thermal equilibrium from the path integral molecular
dynamic (PIMD) simulations of 1000 H2 molecules in the case of (A) no cavity
(orange), (B) cavity-modified one-body term but no cavity two-body term
(green), and (C) cavity-modified one-body and two-body terms (blue). For these
three cases, the (D) molecular bond axis of molecule $A$ to molecular bond
axis of molecule $B$ ($\theta_{AB}$) angular probability distribution
function, $P\left(\theta_{AB}\right)$ and (E) molecular bond axis to cavity
polarization vector ($\theta_{A\varepsilon}$), angular probability
distribution function, $P\left(\theta_{A\varepsilon}\right)$, are shown. (F)
molecular bond axis to cavity polarization vector ($\theta_{A\varepsilon}$),
angular probability distribution function,
$P\left(\theta_{A\varepsilon}\right)$, are shown for four different
simulations containing different numbers of H2 molecules. All PIMD simulations
shown in this figure were performed using neural networks trained with CCSD
(no cavity) or QED-CCSD-12-SD1 with $\lambda=0.1$ a.u. (cavity) calculated
energies.

## References

  * Hobza and Šponer (2002) P. Hobza and J. Šponer, J. Am. Chem. Soc. 124, 11802 (2002).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, aac9439 (2016).
  * Sternbach  _et al._ (2021) A. J. Sternbach, S. H. Chae, S. Latini, A. A. Rikhter, Y. Shao, B. Li, D. Rhodes, B. Kim, P. J. Schuck, X. Xu, X. Y. Zhu, R. D. Averitt, J. Hone, M. M. Fogler, A. Rubio, and D. N. Basov, Science 371, 617 (2021).
  * Maitland  _et al._ (1981) G. C. Maitland, G. D. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham,  _Intermolecular Forces: Their Origin and Determination_ (Oxford University Press, USA, 1981).
  * Stone (2013) A. Stone,  _The Theory of Intermolecular Forces_, 2nd ed. (Oxford University Press, Oxford, 2013) p. 352.
  * London (1937) F. London, Trans. Faraday Soc. 33, 8b (1937).
  * Halgren (1992) T. A. Halgren, J. Am. Chem. Soc. 114, 7827 (1992).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Thirunamachandran (1980) T. Thirunamachandran, Mol. Phys. 40, 393 (1980).
  * Milonni and Smith (1996) P. W. Milonni and A. Smith, Phys. Rev. A 53, 3484 (1996).
  * Sherkunov (2009) Y. Sherkunov, J. Phys. Conf. Ser. 161, 012041 (2009).
  * Fiscelli  _et al._ (2020) G. Fiscelli, L. Rizzuto, and R. Passante, Phys. Rev. Lett. 124, 013604 (2020).
  * Haugland  _et al._ (2021) T. S. Haugland, C. Schäfer, E. Ronca, A. Rubio, and H. Koch, J. Chem. Phys. 154, 094113 (2021).
  * Ribeiro  _et al._ (2018) R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem. Sci. 9, 6325 (2018).
  * Rivera  _et al._ (2019) N. Rivera, J. Flick, and P. Narang, Phys. Rev. Lett. 122, 193603 (2019).
  * Thomas  _et al._ (2019) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Science 363, 615 (2019).
  * Li  _et al._ (2020) T. E. Li, J. E. Subotnik, and A. Nitzan, Proc. Natl. Acad. Sci. U. S. A. 117, 18324 (2020).
  * Garcia-Vidal  _et al._ (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, eabd0336 (2021).
  * Li  _et al._ (2021a) T. E. Li, A. Nitzan, and J. E. Subotnik, Angew. Chemie 133, 15661 (2021a).
  * Vahala (2003) K. J. Vahala, Nature 424, 839 (2003).
  * Cortese  _et al._ (2017) E. Cortese, P. G. Lagoudakis, and S. De Liberato, Phys. Rev. Lett. 119, 043604 (2017).
  * Joseph  _et al._ (2021) K. Joseph, S. Kushida, E. Smarsly, D. Ihiawakrim, A. Thomas, G. L. Paravicini-Bagliani, K. Nagarajan, R. Vergauwe, E. Devaux, O. Ersen, U. H. F. Bunz, and T. W. Ebbesen, Angew. Chem. Int. Ed. 60, 19665 (2021).
  * Fukushima  _et al._ (2022) T. Fukushima, S. Yoshimitsu, and K. Murakoshi, J. Am. Chem. Soc. 144, 12177 (2022).
  * Sandeep  _et al._ (2022) K. Sandeep, K. Joseph, J. Gautier, K. Nagarajan, M. Sujith, K. G. Thomas, and T. W. Ebbesen, J. Phys. Chem. Lett. 13, 1209 (2022).
  * Galego  _et al._ (2015) J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 5, 41022 (2015).
  * Lacombe  _et al._ (2019) L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys. Rev. Lett. 123, 083201 (2019).
  * Fregoni  _et al._ (2022) J. Fregoni, F. J. Garcia-Vidal, and J. Feist, ACS Photonics 9, 1096 (2022).
  * Haugland  _et al._ (2020) T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio, and H. Koch, Phys. Rev. X 10, 041043 (2020).
  * Zhong  _et al._ (2016) X. Zhong, T. Chervy, S. Wang, J. George, A. Thomas, J. A. Hutchison, E. Devaux, C. Genet, and T. W. Ebbesen, Angew. Chem. Int. Ed. 55, 6202 (2016).
  * Du  _et al._ (2018) M. Du, L. A. Martínez-Martínez, R. F. Ribeiro, Z. Hu, V. M. Menon, and J. Yuen-Zhou, Chem. Sci. 9, 6659 (2018).
  * Xiang  _et al._ (2020) B. Xiang, R. F. Ribeiro, M. Du, L. Chen, Z. Yang, J. Wang, J. Yuen-Zhou, and W. Xiong, Science 368, 665 (2020).
  * Herrera and Spano (2016) F. Herrera and F. C. Spano, Phys. Rev. Lett. 116, 238301 (2016).
  * Yang and Cao (2021) P. Y. Yang and J. Cao, J. Phys. Chem. Lett. 12, 9531 (2021).
  * Li  _et al._ (2021b) X. Li, A. Mandal, and P. Huo, Nat. Commun. 12, 1315 (2021b).
  * Simpkins  _et al._ (2021) B. S. Simpkins, A. D. Dunkelberger, and J. C. Owrutsky, J. Phys. Chem. C 125, 19081 (2021).
  * Philbin  _et al._ (2022) J. P. Philbin, Y. Wang, P. Narang, and W. Dou, J. Phys. Chem. C 126, 14908 (2022).
  * White  _et al._ (2020) A. F. White, Y. Gao, A. J. Minnich, and G. K. L. Chan, J. Chem. Phys. 153, 224112 (2020).
  * Eisenschitz and London (1930) R. Eisenschitz and F. London, Zeitschrift für Phys. 60, 491 (1930).
  * London (1930) F. London, Zeitschrift für Phys. 63, 245 (1930).
  * Dahlke and Truhlar (2007) E. E. Dahlke and D. G. Truhlar, J. Chem. Theory Comput. 3, 46 (2007).
  * Barron (2017) J. T. Barron, Continuously differentiable exponential linear units (2017), arXiv:1704.07483 .
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014), arXiv:1412.6980 .
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Bussi and Parrinello (2007) G. Bussi and M. Parrinello, Phys. Rev. E 75, 056707 (2007).
  * Ceriotti  _et al._ (2009) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 103, 030603 (2009).
  * Ceriotti  _et al._ (2010a) M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010a).
  * Ceriotti  _et al._ (2011) M. Ceriotti, D. E. Manolopoulos, and M. Parrinello, J. Chem. Phys. 134, 084104 (2011).
  * Ceriotti  _et al._ (2010b) M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 133, 124104 (2010b).
  * Ceriotti  _et al._ (2014) M. Ceriotti, J. More, and D. E. Manolopoulos, Comput. Phys. Commun. 185, 1019 (2014).