Chemistry in Quantum Cavities: Exact Results, the Impact of Thermal Velocities, and Modified Dissociation

In recent years tremendous progress in the field of light–matter interactions has unveiled that strong coupling to the modes of an optical cavity can alter chemistry even at room temperature. Despite these impressive advances, many fundamental questions of chemistry in cavities remain unanswered. This is also due to a lack of exact results that can be used to validate and benchmark approximate approaches. In this work we provide such reference calculations from exact diagonalization of the Pauli–Fierz Hamiltonian in the long-wavelength limit with an effective cavity mode. This allows us to investigate the reliability of the ubiquitous Jaynes–Cummings model not only for electronic but also for the case of ro-vibrational transitions. We demonstrate how the commonly ignored thermal velocity of charged molecular systems can influence chemical properties while leaving the spectra invariant. Furthermore, we show the emergence of new bound polaritonic states beyond the dissociation energy limit.


.1 Center of Mass (COM) Separation
In order to tackle the coupled light-matter problem defined by the Pauli-Fierz Hamiltoninan, given in Eq. (1) of the letter, in dipole approximation, it is useful to switch to a centre-ofmass (COM) coordinate system. This means that we define r i = R c + r ci , with the COM explicitly given as R c := i m i r i i m i . The Hamiltonian can then be re-written in the COM frame as,Ĥ (S1) Note that a priori the dipole operator still contains a COM dependence at this stage. In a next step, we apply the unitary Power-Zienau-Woolley (PZW) transformation shifting the COM,Û P ZW := e i Q tot λα·Rc ωαp α , with Q tot := N i=1 eZ i . By using the Baker-Campbell-Hausdorf formula, one can show that our PZW transformation obeys the following properties, S2 which allow to write the shifted Hamiltonian as, by applying the PZW transformation for each mode α. Afterwards, the original eigenvalue problem Hψ = Eψ can be simplified to Eq. (2) given in the letter by using a wave function Ansatz of the form ψ (R c , r c , q α ) = e ikRc Φ (r c , q α ). For one mode α and three bodies, one obtains the following expression, In case of neutral systems (i.e. Q tot = 0), the additional interaction of the COM motion with the quantized field vanishes. Note that in Eq. (S9) it was used that for one mode (i.e. M = 1) the PZW-shifted eigenvalue problem can be additionally simplified by absorbing 1 with the corresponding momentâ p = i ω 2 (â † −â) and position operatorsq = 1 2ω (â † +â). They obey the usual canonical commutation relations [q ,p ] = i. This shift will only be relevant for relatively light charged particles at low resonance frequencies ω (e.g. fundamental L0-L1 transition of HD+).

Observables
When calculating observables of the coupled system (λ = 0), a priori one cannot neglect any involved coordinate. However, in practise not always all integrals have to be solved explicitly. For example, the dipole oscillatory strengths can be calculated from COM relative S3 coordinates r ci as follows: where in the last step it was used that ψ|R In contrast, photonic observables have to be transformed back to the length gauge to be consistent. Hence, an integration over the COM position has to be performed explicitly.
However, due to the cylindrically symmetric setup with respect to the z-axis of the lab frame, i.e.
the COM integration is reduced to one dimension only.
The general expectation-value integral in our chosen spherical-cylindrical coordinate system (see Sec. 1.3 below) is given as, where in the last step it was assumed thatÔ does not explicitly depend on the COM coordinate and positions. Note, if λ = 0 or if Q tot = 0 or ifÔ does not depend onq and P cz , S4 the R cz integral is unity as it was the case for a matter only system.

