Quantum-Electrodynamical Time-Dependent Density Functional Theory Description of Molecules in Optical Cavities

TL;DR

This paper introduces a new computational method to study molecules in optical cavities, showing how cavity effects can change molecular binding and structure.

Contribution

A new quantum-electrodynamical time-dependent DFT method (QED-TDDFT-TP) is introduced for modeling molecules in optical cavities.

Findings

QED-TDDFT-TP achieves good agreement with high-level methods for ground-state energies and polaritonic spectra.

Cavity confinement significantly alters binding energies and geometries of weakly bound dimers.

The method offers a computationally efficient and accurate tool for studying cavity-modified molecular interactions.

Abstract

We introduce a quantum-electrodynamical
time-dependent
density
functional theory with a tensor-product representation (QED-TDDFT-TP)
to model molecules strongly coupled to quantized cavity fields. By
combining real-space electronic wave functions with truncated Fock-space
photon states, the method captures light–matter correlations
at a computational cost close to standard DFT. Benchmark calculations
show good agreement with QED-FCI and QED-CASCI for ground-state energies
and polaritonic spectra. Applications to weakly bound dimersincluding
(H2)2, Ar2, (H2O)2, and HFdemonstrate that cavity confinement can significantly
alter binding energies and geometries in a polarization-dependent
manner. The framework provides an accurate and scalable tool for studying
cavity-modified molecular structure and interactions.

Figures

Comparison of LiH PES calculated using QED-DFT-TP with PN-QED-FCI. All graphs use λ = 0.05 and ω = 0.121. The PN-QED-FCI-N F calculation used N F = 2 and N F = 10, and the QED-DFT-TP results were calculated using N F = 4.

LiH: Ground-state energy vs bond length for ω = 0.121 for λ = 0.01, 0.02, 0.03, 0.04, and 0.05 (left). Ground-state energy vs bond length for ω = 0.03025, 0.0605, 0.121, 0.1815, and 0.242 while keeping λ = 0.05 (right).

Dependence of the occupation of photon space |1⟩ on λ and ω for different bond lengths.

Left: Dependence of the polaritonic energies on λ. The upper and lower curves show the energies of the upper and lower polaritons; the middle curve shows the energy dependence due to the diamagnetic term without coupling to the light. Right: Upper and lower polaritonic energies after subtracting the diamagnetic term’s effect.

Comparison of the polaritonic energies of BH3 using QED-TDDFT and QED-FCI. The QED-TDDFT graph was shifted upward by 0.15 eV to fit it into the same scale. The graph begins at λ = 0.02 because for lower coupling strengths we do not resolve a splitting in QED-TDDFT.

Potential energy surface for (H2)2: λ = 0.05, ω = 12.7 eV. Each H2 molecule is aligned along the x-axis with a H–H distance of 1.4 au, centered at z = ±R/2.

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for Ar2. To enable comparison within a single figure, the occupation probabilities corresponding to λ = 0.05 have been offset by adding 0.38. ω = 0.467 au was used in the calculations.

Potential energy surface for Ar2. λ = 0.1 and ω = 0.0375 au used in the calculations.

Two water dimer configurations used in the calculations.

Potential energy surfaces for (H2O)2 systems.

Two HF dimer configurations used in the calculations.

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for the HF dimer (dimer 1 in Figure ).

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for the HF dimer (dimer 2 in Figure ).

Full text

Introduction

Cavity quantum electrodynamics (cavity QED) is a fundamental framework in quantum physics that studies the interaction between light and matter in confined electromagnetic environments. Its importance spans both theoretical understanding and practical applications. Cavity QED provides precise control over light–matter interactions by confining photons in optical cavities alongside atoms or other quantum emitters. This confinement modifies the electromagnetic environment, leading to phenomena such as the Purcell effect,? where spontaneous emission rates are enhanced or suppressed depending on the cavity properties. These controlled interactions reveal fundamental aspects of quantum mechanics, including entanglement between light and matter ?−? ? ? and quantum superposition states. ?,? Light–matter coupling as a means to tune physical and chemical properties has become a major focus of experimental research. ?−? ? ? ? ? ? ? ? ? ? ? ? ? Theoretical investigations have also developed alongside the experimental progress. ?−? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Several excellent review articles survey the current state of experimental and theoretical approaches to cavity light–matter interactions. These encompass reviews on hybrid light–matter states, ?−? ? ?

ab initio computational methods, ?,?,? and molecular polaritonics. ?−? ? Describing coupled light–matter systems theoretically and computationally presents significant challenges. The quantum many-body problem involving electron–nuclear interactions is already complex, and incorporating photon degrees of freedom makes it even more demanding. Recent years have seen numerous approaches developed ?,?,?,?,?,?−? ? ? ? ? ? ? ? ? ? ? that extend beyond the basic two-level atom model.? These methods typically build upon established many-body quantum techniques, adapting them to account for photon interactions.

