Numerically Exact Solution for a Real Polaritonic System under Vibrational Strong Coupling in Thermodynamic Equilibrium: Loss of Light–Matter Entanglement and Enhanced Fluctuations

The first numerically exact simulation of a full ab initio molecular quantum system (HD+) under strong ro-vibrational coupling to a quantized optical cavity mode in thermal equilibrium is presented. Theoretical challenges in describing strongly coupled systems of mixed quantum statistics (bosons and Fermions) are discussed and circumvented by the specific choice of our molecular system. Our numerically exact simulations highlight the absence of zero temperature for the strongly coupled matter and light subsystems, due to cavity-induced noncanonical conditions. Furthermore, we explore the temperature dependency of light–matter quantum entanglement, which emerges for the ground state but is quickly lost already in the deep cryogenic regime. This is in contrast to predictions from the Jaynes–Cummings model, which is the standard starting point to model collective strong-coupling chemistry phenomenologically. Moreover, we find that the fluctuations of matter remain modified by the quantum nature of the thermal and vacuum-field fluctuations for significant temperatures, e.g., at ambient conditions. These observations (loss of entanglement and coupling to quantum fluctuations) have implications for the understanding and control of polaritonic chemistry and materials science, since a semiclassical theoretical description of light–matter interaction becomes reasonable, but the typical (classical) canonical equilibrium assumption for the nuclear subsystem remains violated. This opens the door for quantum fluctuation-induced stochastic resonance phenomena under vibrational strong coupling, which have been suggested as a plausible theoretical mechanism to explain the experimentally observed resonance phenomena in the absence of periodic driving that has not yet been fully understood.


S1 Numerical Solution of Pauli-Fierz Hamiltonian and Ensemble Averaging
The exact solution of our quantized 3-body matter system coupled to a cavity mode relies on the choice of a highly optimized, i.e. problem specific, coordinate representation, S1,S2 in combination with Gauss-Laguerre numerical quadrature (see Supporting Information of Ref. S3 for numerical details and S4 of this work for specific computational parameter choices).
When evaluating ensemble averages numerically, we assume V → ∞, which implies a continuum of quantum numbers k z , which allows an efficient numerical evaluation, i.e. we replace ∞ dk z .Hence, we can apply the Gauss-Hermite quadrature procedure to approximate the infinite integral along k z numerically, i.e. by using the following relation, S4 where x i are the roots of the Hermite polynomial H s (x) and the weights w In our case, we map the ensemble averages of observable Â as, and accordingly, ) This assumes that the transformed Â does not depend on k z , k y , Rc and Pc , which is indeed the case for all our observables.The COM reduced Energy eigenvalue is defined as it is assumed that the ground-state contribution pα 0 to the ensemble average dominates.However, from the symmetry argument in Eq. ( 19) we immediately notice µ = 0, since pα = 0 for all states in absence of external driving currents, which makes accurate integration possible for few discrete k z -evaluations only.

S2 Numerics for Logarithmic Negativity
The canonical density matrix given in Eq. ( 9) can be written for our system in terms of COM, relative and photon basis explicitly as, which accounts for the blockdiagonal nature in terms of COM coordinates.For our subsytem choices W ∈ r ci , qα the partial transpose follows from Eq. (S6) by either c * i,s c j,t → c * j,s c i,t or c * i,s c j,t → c * i,t c j,s , leading e.g. to the partial transpose of the photon subsystem, which can be diagonalised yielding the corresponding eigenvalues w and eigenfunctions e −ik R c ∑ i,s d Γ q α i,s (k z , l) |r ci i |q α s .For a fixed k z -value, this diagonalization can efficiently be computed numerically for our HD + molecule and the choice of our subsystems.
The logarithmic negativity entanglement measure then follows immediatly from: where in the last step the summation over k were approximated by an integral, as it was previously the case for the ensemble averages.
S3 Jaynes-Cummings Light-Matter Entanglement in Canonical Equilibrium The three energetically lowest eigenfunctions of the JC-model are with partial transpose The four eigenvalues of the matrix ρ Γ m ,JC can be determined and we find three strictly positive eigenvalues and only one negative, which then solely enters the logarithmic negativity measure.

S4 Simulation Parameters
For the numerical solution, parameters for HD + in a cavity were taken from Ref. S3 (e.g.particle mass or radial grid scaling) if not explicitly mentioned otherwise.The number of radial grid points, was slightly reduced from N matter = 12 to N matter = 10, since we required additional resources to explicitly account for the ensemble averaging in k z (i.e.thermal effects acting on the charged s = 9 (subsystem temperatures and entanglement) grid points, which effectively required either 3 or 5 evaluations of the computationally expensive h(k z ) expression only.This simplification arises to the ensemble symmetry of k z with respect to the origin k z = 0, (i.e.µ = 0).The number of Fock states to represent the quantized field was chosen to be N pt = 4 for the fluctuation analysis, N pt = 3 for the subsystem temperature and N pt = 2 for the entanglement measurements.For our chosen coupling strengths a choice N pt > 2 only becomes relevant if one is interested in small deviations of the field fluctuations.Numerical convergence was ensured for all our results.
Figure S1: Light-matter entanglement measured by logarithmic negativity η with respect to the temperature T for a coupling constant λ = 0.005 at different cavity frequencies ω c (given in atomic units).On the top, the cavity frequency is tuned close to the first ro-vibrational excitation, whereas in the bottom plot the IR to visible regime is covered.In both cases, light-matter entanglement is lost quickly in the deep cryogenic regime.However, while the numerical results for the upper case are certainly converged, deviations may occur in the lower figure due to the highly optimized grid-representation, which was designed to reproduce groundstate and ro-vibrational matter properties with high accuracy, but not simultaneously (!) the higher electronic and vibrational excitations.These states cannot be populated anyways in the chosen temperature regime, but still the accuracy of the matter basis contributions, which are mixed into the hybrid groundstate, may still be reduced in principle.Nevertheless, the loss of light-matter entanglement at low cryogenic temperatures seems to be a generic property for thermal ensembles under strong coupling conditions, roughly independently of the chosen cavity frequency.This is in line with typical experimental evidence, which does not allow for quantum computing devices at sizable temperatures.
) which are composed from the bare matter {|g , |e } and photon {|0 , |1 } eigenstates, assuming the rotating wave approximation.Furthermore, we have assumed for simplicity that the cavity is tuned exactly on resonance with the first bare matter excitation energy ω α = E e − E g .Expressing the corresponding canonical equilibrium density in the light-matter basis {|g ⊗ |0 , |g ⊗ |1 , |e ⊗ |0 , |e ⊗ |1 } leads to the following matrix representation,