Spatial Coordinate Representation
In order to make an exact solution of a quantized three-body system coupled to one cavity mode feasible, a smart choice of the spatial representation of the matter is pivotal. For this purpose we follow the approach in Refs. 1,2, which is known to reach excellent numerical accuracy for matter only calculations of three bodies, with only moderate computational costs. This means, we represent the COM relative coordinates r ci of our three bodies in a 6-dimensional combined spherical and cylindrical coordinate system. Its angular coordinates {θ, φ, ψ} will eventually be treated analytically, whereas the radial coordinates {R, ρ, ζ} will be treated numerically after an additional transformation into perimetric coordinates.
To obtain the combined spherical and cylindrical coordinates, one expresses the COM relative coordinates r ci in terms of two vectors X, Y, which are chosen such that r c2 −r c1 = X, Moving the origin of X, Y to the COM, results in, In a next step, one expresses X in spherical coordinates and Y in cylindrical coordinates S5 with respect to X. The resulting expressions are, 2 Using Eqs. (S15)-(S17), the COM relative coordinates and the corresponding transformed volume element becomes A sketch of the combined coordinate system is given in Fig. S1. Figure S1: Schematics of the coordinate system. We assume the relevant (green) cavity mode polarized along the z-direction, i.e. λ α z. The spherical coordinates of vector X are shown in red, whereas the cylindrical coordinates of vector Y are displayed in blue.

Basis Set Expansion
First of all, λ-coupling is assumed along the z-axis only, i.e.
In a next step, we employ the Ansatz wave function defined in Eq.
The resulting representation of the Pauli-Fierz Hamiltonian takes the following block-diagonal form: where the first term corresponds to the matter-only problem promoted to the coupled matter- where β, γ, are defined in spherical-cylindrical coordinates as, The coefficients contain the analytical evaluation of the angular integrals of the λ 2 -term, S11 which amount to the following non-zero values, For the λ angular integrals, the resulting ± √ 3 3 was already included in H SP 0 and H SP 1 , respectively. Analysing the matrix given in Eq. (S29) in terms of S, P even and P odd , one notices that the block-diagonal nature of the non-interacting terms remains preserved by the λ 2 -term, only broken by mixing of S and P odd states due to the λ-term. Note that one can show that S-states do not mix via the λ-term for any excited angular momentum states beyond l = 1. However, this is not necessarily true for λ 2 contributions.

Gauss-Laguerre Quadrature for Radial Integrals
So far it was only stated that there is a matrix representation of the coupled Hamiltonian, but it was not yet specified how to treat the radial coordinates R, ζ, φ numerically. For this purpose, a coordinate transformation into a h i -scaled perimetric coordinate system of the following form is performed in a first step, where with new volume element 2 dV = h 1 h 2 h 3 sin(θ)(x +ỹ)(x +z)(ỹ +z)dR cx dR cy dR cz dxdỹdzdφdθdψ, andx := h 1 x,ỹ := h 2 y,z := h 3 z. The scaling factors h i will later be used to adjust the radial grid to the spatial extend of simulated syste. In a next step, the orthonormal basis given in Eqs. (S25)-(S28) is rewritten as, where a regularization factor R(x, y, z) = ρR = √ xyz(x+y+z) 2 was introduced in agreement with the literature. 3 It suppresses singularities of the matter-only Hamiltonian, which may cause numerical difficulties. However, its only practical relevance is restricted to radial momentum operators, which do not appear in the coupling Hamiltonian and thus R eventually cancels. The newly introduced scaled Lagrange function F ijk (x,ỹ,z) is defined as, with The Lagrange-Laguerre functions are defined as, with L N the Laguerre polynomial of degree N with roots u i and Lagrange property f i (u j ) = (λ N i ) −1/2 δ ij . The coefficients λ N i can be chosen to fulfill the Gauss-Laguerre quadrature by using the orthonormality property of F ijk for the perimetric volume element in Gauss- S14 N ijk = 1 is implied from the normalisation condition. Hence, the original eigenvalue problem given in Eq. (S9) is now discretized and numerically accessible by solving for E and c.