The Pauli–Fierz (PF) nonrelativistic quantum electrodynamics Hamiltonian has emerged as the most practical framework ?,?,?,?,? for computational applications. The Pauli–Fierz Hamiltonian is the fundamental theoretical framework for describing nonrelativistic quantum electrodynamics (QED), originally developed by Wolfgang Pauli and Markus Fierz in 1938. ?−? ? The PF Hamiltonian consists of three components: the electronic Hamiltonian, the photonic Hamiltonian, and an interaction term that couples electrons and photons. The presence of this coupling term necessitates the use of a combined electron–photon wave function, where electronic states are represented through an appropriate basis set while photonic states are expressed using the Fock-space representation. The Pauli–Fierz Hamiltonian is characterized by several distinctive features. First, it operates within a nonrelativistic framework, which differs from complete relativistic QED by treating matter particles nonrelativistically while preserving the quantum nature of the electromagnetic field. Second, it employs minimal coupling ?,? to describe light–matter interactions, achieved by substituting the canonical momentum p with p – e A, where A represents the electromagnetic vector potential. Third, many practical implementations, especially in cavity QED and molecular physics, utilize the long-wavelength or dipole approximation,? which considerably reduces computational complexity.

To solve the Pauli–Fierz Hamiltonian, one must include both the matter degrees of freedom and the quantized electromagnetic field. This is because the Hamiltonian contains products (couplings) of matter operators with photon-field operators, so the appropriate state space is the tensor product $[eqn]$ , where $[eqn]$ is the Hilbert space of the matter degrees of freedom and $[eqn]$ is the bosonic Fock space of photons. Concretely, a state is written as a sequence of n-photon components Ψ = (ψ^(0)^, ψ^(1)^, ψ^(2)^, ...), where each ψ^(n)^ depends on the particle coordinates and spin and on n photon variables (momenta and polarizations), and is symmetric under exchange of the photon labels.

Similar to the Schrödinger equation, the Pauli–Fierz Hamiltonian lacks analytical solutions for multielectron atoms. For a single-electron atom or ion, the problem becomes tractable by constructing a product basis from hydrogenic eigenfunctions and Fock basis states, allowing exact diagonalization to yield the solution.? However, systems with more than one electron require numerical methods for their solution.

As is typical in electronic structure calculations, methodologies can be broadly categorized into two distinct families: wave function–based methods and density–based approaches. Wave function–based methods ?,?−? ? ?,?,? characteristically employ coupled electron–photon wave functions, and their product structure leads to a substantial increase in computational dimensionality. The coupled electron–photon wave function can be written as

$$[eqn]$$

where Φ_ nm _ is a many-body basis function representing the electrons (and nuclei, if present), |n⟩ is a Fock-space basis for photons, and C _ nm _ is the linear expansion coefficient. The Fock-space basis can represent a single mode or multiple photon modes. The Φ_ nm _ notation emphasizes that the spatial basis functions can be different for different photon sectors if needed.

The simplest approach, the cavity QED Hartree–Fock, ?,? extends the traditional Hartree–Fock method to include quantized electromagnetic field modes within optical cavities. This approach treats the coupled electron–photon system using a mean-field approximation, where electrons experience an effective field created by all other electrons and the cavity photon modes. The method employs a polaritonic wave function ansatz that is typically written as a product of an electronic Slater determinant and a photon state. refs. ?,?,?,? employ a coupled-cluster (CC) methodology, ?,?,? that constructs a reference wave function from the direct product of a Hartree–Fock Slater determinant and the photon vacuum state. The ground-state QED-CC wave function is then defined by applying an exponentiated cluster operator to this product state. The primary advantage of this method lies in its systematic improvability. The traditional Complete Active Space Configuration Interaction (CASCI) approach has been also extended to include quantized electromagnetic field modes.? In the CASCI ansatz for the electronic subspace, a subset of active electrons and orbitals are identified, where a full CI expansion is performed within that active space. In QED-CASCI, this framework is generalized to simultaneously treat both electronic correlation within the active space and the coupling to cavity photon modes.

The stochastic variational method (QED-SVM) ?,?,?,? similarly employs a product form combining matter and photonic wave functions, but differs in its treatment of the matter component through explicitly correlated Gaussian basis states. Variational parameters are optimized via stochastic selection procedures, yielding highly precise energies and wave functions. Due to the N! scaling of explicit antisymmetrization of the N-particle basis functions, the practical application of the QED-SVM approach is restricted to small atomic and molecular systems.