Convergence and Numerical Tests
Multiple explicit convergence/ sanity checks were perform to ensure that finite basis set errors or implementation mistakes do not spoil the results. Matter-only energy eigenvalues were compared with reference calculations from literature given in Tab. S1. For He, the λ 2 -scaling of the ground-state could be compared with QEDFT calculations with photon OEP accuracy in the Born-Oppenheimer limit, which indicates an agreement on the same accuracy level as one expects from the previous matter-only considerations (see Fig. S1). All QEDFT simulations were performed with the OCTOPUS code. 6 Note that highly accurate simulations of ground-state nuclear contributions in QEDFT would be very challenging to obtain. 7 Last but not least, the Thomas-Reiche-Kuhn sum rule was well preserved for He, HD+ and H + 2 , which implies that there is no fundamental implementation error present. Moreover, all observables were consistent with theory and particularly in agreement with the JC model.
This offers an additional sanity check of the numerical results (see main section of the manuscript). S16 -0.58715567914 -0.585948 2P (odd) * -0.4990065652928 * -0.498039 * Figure S2: Cutting the basis set expansion for angular momentum quantum numbers l > 1, may introduce significant numerical errors for stronger couplings. In order to check the validity by allowing e.g. l = 2, the code complexity of the implementation would be more than doubled. Therefore, we decided to use an alternative route. The comparison of the the He ground-state energy shift with results from QEDFT simulations with photon OEP 9 indicates that inaccuracies from l < 2, are on the order of milihartree or below, which is in line with the accuracy reached for the absolute matter-only energies given in Tab. S1. As it is the case for matter-only values, one expects considerable smaller relative errors in terms of energy differences.

HD+
In Fig. S4 a zoom of the dispersion relations given in the main section is shown to visualize dressing effects, caused by dipole self-interaction, and the shift of the resonance frequency due to the non-zero net-charge. In Fig. S5a Figure S5: In (a) we consider the HD+ resonant case for ∆E at λ = 0.01 with respect to the kinetic energy E kin of the COM. While the spectrum is not changed (up to numerical inaccuracies for higher-lying states) the COM motion (a) redistributes the oscillator strengths. Notice that the grey area indicates less reliable eigenvalues, which are not converged for the chosen photon number basis with N pt = 6. In (b) and (c) the HD+ polaritonic dispersion curves for k z = 0 and k z = 1 are shown with respect to the wave-function overlap between our exact calculations and the corresponding JC model. Second, bound excitations are invariant with respect to changes of the radial basis set and the corresponding scaling parameters, whereas the discrete continuum reacts very sensitively.
Based on these considerations, we could distinguish the bound states, which are shown in the main section, from continuum states composed of a H-atom and a free proton coupled to the cavity. With the latter method one can also identify the dissociation energy limits, i.e. when increasing the scaling factors of the radial grids one observes an accumulation of eigenvalues at the dissociation energies, 2 while the newly discovered bound polaritonic states remain invariant. Similarly to the main section, in Fig. S7 the proton-proton distance vs. energy plot is shown for the previously identified bound states. However, for this figure we assumed an infinite mass for the nuclei. Qualitatively, the system behaves very similar compared with the results for finite proton masses. Overall, the proton-proton distances of the bound states are slightly reduced and some of the excitation energy differences are moderately shifted compared with the finite mass reference. Nevertheless, the appearance of bound states below and above the dissociation limit remains preserved in the infinite-mass limit.
Aside from that, one can also investigate the bound 1S-iS transitions below the dissoci-S22 ation limit with respect to different COM velocities. Our simulations for H + 2 confirm that a finite COM motion of a charged molecule indeed leads to an increased dipole-transition probability, as it was discussed for HD+ in the main section.
(a) (b) Figure S6: Dispersion relations for the H + 2 molecule with singlet nuclear spin configuration in a cavity. The frequencies are centered around the dissociation energy of the bare matter system. It is assumed that √ ω λ = const, i.e. the coupling strength ∝ ω, and λ = 0.057 at resonance with the dissociation energy. The oscillator strength color bar is chosen to visualize changes arising from finite COM velocities in a relatively weak regime. The COM motion was set to E kin = 0 for (a), whereas for (b) a non-vanishing E kin = 0.37 [eV] was chosen, which is still in the non-relativistic limit, i.e. k ≈ 0.072c. The vertical dotted line indicates the dissociation resonance condition. Figure S7: Quantized (i.e. bound) proton-proton distances for H + 2 with respect to groundstate energy differences ∆E and corresponding oscillator strength in the infinite proton mass limit. The blue dots correspond to dressed bare matter states whereas red dots indicate the emerging bright photon replicas, which are absent without a cavity. The cavity frequency is ω = 2.15 eV (dashed vertical line) with λ = 0.051 and zero COM motion. The green area indicates energy ranges beyond the p + -dissociation limit according to matter-only simulations.