The density-based approach, namely quantum electrodynamical density functional theory (QED-DFT), is an extension of traditional density functional theory (DFT). ?,?−? ? ? ? ? ? QED-DFT bridges the gap between quantum optics and electronic-structure theory, making it possible to describe phenomena where light and matter interact strongly, such as in optical cavities. QED-DFT is an exact reformulation of the PF Hamiltonian, based on many-body wave theory. In QED-DFT, the complex coupled electron–photon system is represented by two uncoupled, yet nonlinear, auxiliary quantum systems. The electrons are described by the usual DFT equation which now contains potentials describing the interaction of light and matter. A separate Maxwell-like equation is used for the photons. The QED-DFT calculations mostly use real-space bases but extensions to Gaussian basis representation also exist. ?,? Combinations of QED-DFT with macroscopic QED, ?,? and extensions to Dicke ?,? and Rabi models? have also been developed.

Standard electronic exchange-correlation (XC) functionals are inadequate for QED-DFT because they fail to account for electron-photon correlations that emerge under strong light-matter coupling. This limitation results in inaccurate predictions of polaritonic energy levels, ground-state modifications, and photon distribution statistics, making the creation of specialized QED exchange-correlation (QED-XC) potentials crucial for these systems. ?,?,? Recent theoretical advancesincluding optimized-effective-potential formulations and exact-model benchmarkshave established the formal structure of QED-XC and demonstrated the role of photon-mediated electron–electron interactions in modifying electronic structure and dynamics. ?,?−? ? These developments enable ab initio predictions of cavity-induced modifications to excitons, charge transfer, and chemical reactivity, successfully reproducing observed Rabi splittings and vacuum-field effects in molecules and materials. ?,?,?

Our QED-DFT methodology employs a coupled electron–photon wave function analogous to those utilized in wave function–based methods. This wave function is constructed on a tensor product combining a spatial grid with a Fock-state representation. To differentiate this framework from previously discussed QED-DFT approaches, we designate the current method as QED-DFT-TP and QED-TDDFT-TP. The QED-DFT-TP represents a specific implementation of QED-DFT that adopts an alternative ansatz through the use of a coupled electron–photon wave function. While the tensor product formulation elevates the computational dimensionality, it maintains the discrete nature of quantized photon states. The coupled electron–photon wave function offers an enhanced characterization of light–matter interactions through the calculation of spatial wave functions within individual photon sectors. In this approach, each molecular orbital is paired with distinct Fock basis states representing quantized photon modes. The light–matter interaction component of the Hamiltonian governs the coupling between orbital elements across various photon states. The orthogonality of Fock states maintains the sparse structure characteristic of real-space DFT Hamiltonians. This sparsity enables the implementation of computationally efficient iterative diagonalization techniques commonly employed in conventional real-space DFT methodologies.

This paper aims to use the QED-DFT-TP methodology for computing various physical properties of molecules within optical cavities, investigate the influence of cavity parameters on these properties, and benchmark the results against established theoretical approaches. Small molecules, including LiH, BH_3_, Ar_2_, H_2_, HF and water dimers, will be used as examples.

Formalism

The systems we consider in this paper are all nonrelativistic and as a result the light–matter coupling can be consistently described by the Pauli–Fierz nonrelativistic QED Hamiltonian. ?,?,? In addition, since we are working with small-sized systems, we assume that the spatial variation of the cavity field is negligible over the dimension of the system, i.e., we will use the dipole approximation. The PF Hamiltonian in the velocity gauge can be written as a sum of the kinetic energy, Kohn–Sham potential, and the photonic Hamiltonian,

$$[eqn]$$

where V KS(r) refers to the Kohn–Sham (KS) noninteracting potential adapted from the KS-TDDFT scheme? and is given by

$$[eqn]$$

where ρ is the electron density, V H is the Hartree potential, V XC is the exchange–correlation potential, and V ion is the external potential due to the ions. The exchange–correlation potential V XC is approximated using the generalized gradient approximation (GGA), developed by Perdew et al.?

In the long-wavelength limit, the vector potential is spatially uniform over the matter extent,

$$[eqn]$$

with polarization ε α, quantization volume V, and frequency ω_α_. The expansion of the kinetic term in (?) contains the paramagnetic coupling $[eqn]$ and the diamagnetic (seagull) term $[eqn]$ . By introducing

$$[eqn]$$

the paramagnetic term becomes

$$[eqn]$$

and the diamagnetic term

$$[eqn]$$

The paramagnetic interaction links photon states that differ by one quantum number (Δn = ±1), while the diamagnetic interaction connects photon states with quantum number changes of Δn = 0, ±2. In the special case where the diamagnetic term couples photon states with identical quantum numbers (Δn = 0), it behaves analogously to the dipole self-interaction (DSI) found in the length-gauge formulation, introduced below.

The Hamiltonian can also be transformed into the length gauge (see Appendix):

$$[eqn]$$

where D is the dipole moment and Ê α the transverse electric field of mode α. The first term of the second line couples photon states with Δn = ±1. The last term of the second line is the DSI, coupling only states with Δn = 0. This term is always present and, contrary to what was believed previously, plays a crucial role in the variational formulation of the eigenvalue problem.? The relation of length and the velocity-gauge Hamiltonians are further discussed in the Appendix.

The coupled system is described by orbitals defined on a tensor product of a real-space and a Fock-space. At the KS level, we can represent the orbitals as

$$[eqn]$$

where |n⟩ is the Fock-space basis for the photons, N F is the dimension of the Fock-space, and N occ is the number of orbitals. For this paper, we will assume that there is one dominant mode and we can ignore all others, i.e., N p = 1. Our system is, thus, described by a four-dimensional (4D) grid: N _ x _ × N _ y _ × N _ z _ × N F, where N _ x _, N _ y _, and N _ z _ are the number of grid points in Cartesian real space, and N F refers to the size of the truncated Fock-space, i.e., the vacuum state |0⟩ has N F = 1. Because the Fock-basis states are orthogonal, the elements of the overlap matrix are given by

$$[eqn]$$

where the round bracket stands for integration over both real and Fock-space and the angle bracket is integration over only the real part,

$$[eqn]$$

The calculation of the matrix elements can be simplified further by orthogonalizing the real part of the orbitals for each Fock-state using the Gram–Schmidt method. This new orthogonal set can then be normalized

$$[eqn]$$

where ϕ̂ _ mn _ with (m = 0, ..., N occ) represents the orthogonalized set of basis for the same photon state |n⟩. In the present work, the minimization for the coupled light–matter orbitals is carried out by the conjugate-gradient method. The construction of the Hamiltonian matrix in the coupled basis is described in detail in our previous work.?

The ground-state calculation follows conventional DFT approaches using steepest descent or conjugate gradient approaches to calculate the energies and the orbitals.

The ground state orbitals will be used to initialize the time propagation

$$[eqn]$$

Any time propagation method typically used for TDDFT can be used here, and in this work we use Taylor time propagation?

$$[eqn]$$

where Δt is chosen to be sufficiently small to conserve the norm of the orbitals during propagation.

To compute the absorption spectrum we apply a weak instantaneous perturbation (a “delta kick”) to the ground-state wave function,

$$[eqn]$$

where κ is small and x̂ is the dipole operator. The time-dependent Schrödinger equation is then propagated, and the dipole moment d(t) = ∑_ m ⟨Φ m (t)|x̂|Φ m _(t)⟩ is recorded. The absorption spectrum is obtained from the imaginary part of the Fourier transform of the dipole response,

$$[eqn]$$

Results

The primary objective of these calculations is to benchmark the QED-DFT-TP approach against other established methods across various molecules and molecular dimers. This comparison is particularly timely given the rapid development of numerous new approaches in this field. To ensure a correct comparison, we have employed identical parameters across all methods, including coupling strengths, cavity frequencies, and molecular geometries. We represent λ = λε, where ε is a unit vector describing the polarization of the cavity mode, e.g., (1, 0, 0), and λ is the coupling strength. In this work, we will use λ ≤ 0.1, which, according to $[eqn]$ , corresponds to subnm^3^ effective volumes. ?,? This range coincides with volumes achieved in picocavity experiments. ?,? The velocity-gauge Hamiltonian is employed in all calculations. The velocity gauge allows more straightforward analysis of photonic Fock state populations compared to most ab initio cavity QED approaches. The latter typically use the length gauge and therefore require Power–Zienau–Woolley transformation of the ladder operators to extract physical photonic Fock state populations. The differences of results in velocity and lengths gauges are discussed in Appendix A. Since this calculation employs a finite-difference grid to represent the wave functions, the computed total energies are sensitive to the alignment between grid points and ionic positions. To maintain consistency when investigating energy as a function of intermolecular distance, we preserve the relative positioning by ensuring that molecular coordinates remain commensurate with the underlying computational lattice. This constraint limits our ability to position molecules at arbitrary locations; instead, molecular displacements are restricted to shifts with integer multiples of the grid spacing.

Selected Test Cases

This work examines the following 6 molecular systems under cavity quantum electrodynamics conditions. We begin by revisiting the potential energy surface (PES) of LiH, which was previously explored in our earlier study.? Recent developments? now provide additional benchmarks for comparison with our methodology. The second system of interest is BH_3_ for which we calculate the Rabi splitting via time propagation of molecular orbitals to obtain the absorption spectrum, enabling comparison with the recent work of ref ?. Subsequently, we investigate the PES of an H_2_ dimer to characterize cavity-induced modifications to intermolecular forces, with results benchmarked against those reported in ref ?. A parallel study of the Ar dimer allows us to compare our approach with QED-DFT calculations from ref ?, providing a valuable cross-validation between distinct QED-DFT formulations. We then examine cavity effects on hydrogen bonding by computing the PES of a water dimer, comparing our findings with ref ?. Finally, we explore the orientational dependence of the cavity-modified PES for an HF dimer, a system not previously characterized in the literature. The molecular geometries and other parameters of the calculation are listed in the Supporting Information for each system.

LiH

In this section, we model the ground-state potential energy surface of LiH and compare our results with Photon-Number Quantum Electrodynamical Full Configuration Interaction (PN-QED-FCI) calculations from ref ?. The LiH molecule is oriented with its internuclear axis parallel to the polarization vector of a cavity mode with frequency ω = 0.121 au This frequency is chosen to be resonant with the molecule’s lowest singlet excitation.?

Figure presents a comparison of PN-QED-FCI, and QED-DFT-TP results. The QED-DFT-TP calculation fully converges using N F = 4. Only the lowest Fock states exhibit significant coupling to the electronic subsystem, with the |0⟩ state accounting for approximately 98% of the total photonic probability. To facilitate visual comparison, the QED-DFT-TP energies have been uniformly shifted by Δ = −7.264 hartree. All computational parameters were selected to match those employed in the PN-QED-FCI benchmarks of ref ?. The agreement between methodologies is excellent, with nearly identical predictions for the equilibrium bond length. The primary discrepancy emerges in the dissociation limit, which may stem from either the compact 6-311G basis set employed in the wave function methods or limitations of the GGA exchange-correlation functional in QED-DFT-TP.

Comparison of LiH PES calculated using QED-DFT-TP with PN-QED-FCI. All graphs use λ = 0.05 and ω = 0.121. The PN-QED-FCI-N F calculation used N F = 2 and N F = 10, and the QED-DFT-TP results were calculated using N F = 4.

We quantified the agreement between potential energy surfaces obtained with different light–matter electronic-structure methods using the nonparallelity error (NPE), defined as the difference between the maximum and minimum absolute energy deviations along a given reaction coordinate. This metric isolates differences in shape between two PESs, independent of constant vertical energy offsets, and therefore provides a stringent test of whether two approaches predict the same underlying physics across configuration space. For the LiH case, the calculated NPE between the PN-QED-FCI and QED-DFT-TP is approximately 0.3 eV, indicating that the curves are nearly parallel over the full bond-length range. The residual nonparallelity arises from a small curvature mismatch rather than from systematic distortions of the PES.

Figure presents the ground-state energy of LiH as a function of bond length for various coupling strengths λ and cavity frequencies ω. The molecular system is coupled to both the |0⟩ and |1⟩ photon Fock states.

LiH: Ground-state energy vs bond length for ω = 0.121 for λ = 0.01, 0.02, 0.03, 0.04, and 0.05 (left). Ground-state energy vs bond length for ω = 0.03025, 0.0605, 0.121, 0.1815, and 0.242 while keeping λ = 0.05 (right).

The results demonstrate that increasing λ elevates the ground-state energy due to the diamagnetic contribution in the Pauli–Fierz Hamiltonian. Notably, the overall shape of the potential energy curves remains qualitatively unchanged across different coupling strengths, exhibiting primarily a vertical shift that scales approximately as 2λ^2^. This behavior indicates that the diamagnetic self-energy contributes a nearly constant offset across all bond lengths.

In contrast, increasing ω reduces the ground-state energy (Figure), producing an opposite effect on the system’s energetics. This inverse relationship arises because the diamagnetic term, which dominates the energy modification, exhibits an inverse dependence on the cavity frequency. Figure illustrates the variation in |1⟩ state occupation, which quantifies the photonic excitation from the |0⟩ state. The bond-length dependence of the |1⟩ occupation follows a trend similar to that observed for the energy in Figure. The occupation increases substantially with larger λ values; note that the occupations in Figure are scaled by different multiplicative factors to facilitate visual comparison on a single plot. The cavity frequency dependence of the |1⟩ occupation (Figure) exhibits behavior similar to the energy variation for ω = 0.121 au and ω = 0.1815 au However, for larger ω values, the occupation decreases with increasing internuclear separation. This is likely due to the fact that increasing the frequency ω raises the energy of the nn nth photonic mode according to ℏω. This increased energy separation between states diminishes the effective coupling strength between the |0⟩ state and higher photon number states |n⟩, consequently reducing the occupation probability of these higher-lying states.

Dependence of the occupation of photon space |1⟩ on λ and ω for different bond lengths.

BH3

In this section, we compare our calculated Rabi splitting for BH_3_ with QED-FCI results from ref ?. While the QED-FCI approach determines the Rabi splitting by directly diagonalizing the Hamiltonian to obtain the energy separation between upper and lower polariton states, our methodology employs an alternative strategy. We compute the absorption spectrum via real-time propagation and extract the Rabi splitting from the spectral positions of the polaritonic peaks.

Our TDDFT time-propagation approach evaluates the absorption spectrum by applying a brief delta-function electric field pulse to the ground-state system, then propagating the time-dependent Kohn–Sham equations to track the subsequent evolution of the electronic density. The time-dependent induced dipole moment is recorded throughout the propagation, from which we obtain the frequency-dependent polarizability through Fourier transformation. The absorption spectrum is then derived from the imaginary component of the polarizability, which corresponds directly to the optical absorption cross-section. The Rabi splitting is identified as the frequency separation between the upper and lower polariton peaks in this spectrum.

The cavity-free case exhibits an excitation peak at 0.4732 au (see Figure), which closely matches the molecule’s third singlet excited state reported by Vu et al. (2024).? Initially, we couple the molecule to the |0⟩ Fock-state. Although this configuration lacks light–matter coupling, the excitation peak energy increases due to the positive diamagnetic term, as demonstrated in Figure. To determine the upper and lower peak positions, we set the cavity frequency ω equal to E/ℏ, where E represents the absorption peak energy (middle curve in Figure). Since the diamagnetic term causes the peak position to shift with varying λ, we correspondingly adjust ω in our calculations. When light–matter coupling is introduced, the single absorption peak splits into two distinct peaks. The energies of these split peaks define the upper and lower polaritons.

Left: Dependence of the polaritonic energies on λ. The upper and lower curves show the energies of the upper and lower polaritons; the middle curve shows the energy dependence due to the diamagnetic term without coupling to the light. Right: Upper and lower polaritonic energies after subtracting the diamagnetic term’s effect.

Figure illustrates the λ dependence beginning at λ = 0.015. Below this threshold value, peak splitting is not observed. This limitation arises because QED-TDDFT-TP produces absorption spectra with peaks of finite width, causing overlap at small λ values. While extending the propagation time could potentially resolve this issue, such calculations would be computationally impractical. As anticipated, the upper polariton energy increases with λ. However, the lower polariton energy exhibits minimal variation. This behavior occurs because the peak frequencies shift upward with increasing λ, creating the appearance that E LP remains constant. To provide a clearer visualization, the right panel of Figure shows the polariton energies after removing the diamagnetic contribution (by subtracting the middle curve from both the upper and lower curves). With this adjustment, both upper and lower polaritonic energies display the expected λ-dependent behavior.

Figure presents a comparison of polaritonic energies for BH_3_ calculated using QED-TDDFT and QED-FCI methods for the same excited state. The QED-FCI results are taken from Vu et al. (2024).? The graph reveals that while both methods demonstrate increasing relative energy with rising λ values, QED-TDDFT-TP consistently underestimates the splitting magnitude and exhibits greater nonlinear behavior compared to QED-FCI. Additionally, our approach underestimates the upper polariton energy, but the overall agreement between the QED-FCI and QED-TDDFT-TP is encouraging.

Comparison of the polaritonic energies of BH3 using QED-TDDFT and QED-FCI. The QED-TDDFT graph was shifted upward by 0.15 eV to fit it into the same scale. The graph begins at λ = 0.02 because for lower coupling strengths we do not resolve a splitting in QED-TDDFT.

(H2)2

Intermolecular forces can be fundamentally understood as electronic interactions mediated by transverse electromagnetic fields. ?,? Consequently, modifying electromagnetic boundary conditions through cavity confinement can substantially alter intermolecular interactions.? Moreover, in the strong coupling regime, molecules can interact with one another via delocalized cavity photons even at separations where direct Coulombic interactions become negligibly weak. Figure presents the potential energy surface of (H_2_)2 computed using QED-DFT-TP, compared with QED-FCI results from ref ?. The QED-DFT-TP method reproduces the overall topology of the QED-FCI surface and correctly predicts that cavity modes with ϵ_ z _ polarization (parallel to the intermolecular axis) yield lower binding energies than those with ϵ_ x _ polarization (perpendicular). However, QED-DFT-TP systematically overestimates the binding energy across all configurations. The comparison between QED-DFT and QED-FCI reveals substantially larger shape discrepancies. For x-polarized coupling the NPE is approximately 6 meV, while for z-polarization it increases slightly to about 6.5–7 meV. In both polarizations, the two methods coincide in the dissociation limit but differ markedly near the equilibrium region, where QED-DFT predicts a significantly deeper minimum than QED-FCI. Consequently, the NPE is dominated by the depth mismatch at short and intermediate separations, indicating that the PESs are not parallel. This behavior reflects an overbinding tendency of QED-DFT relative to the wave function reference in the cavity-modified van der Waals regime. For deeper understanding of the differences convergence studies of the FCI approach would be valuable.

Potential energy surface for (H2)2: λ = 0.05, ω = 12.7 eV. Each H2 molecule is aligned along the x-axis with a H–H distance of 1.4 au, centered at z = ±R/2.

A notable discrepancy is that QED-DFT-TP predicts the ϵ_ x _ polarization to produce a potential energy surface nearly identical to the cavity-free case, whereas QED-FCI shows substantial cavity modification for this polarization. It is important to recognize that in the cavity-free limit, our method reduces to conventional DFT, which cannot accurately describe van der Waals interactions despite yielding qualitatively reasonable curves.? We also note that the QED-DFT approach employed in ref ? does not appear to reproduce a bound potential energy minimum for this system.

Several factors may account for the stronger binding predicted by QED-DFT-TP relative to QED-FCI. First, there is a difference in reference energies: QED-DFT-TP employs a reference state with the two water clusters separated by 200 Å, whereas such large separations are computationally intractable for QED-FCI, necessitating a 30 Å reference separation instead. Additional discrepancies could arise from the approximate exchange-correlation functional used in the DFT approach or from basis set truncation effects in the FCI calculation.

Ar2

We have also investigated cavity-modified intermolecular interactions in Ar_2_. The two argon atoms are positioned along the z-axis with varying separation distance R. In the cavity-free case (λ = 0), our results closely match those of ref ?, predicting an equilibrium separation of 3.9 Å, although our binding energy is approximately 5 meV stronger. As shown in Figure, the binding energy remains nearly unchanged when the cavity polarization ϵ is perpendicular to the dimer axis and the coupling strength is weak (λ = 0.05). This behavior can be understood by examining the occupation probability of the |1⟩ Fock state (Figure), which exhibits minimal variation with internuclear separation, leaving the intermolecular potential essentially unmodified relative to the cavity-free case. Increasing the coupling strength to λ = 0.1 introduces a stronger distance dependence in the |1⟩ state occupation (Figure), resulting in reduced binding energy compared to the cavity-free reference.When the cavity polarization is parallel to the dimer axis, the occupation of the |1⟩ state exhibits pronounced sensitivity to the interatomic separation, and the binding energy decreases substantially with increasing λ.

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for Ar2. To enable comparison within a single figure, the occupation probabilities corresponding to λ = 0.05 have been offset by adding 0.38. ω = 0.467 au was used in the calculations.

Figure displays the Ar dimer energy as a function of interatomic distance using parameters λ = 0.1 and ω = 0.0375 atomic units, identical to those employed in ref ?. When compared to Figure, the perpendicular configuration exhibits enhanced binding energy that approaches the cavity-free binding energy more closely, which results from the 10-fold reduction in frequency. Conversely, the parallel configuration shows reduced binding relative to the previous case. ref ?. predicted that parallel configurations would be less bound than the cavity-free case, while perpendicular configurations would be more strongly bound than the cavity-free case. Our findings are consistent with ref ?. regarding the parallel configuration but disagree concerning the perpendicular configuration. This difference can be attributed to several methodological distinctions. While we employ GGA functionals, the alternative QED-DFT approach utilizes hybrid functionals. Additionally, our method incorporates photonic effects through a tensor product representation, whereas the other approach implements photon many-body dispersion within the Hamiltonian.

Potential energy surface for Ar2. λ = 0.1 and ω = 0.0375 au used in the calculations.

(H2O)2

We investigated the influence of cavity confinement on hydrogen bonding in a water dimer by varying the oxygen–oxygen separation distance RR R. This system was recently examined in ref ?. Since the exact atomic coordinates from that study were not available, we constructed two geometrically similar configurations based on the structures depicted in ref ?, as shown in Figure. The cavity polarization is oriented along the O–O axis, and the cavity frequency and coupling strength are set to match those of ref ?: ω = 7.86 eV and λ = 0.1.

Two water dimer configurations used in the calculations.

Figure shows that our potential energy curve exhibits excellent agreement with the Coupled-Cluster Singles and Doubles (CCSD) results from ref ? in the cavity-free case. Both CCSD and QED-DFT-TP predict an equilibrium separation of 2.9 Å. The binding energies are also in good agreement: CCSD yields 214 meV (Figure 6 in ref ?), while QED-DFT-TP predicts 240 meV. Under cavity confinement, CCSD calculations indicate a modest 30 meV reduction in binding energy, whereas QED-DFT-TP predicts a substantially larger decrease of approximately 130 meV. This more pronounced modification is physically reasonable given the strong light–matter coupling regime accessed by the chosen parameters. The quantitative discrepancy between the two approaches likely stems from differences in the water dimer geometries employed. This interpretation is supported by the substantial variation in cavity-induced binding energy changes observed between our two distinct water dimer configurations, which underscores the geometric sensitivity of cavity effects on hydrogen bonding.

Potential energy surfaces for (H2O)2 systems.

HF Dimer

We examined the distance-dependent binding energy of HF dimers in two distinct molecular arrangements: parallel orientation (HF–HF) and antiparallel orientation (HF–FH), as illustrated in Figure. The HF molecular axes are oriented along the x-direction. Two cavity polarization configurations were investigated: ϵ = (1, 0, 0), designated as “perpendicular” since it is orthogonal to the intermolecular z-axis, and ϵ = (0, 0, 1), termed “parallel” as it aligns with the axis connecting the two molecules.

Two HF dimer configurations used in the calculations.

For the parallel HF–HF configuration, Figure demonstrates that the cavity-free system exhibits no intermolecular binding, and light–matter coupling in either polarization fails to induce molecular association. The parallel polarization yields higher total energies than the perpendicular case, while both polarizations produce minimal |1⟩ photon-state occupation across all intermolecular separations.

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for the HF dimer (dimer 1 in Figure ).

In contrast, the antiparallel HF–FH arrangement (Figure) exhibits substantial intermolecular binding with an energy minimum at approximately 2.7Ånotably about three times the HF monomer bond length of 0.92 Å. Light–matter coupling reduces the binding strength for both polarization directions relative to the cavity-free case. Parallel polarization decreases the equilibrium separation, while perpendicular polarization increases it. The |1⟩ photon-state occupation in the antiparallel configuration is significantly enhanced compared to the parallel arrangement. The occupation profiles exhibit polarization-dependent behaviors: at short intermolecular distances they differ substantially but converge at large separations. Notably, the parallel polarization exhibits an occupation minimum near 3 Å, whereas the perpendicular polarization shows monotonic behavior with no such feature across the investigated distance range.

Potential energy surface (left) and occupation probability of the |1⟩ Fock-space (right) for the HF dimer (dimer 2 in Figure ).

Summary

This work examines the effectiveness of combining quantum electrodynamics with time-dependent density functional theory to model molecular systems under strong coupling conditions within optical cavities, benchmarking the results against existing theoretical approaches. Conventional density functional theory is extended by incorporating the Pauli–Fierz nonrelativistic quantum electrodynamics Hamiltonian, with the coupled electron–photon system represented on a tensor product of real-space and Fock-space. The method is benchmarked against established wave function approaches including photon-number quantum electrodynamical full configuration interaction and complete active space configuration interaction for small molecular systems. Potential energy surfaces for LiH are shown to demonstrate excellent agreement with reference calculations, while Rabi splitting calculations for BH_3_ exhibit qualitative agreement despite some quantitative differences. The influence of cavity parameters on intermolecular interactions is investigated through studies of hydrogen-bonded and van der Waals dimers, including (H_2_)2, Ar_2_, (H_2_O)2, and (HF)2 systems. Cavity coupling is found to significantly modify binding energies, equilibrium distances, and photon occupation numbers, with strong dependence on molecular orientation relative to field polarization observed. For HF dimers, parallel (HF–HF) configurations are shown to remain unbound in cavities while antiparallel (HF–FH) arrangements exhibit substantial binding that decreases upon light–matter coupling.

The QED-TDDFT-TP approach demonstrates computational tractability while maintaining sufficient accuracy to describe cavity QED effects and polariton physics, reproducing the qualitative features of strong coupling regimes. Residual quantitative discrepancies relative to wave function calculations warrant further analysis of finite photon-basis effects and density-functional approximation quality.

The QED-DFT-TP method presented in the paper is orders of magnitude cheaper than both QED-CC and especially QED-FCI. It retains explicit quantized photon states while keeping computational cost close to standard DFT, enabling simulations of realistic molecules and dimers under strong light–matter coupling.

The present work assumes a bare matter exchange correlation functional which is the same in each Fock-space. This approximation could be a source of error relative to the exact solution. Future work is needed to develop suitable exchange-correlation functionals for the coupled electron orbitals Fock space approach.

Supplementary Material