Polaritonic response theory for exact and approximate wave functions

Polaritonic chemistry is an interdisciplinary emerging field that presents several challenges and opportunities in chemistry, physics, and engineering. A systematic review of polaritonic response theory is presented, following a chemical perspective based on molecular response theory. We provide the reader with a general strategy for developing response theory for ab initio cavity quantum electrodynamics (QED) methods and critically emphasize details that still need clarification and require cooperation between the physical and chemistry communities. We show that several well-established results can be applied to strong coupling light-matter systems, leading to novel perspectives on the computation of matter and photonic properties. The application of the Pauli-Fierz Hamiltonian to polaritons is discussed, focusing on the effects of describing operators in different mathematical representations. We thoroughly examine the most common approximations employed in ab initio QED, such as the dipole approximation. We introduce the polaritonic response equations for recently developed ab initio QED Hartree-Fock and QED coupled cluster methods. The discussion focuses on the similarities and differences from standard quantum chemistry methods, providing practical equations for computing the polaritonic properties.


INTRODUCTION
2][3] In 1963, Jaynes and Cummings described with a phenomenological quantum picture the formation of polaritons from the interaction of photons and atomic states. 4Their predictions were later confirmed experimentally for Rydberg atoms in a quantum microcavity. 5,6The use of quantum optical devices to confine the fields, such as Fabry-Pérot cavities or metal nanostructures, [7][8][9][10][11][12] increase the coupling between matter and light, and polaritons have been observed since the 1970s for inorganic semiconductors 13,14 and from the 1990s for organic molecules in optical cavities. 15,169][20] Several experiments have now shown that electromagnetic confinement leads to essential modifications of processes such as chemical reactions, [17][18][19][20][21][22] singlet fission, [23][24][25] intersystem crossing, [26][27][28] crystallization and assembly, 29,30 as well as optical properties like absorption, scattering and emission, [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] although the reproducibility of the chemical reactivity modifications is not always straightforward. 47Crucial experimental aspects are still to be clarified, and cooperation between experimentalists and theoreticians can be of great assistance in unravelling the conditions and the consequences of vibrational and electronic strong coupling.9][50][51][52][53][54][55] Following well-established electronic structure a) Electronic mail: enrico.ronca@unipg.itb) Electronic mail: henrik.koch@ntnu.notheory, several ab initio quantum electrodynamics (QED) approaches have been developed recently.These include QED density functional theory (QEDFT), 56 reduced density matrix approaches for matter-photon systems, 57,58 polaritonic coupled cluster, 59 QED Hartree-Fock (QED-HF), 60,61 QED coupled cluster (QED-CC), 60,61 QED full configuration interaction (QED-FCI), 60,61 polarized Fock states approach (PFSs), 62 polaritonic unitary coupled cluster (QED-UCC), 63 strong coupling (SC)-QED-HF, 64 QED Møller-Plesset perturbation theory of second-order (QED-MP2), 65 and second-order QED algebraic diagrammatic construction scheme for the polarization propagator (QED-ADC(2)). 655][86][87][88][89][90][91][92][93] While several of these results are still valid for QED systems, essential and subtle differences arise due to the explicit modelling of the photon field, highlighting additional challenges and opportunities.Although some response schemes for QED methods have been proposed, [94][95][96][97][98][99][100][101] a systematic discussion of polaritonic response theory and its connections to molecular response theory is still not available.This paper aims to fill this gap by proposing a systematic discussion of wave function based methods for computing polaritonic properties.We will review several results from molecular response theory that can be applied to polaritonic wave functions and critically discuss the limitations and future challenges such methods have to tackle.In section (2), we briefly introduce the most common optical devices employed to achieve strong coupling.In section (3), we present the Pauli-Fierz Hamiltonian and focus on how this framework is used to model the strong coupling regime.We review the most common approximations for ab initio QED in section (4).In section (5), we discuss polaritonic response theory for exact states, focusing on similarities and differences with molecular response theory, and provide several examples of polaritonic response properties.The derived equations are equivalent to QED-FCI response theory.In section (6), we discuss how to develop response theory for approximate models and derive linear and quadratic response equations for QED-HF and QED-CC.Finally, section (7) recaps the main results described in this paper and contains some concluding remarks and future perspectives.

OPTICAL CAVITIES
This section provides a brief introduction to quantum optical devices.The discussion, although far from being exhaustive, aims at reviewing the most common experimental setups employed in polaritonic chemistry and provides the reader with phenomenological basics of light-matter strong coupling.
The use of quantum optical devices provides a non-invasive way to engineer a material's properties by exploiting the confinement of the electromagnetic field.The geometry and the material of the device define the resonator eigenmodes.In a real cavity, the confinement is always defective and depends specifically on the frequency of the field and the cavity material.The so-called Q-factor measures the quality of a cavity 102 Q(ω 0 ) = ω 0 energy stored power loss , where ω 0 is the resonance frequency of the device.This equation implies that the electric field will exhibit an exponential decay such that the energy spectrum will show a Lorenzian lineshape ( From Eq. ( 3), we have a more practical definition of the Qfactor as the ratio between the frequency of the hosted field mode ω 0 and the full width at half maximum (FWHM) ∆ω of the cavity spectrum The coupling λ between matter and the confined cavity modes is proportional to the inverse square root of the effective confinement volume The efficiency of a device in enhancing the interaction between light and matter can then be measured by the Q/V e f f ratio.Strong coupling conditions are obtained once the coupling strength exceeds the material and cavity losses (i.e. the coherent energy exchange rate between the electromagnetic field and the molecule exceeds the dissipation processes).Decoherence interactions with the environment are suppressed, and hybrid states between the cavity and the molecular eigenstates naturally describe the system.In this regime, the molecular properties are therefore affected even in the absence of an external photon pumping of the device.
4][105][106][107][108] These resonators are often produced by sputtering gold or silver on a substrate (often SiO 2 , BaF 2 or ZnSe), and direct contact with the sample is prevented by using a thin polymer film.The sample is then coated by this multilayer structure or injected into the hollow cavity.The cavity eigenmodes are linearly-polarized standing waves whose wavelength λ n depends on the distance L between the mirrors The spatial distribution of the optical mode can be observed by changing the position of a thin sample slab inside the cavity.The strong coupling features of this system are a function of the slab position, which gives us the field strength distribution. 46,1091][112][113][114][115][116] The coupling strength λ in Fabry-Pérot cavities is limited by the diffraction limit, which restrains the field confinement.Devices which bound the fields on a scale smaller than the diffraction limit are called subwavelength cavities and can facilitate ultrastrong coupling.To this end, using nanostructured metamaterials in a suitable geometric arrangement can allow efficient confinement and overcome the diffraction limit, leading to efficient subwavelength cavities.In Fig.

PAULI-FIERZ HAMILTONIAN
The starting point of molecular modelling is the nonrelativistic molecular Hamiltonian In Eq. ( 7), the capital indices M and N refer to nuclei of mass and charge m M and Z N , i and j to electrons with mass m = 1 and charge q = −1.The electronic and nuclear positions are labelled r and R respectively, and p is the linear momentum of the particles.The first line of Eq. ( 7) contains the nonrelativistic kinetic energy operators for nuclei and electrons, while the second line describes the Coulomb interaction.
Although it is often assumed that the Hamiltonian in Eq. ( 7) is sufficient to model a molecular system, this is not the case.This is apparent if we realise that the excited states of the system have an infinite lifetime, i.e. they do not radiate.In Eq. ( 7), we disregard the electromagnetic degrees of freedom, retaining only the electrostatic interaction between particles while neglecting retardation effects.The motion of electric charges generates radiation fields which in turn influence the particles' motion: this back-reaction mechanism is the origin of spontaneous decay.The inclusion of the electromagnetic degrees of freedom is also fundamental when optical devices are employed to modify the photon environment.3][164] The particle motion is determined by the Lorenz force where c is the speed of light, v and a are the velocity and the acceleration of the particle, q is the charge and m is the mass.The electric E and magnetic B fields are determined by Maxwell's equations, here expressed in cgs units where r is a vector defining a point in space, ρ(r,t) is the charge density and j(r,t) is the current density.The fields are naturally described using potentials (φ , A), such that We notice that, unlike the fields B and E, the vector potential A and the scalar potential φ are not measurable quantities and must, therefore, only be seen as mathematical tools.In Eqs.(10) and (11), the potentials are identified only through their derivatives.This implies that the same fields, i.e. the same system, can be modelled using different potentials.The transformation known as gauge transformation, leads to unchanged electromagnetic fields for any well-behaved scalar function χ(r,t).This gauge freedom can be exploited to simplify the equations.We will make use of the so-called Coulomb gauge which allows us to account for the Coulomb interactions among the particles explicitly.The scalar potential is, in this gauge, determined by the Poisson equation as in the electrostatic case while the vector potential is obtained from the following equation: where j ⊥ is the solenoidal (divergence-free) part of the density current.
The Hamiltonian for the field-particle system is described in terms of the position of the particles, the vector potential and their conjugate momenta.The self-consistent interaction between molecules and light is described by the Pauli-Fierz Hamiltonian, for which Hamilton's equations can be shown to be equivalent to Eqs. ( 8) and ( 9) 162,163 In the first line of Eq. ( 17), we have the vector potential A computed at the positions of nuclei and electrons, and the second line describes the Coulomb interactions between the particles.The explicit separation of Coulomb interactions provided by the Coulomb gauge allows for the separation of Hamiltonian in Eq. ( 7), providing a natural chemical description of the system.The last line accounts for the electromagnetic field energy.Expanding the first terms, we obtain an interaction term between the particle Hamiltonian in Eq. ( 7) and the radiation fields mediated by the vector potential A where x runs on both electrons and nuclei.Following the QED prescription, the vector potential in Eq. ( 17) is promoted to an operator where b † kτ and b kτ respectively create and annihilate a photon with wave vector k and frequency ω k = |k| c .The field polarization kτ with τ = 1, 2 spans the 2D plane perpendicular to k.If we take the polarization kτ to be real, the fields are described by a superposition of linearly polarized plane waves.In an equivalent way, we can take the complex orthogonal unit vectors which describe left and right circularly polarized waves.

Remarks
There are now some subtle points to mention.Choosing the Coulomb gauge, the vector potential only has (two) transverse components A = A ⊥ .At the same time, the longitudinal component of the electric field (Coulomb field) is described by the scalar potential.In this gauge, the equations are not manifestly covariant.However, this does not imply a loss of relativistic invariance since the electromagnetic field prediction agrees with special relativity.In fact, the retardation of the interaction between particles is obtained by a cancellation of the instantaneous interactions from the Coulomb term and the transverse field.For a different gauge, the Coulomb interactions would also be mediated by the longitudinal component of the vector potential. 162,163he canonical momentum p i , contrary to the standard quantum mechanical formulation, is not the linear momentum of the particle but here has a field-dependent component 162,163 This means that, although in Eq. ( 17) we can identify the original particle Hamiltonian in Eq. ( 7), the physical meaning of these two operators is different.The Pauli-Fierz Hamiltonian is self-adjoint and bounded from below 165,166 such that a stable ground state exists.][169][170] In Eq. ( 17), we have disregarded the nuclear and electronic spin, which could be introduced in Eq. ( 17) with the addition of the following contribution where the first term of Eq. ( 22) accounts for the interaction of the electronic spin, described by the Pauli matrices σ i , with the (internal) magnetic field B, while the second term describes the energy associated with the nuclear spin I N with magnetogyric ratio γ N .Finally, in Eq. ( 17) particles can in principle interact with infinitely high-frequency modes.However, the description of the kinetic energy is nonrelativistic.When the field frequency is comparable with the rest energy of the particles ω ∼ mc 2 , relativistic effects (e.g.creation of electron-positron pairs) not included in Eq. 17 become relevant. 162It is then necessary to introduce a cutoff as the field's momentum k increases (ultraviolet cutoff).Moreover, the masses in Eq. ( 17) are not the physical masses usually employed in quantum mechanics, which include the contribution of the electromagnetic energy created by the particle, but they are their bare masses. 163

APPROXIMATIONS FOR AB INITIO POLARITONIC CHEMISTRY
While the Hamiltonian in Eq. ( 17) allows for a consistent treatment of light and matter, approximations are needed to perform computational studies.The problem of the continuum of photonic modes could be overcome by employing a fine discretization of the spectrum.Moreover, different approaches can be developed depending on the treatment of the photonic degrees of freedom.In the Cavity Born-Oppenheimer approach (CBO), [171][172][173] the photon coordinates q α are embedded in the nuclear wave function χ(R, q,t) and separated from the electronic degrees of freedom, described by and electronic wave function φ (r; R, q).The complete wave function can then be expanded as In a polaritonic approach, the electronic and photon coordinates are treated on the same footing and described by a polaritonic wave function which depends parametrically on the nuclear coordinates only 56,60,81,82,[174][175][176] While these Born-Huang expansions are in principle equivalent, approximations will lead to different results.In the simplest case, we take into account only a single term of these expansions as in standard electronic Born-Oppenheimer approximation.In the discussion of exact response theory in Section (5) we do not refer to any explicit form of the Hamiltonian.In Section (6), we focus on the polaritonic approach and describe the QED-CC and QED-HF response theory.8][179] Nevertheless, we will not explore such effects, and we will work in a fixednuclei framework.Although these approximations already reduce the complexity of the problem, further simplifications are needed to deal with the electromagnetic environment and treat the molecule and its interaction with the fields with chemical accuracy.

Cavity QED
The Hamiltonian in Eq. ( 17) allows for simultaneous treatment of light and matter, where we could, in principle, include any optical device that contributes to a modification of the electromagnetic environment.Nevertheless, the molecular complexity and the presence of the optical apparatus make the problem impractical from a computational point of view.
The explicit modelling of the photonic device is avoided by defining an effective change in the structure of the electromagnetic modes.To this end, we impose boundary conditions on the fields (e.g.vanish on the surface of the optical device) to model the electromagnetic confinement induced by the cavity.The vector potential is expanded in terms of a complete orthonormal set of suitable electromagnetic normal mode functions {S α (r)}. 102,180For a rectangular parallelepiped with perfectly conducting walls of sides L x , L y and L z , the transverse mode functions are standing waves 180 where V c = L x × L y × L z is the total volume of the paral-lelepiped and k is the wave vector which is now quantized: In the absence of external charges, if we assume that the walls of the parallelepiped are grounded, the scalar potential vanishes identically φ = 0.The quantization is then performed by promoting the expansion coefficients to quantum operators, and the radiation Hamiltonian for the free fields inside the perfect optical device is discretized in terms of the cavity eigenmodes where b † α is the creation operator for the α-mode of frequency ω α .
In the presence of molecules, which can be considered as a source of electromagnetic fields, we must also ensure that the scalar potential in Eq. ( 15) is consistent with the boundary conditions.We can enforce them by adding an auxiliary scalar field F to the particles' Coulomb potential.This potential can be thought of as a potential generated by image charges 102,180 placed outside the volume of quantization (i.e.outside the volume of interest of the system).This potential can be expressed in terms of a response kernel Θ(s; r) and the real charge distribution ρ inside the cavity 102,180 where s is a point outside the cavity volume.The Hamiltonian describing the molecular system then reads 180 As seen from ( 27), the last two terms introduced by F correspond to a modified Coulomb interaction kernel due to the boundary conditions, i.e. they lead to a modified longitudinal interaction among the particles.The factor 1/2 appears because it describes an interaction with an image charge distribution. 102While some authors suggest that this contribution is the major difference from free space, [181][182][183] the discussion of these terms for ab initio QED is often neglected.Nevertheless, we will see that they can be handled in a more practical way by a suitable unitary transformation. 72,180

Remarks
The Hamiltonian in Eq. ( 28) now acts on a Hilbert space V which is the direct product of the full-CI Fock space for the electrons V FCI and the photon space spanned by the cavity eigenmodes However, each cavity setup corresponds to specific boundary requirements, leading to different quantum fields: different boundary conditions could substantially modify both the shape of the eigenmodes S α (r) and the Coulomb kernel Θ(s; r).The photonic space V C is usually truncated, considering only a few eigenmodes relevant to the system.Therefore, the excited states of the Hamiltonian in Eq. ( 28) are again true eigenstates with an infinite lifetime.Notice that, in the limit of an infinitely large cavity, if we retain a considerable number of modes we obtain again a fine energy discretization of the photon bath that recovers the decay channels. 48he masses of the particles in (28) are dependent on the electromagnetic mode structure and, in general, are different from their physical and free-space bare masses. 163,184However, in QED ab initio calculations, the physical masses are usually employed, implicitly considering the interaction energy with the continuum photon bath 48 .The effects of the modification of the masses on chemical and physical properties from an ab initio perspective, to the best of our knowledge, is still to explore.Finally, because of the states of the optical device V C , the system's symmetries now differ from the bare molecule ones due to the cavity's modified eigenmodes and potential in Eq. (27).This leads to modifications of selection rules and other symmetry-related properties.

Dipole approximation
The dipole approximation is commonly employed in ab initio polaritonic chemistry.It assumes that the (relevant) elec-tromagnetic modes have a wavelength much larger than the characteristic lengths of the molecules.The complete shape of the modes S α (r) is therefore irrelevant, and we only need to evaluate the fields at a point internal to the molecular struc-ture.In Eq. ( 17), we then set A(r) → A(0), and obtain the Pauli-Fierz Hamiltonian H v PF in dipole approximation and velocity representation where we applied the BO approximation.The dipole approximation allows for the modelling of the photons only by adequately tuning the frequency ω α and the coupling strength λ α for the cavity modes. 48The problem is eventually reduced to a single effective cavity mode.Assuming to work with real field polarization vectors, the final Hamiltonian for our system is where λ α = λ α α is the light-matter coupling vector of the α electromagnetic mode, and we used standard second quantization notation for the electronic Hamiltonian. 185This formulation is suitable for nonchiral cavities, while for helicitypreserving devices, we necessarily need to work with complex polarization vectors (see Eqs (19) and ( 20)).The Hamiltonian in Eq. ( 31) can be recast in a more convenient form, called length representation H l PF , by employing a unitary transformation U such that which will allow us to write the interaction between the molecule and the fields in a manifestly dipolar fashion.The transformation is where d is the electronic dipole operator.By using Eq. ( 33) and a transformation which changes the phases of the creation/annihilation operators 162,186-188 we finally obtain The bilinear light-matter interaction term of this Hamiltonian can be interpreted as the interaction of the molecular dipole with the displacement field. 162The quadratic term is called dipole self-energy and ensures the Hamiltonian leads to welldefined states and properties. 186Moreover, transformation Eq. ( 33) cancels the image charge contribution from F so that no reference to the image charges appears in Eq. ( 35). 72,163,180otice that these Hamiltonians are also often referred to as length or velocity gauge, although in a QED framework this can be misleading, and the term form or representation is more appropriate.The hierarchy of approximations that are employed to study quantum light-matter systems is pictorially summarized in Fig. 2.

Remarks
There are several advantages to employing the Pauli-Fierz Hamiltonian in the length representation.Contrary to the velocity representation Eq. ( 30), Eq. ( 35) has a light-matter interaction term linear in the field with no reference to the image charge distribution.The transformation in Eq. ( 33) implies that the momentum operator p now represents the kinetic momentum of particles so that the first line of Eq. ( 35) truly represents the Coulomb and kinetic energy of the matter subsystem only.
The operator b † is no longer a purely photonic operator since U in Eq. ( 33) mixes electromagnetic and matter degrees of freedom.The original photon creation operator in the length representation after the tranformatoin is 162,186,188 VUb These operators are now connected to the auxiliary fields of the macroscopic Maxwell's equations, where we can recognize the transverse polarization 162,186 The separation into "matter" and "photon" degrees of freedom is, therefore, blurred in this representation. 162

Beyond the dipole approximation
Ab initio QED calculations often rely upon the dipole approximation, which provides a good compromise between accuracy and affordability.0][191] To go beyond the dipole approximation, we could replace the vector potential A(r) with its first-order expansion around the molecular origin.For a sinusoidal field, this reduces to an expansion in terms of the wavevector k 3][194] The corresponding QED Hamiltonian can then be obtained in the velocity representation by the following substitution in Hamiltonian Eq. ( 28) The multipolar Hamiltonian can be obtained similarly to the length representation of the dipole Hamiltonian, applying the Power-Zienau-Woolley (PZW) transformation with the multipolar expansion of the vector potential, 163 i.e. the transformation of Eq. ( 33) with the vector potential of Eq. ( 40).This transformation also leads to electric (dipole and quadrupole) and magnetic self-energy terms, which guarantee the boundedness of the Hamiltonian. 163,186,190Overcoming the dipole approximation is particularly relevant when the molecules interact with chiral fields, which leads to novel phenomena such as cavity-induced circular dichroism or cavity enantiomeric discrimination. 64,190The description of magnetic interactions is then essential.However, the multipolar expansion of the interaction beyond the electric dipole presents some difficulties, already in the semiclassical approximation, [192][193][194][195] as the interaction terms are now origin dependent.The dipole op- 2. Graphic summary of the hierarchy of approximations for computational polaritonic chemistry.The starting point is the open system described by the nonrelativistic Pauli-Fierz Hamiltonian, which includes the molecular system and the optical device immersed in the photon continuum.As a first approximation, the device is assumed to be perfect, and cavity QED formalism with a limited number of effective photon states is employed.Then, in the dipole approximation, the field is assumed to be uniform over the molecular scale.The exact shape of the photon states is irrelevant, and the system is described through an effective coupling.9][60][61][62][63][64][65] Alternatively, there are phenomenological models where the molecular complexity is usually reduced to a few selected reference states, as, for instance, in the Jaynes-Cummings model. 4ator also depends on the choice of the origin for charged systems, but this can be shown to be equivalent to a gauge transformation on Hamiltonian Eq. (31).Manifest origin invariance can then be recovered by a suitable coherent state transformation, even for a finite electronic basis set (see Eqs. ( 201) and ( 202)). 60,61,64At the same time, to the best of our knowledge, a similar solution has yet to be developed for higher-order Hamiltonians, and therefore using such an expansion could lead to unphysical results.This issue could be solved by retaining the complete shape of the field S α (r) instead of relying upon a multipolar expansion.This is considered computationally and theoretically challenging also because of the dependence of the field shape on the boundary conditions.[201]

EXACT POLARITONIC RESPONSE THEORY
In this section, we highlight the similarities and differences with standard molecular response theory following the derivation of Jørgensen and Olsen, 84 with particular emphasis on the quantities that describe the cross-talk between light and matter.
From the exact solution of the time-dependent Schrödinger equation we derive, for an operator Ω, the Ehrenfest theorem 202 where H is here a polaritonic Hamiltonian and |Ψ are polaritonic states of entangled light-matter character.The operator Ω can then be an electronic, photonic or mixed electronphoton operator.The exact eigenstates of the Hamiltonian are defined by the eigenvalue equation and are assumed to form a complete basis for the radiation-matter Hilbert space V .We further assume that at the initial time the system is in its ground polaritonic state |0 .We then apply the following perturbation operator where η is a positive infinitesimal that ensures the perturbation to vanish for t → −∞.The time evolution of the initial ground state is then determined by the time-dependent Hamiltonian and can be modelled by a unitary transformation where we used a tilde to indicate the time-evolved state and P (t) is a Hermitian operator: The ground state phase evolution can be explicitly factorized 85 and if the operator Ω in Eq. ( 42) does not involve time differentiation, it can be disregarded.We will then make use of the phase-isolated wave function and perform an expansion of the state transfer coefficients in orders of the perturbation n + P (2) which leads to the perturbative expansion of the phase-isolated state j P * (1) j + P (1) n ∑ j>0 P (1) By defining P −n ≡ P * n and by using in Eq. ( 42) the state transfer operators Ω ∈ |0 n| ; |n 0| , we obtain the following hierarchy of differential equations, for positive and negative indices 84 n + iθ (kn)V t n P (1) where θ (x) is the Heaviside step function, while ω k , D k and V t n are defined as

D
(3) (1) − j (59) Notice that the solutions for lower-order coefficients are necessary for higher-order equations.We can solve these equations, which all have the structure with the general solution which fulfils the vanishing initial conditions for the coefficients.Expressing the change of the time development average value of an observable A from its zero-order mean value in orders of the perturbation we obtain We then define the linear, quadratic and cubic response functions, respectively These functions express the linear, quadratic and cubic time-variation of an observable A as resulting from a perturbation in the frequency domain and can therefore be interpreted as molecular (hyper)polarizabilities.For instance, for an electric dipole perturbation the first order variation of the molecular dipole reads and the linear dipole-dipole response function d; d ω can then be interpreted as the (negative) time-dependent molecular polarizability at frequency ω.Using the perturbation expansions of the coefficients and the state wave function Eqs. ( 50) and ( 51)-( 54), we obtain By comparing Eqs. ( 64) and ( 73), using Eq. ( 63) for the state transfer parameters, we finally obtain the explicit expressions for the linear, quadratic and cubic response functions in terms of the eigenstates of the Hamiltonian 84 where P(ω 1 , ω 2 ) and P(ω 1 , ω 2 , ω 3 ) sum all the permutations of the frequencies and we set the switching parameter η to zero.
These equations are identical to the expressions obtained in molecular response theory, [84][85][86][87][88][89][90][91][92][93] since they only rely on the Schrödinger equation.However, the explicit treatment of internal electromagnetic fields as dynamical variables in our system leads to novel perspectives.First of all, we must remember that the eigenstates of the system are mixed matter-photon states, so they belong to a larger Hilbert space than the usual molecular one.The properties of these field-dressed states will therefore differ from the bare-molecule ones, leading, for instance, to different excitation energies and transition moments.Moreover, additional possibilities for probing the properties of the system arise since the external perturbation V t can now act on both matter and photon degrees of freedom.Furthermore, considering the internal electromagnetic fields allows us to explore additional observables connected to the radiation fields and their interactions with matter.
The linear response function Eq. ( 74) describes the first-order variation of an observable A and it is connected to the linear polarizabilities of the system.Its poles occur at the polaritonic excitation energies, and the corresponding residues are connected to the transition moments between the ground and the excited states lim The quadratic response function Eq. ( 75) is connected to the first hyperpolarizabilities of the system since it describes the second-order response to an external perturbation.Its residues are connected to the transition moments and matrices between the excited states.These can be obtained once the transition moments between ground and excited states have been computed lim lim These expressions also give us access to excited state properties such as the dipole moment.Finally, the cubic response function describes the third-order variations of an observable as a consequence of an external perturbation, and it is therefore connected to the second hyperpolarizabilities of the system.
The response functions exhibit the following symmetry relations: for Hermitian operators A, B, C, and D.

Time-independent limit
The static response theory can be obtained from the previous derivation by setting the external frequencies ω and the switching parameter η to zero.Since the perturbation is now time-independent, it is possible to relate the static polaritonic properties to the derivatives of the Hamiltonian eigenvalues.The perturbation V acting on the system is time-independent and supposed to be smooth and small, and it is usually expressed as an external field F acting on a molecular multipole x The perturbed Hamiltonian then reads The Hamiltonian in Eq. ( 87) is endowed with well-defined eigenstates and eigenvalues, whose small changes compared to the unperturbed eigenvalues and eigenfunctions are commonly studied via the Rayleigh-Schrödinger (RS) perturbation theory. 203From the eigenfunctions |n of the unperturbed Hamiltonian H 0 (see Eq. ( 43)), the ground-state eigenfunction | 0 and the eigenvalue Ẽ of the perturbed operator are expressed as a power series: The first-and second-order corrections to the energy read: Eq. ( 90) is the well-known Helmann-Feynman theorem 204,205 and only requires the knowledge of the unperturbed wave function, while Eq. ( 91) also includes a sum-over-states contribution from each excited state |k with energy E k .The derivatives of the energy with respect to the external fields provide us with the ground state molecular properties.From Eqs. ( 86) and ( 90), we then identify the derivative of the energy with respect to the field F as the ground state mean multipole value From Eq. ( 91), the second derivative of the energy is then interpreted as the multipole polarizability to the field F , and from Eq. ( 74) we identify the linear response function for ω = 0 The second derivatives of the energy are then connected to the molecular polarizabilities.Higher-order derivatives relate hyperpolarizabilities and nonlinear response functions.Notice that these expressions are the same as standard molecular static response theory, and they also hold for exact polaritonic states since they are based only on the Schrödinger equation.

On the rotational average of computed molecular properties
In this section, we discuss the problem of the orientational average of the computed polarizabilities, which is straightforwardly accomplished in standard molecular response theory to connect the single-molecule calculation to the macroscopic sample response.Since we are now considering states that belong to the extended light-matter Hilbert space in Eq. ( 29), we have to consider the molecule and the electromagnetic environment, which are defined by the optical device.In general, this introduces anisotropy in the system.The cavity field can also break the molecular symmetry, possibly allowing for otherwise symmetry-forbidden transitions.This anisotropy-symmetry breaking introduces further complications in the computation of molecular properties.
The bare molecular states are independent of the spatial orientation of the molecule.The energy is invariant for rigid-body rotations of the system, while the charge density rotates following the molecular structure, and so do properties such as molecular polarizabilities.However, the strong coupling between molecular and cavity states depends on the projection of the molecular dipole operator onto the cavity field, as can be seen in the electron-photon interaction term in Hamiltonian Eq. (31).This also means that there is a non-trivial energy dependence on the relative molecular orientation in the cav-ity field.Excited states will be particularly affected by this FIG.3. Potential energy surfaces for a hydrogen molecule's upper and lower polaritons resonantly coupled with a quantum cavity, as a function of the molecular orientation.The energies are computed by using the time-dependent quantum electrodynamic Hartree-Fock method with a cc-pVDZ basis set and coupling constant set to λ = 0.01 a.u.(see section ( 6)).When the molecule's transition dipole moment is aligned with the cavity field (90°), the Rabi splitting is maximal, while it decreases to zero when the transition moment is perpendicular to the field's polarization (180°).orientation dependence.A bare electronic excitation with a transition moment orthogonal to the cavity field will have a negligible direct coupling with the photonic states.On the other hand, if the transition dipole moment is aligned with the field polarization, there will be mixing, giving rise to intense polaritonic excitations of hybrid light-matter character.This is illustrated in Fig. (3), where we plot the upper and lower polaritonic energies of a hydrogen molecule as a function of the angle between the transition dipole moment and the cavity field.These considerations on the molecular orientation pose some issues for the computation of molecular properties, as the direction of the cavity field is fixed in space by the experimental setup.In contrast, inside the cavity, the molecules are usually randomly or quasi-randomly oriented.Following a chemical approach, we could compute the sample's properties using the frequency-dependent polarizability Eq. ( 74), which must be averaged over the different molecular orientations 193,206 .The molecular orientation will be defined by a general set of parameters Ω.For instance, it can be uniquely defined by the three Euler's angles (φ , θ , χ).Notice that a different set of parameter can be convenient depending on the setup.For instance, for a Fabry-Pérot resonator with two degenerate modes with wavevector k and perpendicular polarizations, a rotation of the molecule along k does not change the Hamiltonian.Then it is sufficient to use only two parameters instead of the three Euler's angles.If we consider the Boltzmann weight of each different orientation, we have where A ; B ω is the averaged response function, E 0 (Ω) is the ground state energy of the polaritonic system as a function of the molecular orientation, k B is the Boltzmann constant, and T is the temperature.Even if we disregard the Boltzmann weight and perform an isotropic average, the result does not simply lead to an isotropic polarization tensor.The excitation energies in the denominator depend on the orientation, as well as the numerator.The polarization will be generally anisotropic because of the introduction of the optical device.
As a concrete example, we can consider the computation of absorption spectra.The linear absorption spectrum is connected to the dipole-dipole polarizability, where A and B in Eq. ( 74) are components of the dipole operator d i 84,85,206 Neglecting the dependence of the polaritons on the molecular orientation, we obtain the usual expression for absorption by randomly oriented molecules in terms of the oscillator strength This corresponds to fixing the relative molecule-field orientation.A fictitious peak-broadening is often applied to the computed spectrum S(ω), e.g. a Lorentzian lineshape with bandwidth broadening parameter ∆ as reported in Refs. 94,95so that the computed spectrum is However, different orientations could play an important role in determining the properties of the sample, and the choice of a fixed relative orientation does not give the full picture.For instance, if we are interested in a specific transition, we could suppose the field polarization to be parallel to the transition dipole moment.However, this could lead to biased results as certain effects could be enhanced while others suppressed.
The study of orientational disorder in ab initio polaritonic chemistry is still to be addressed.A possible approach would be to perform molecular dynamics simulations that sample the different orientations of the molecule inside the cavity and perform an average of the spectra obtained for each snapshot.This corresponds to an incoherent superposition of the spectra of an ensemble.The appearance of collective effects can introduce further complications in the computation of molecular properties.If we could perform an ab initio simulation including a sufficiently large number of molecules randomly or quasi-randomly oriented, the corresponding computed properties would include both collective and orientational effects.Alternatively, collective effects can initially be addressed with a simplified model, such as the Tavis-Cumming (TC) model, 67,68 which considers an ensemble of identical two-level non-interacting systems.Moreover, it disregards the dipole self-energy, the molecular dipoles and the so-called counter-rotating terms of the Hamiltonian.The Rabi splitting is then predicted to scale with the square root of the number of molecules and, therefore, it is sufficient to employ a smaller coupling strength compared to single-molecule calculations.The generalization of the TC model to an energy-broadened set of two-level molecules 207 predicts a Rabi splitting dependent on the quadratic average of the coupling strengths and on the energy broadening.The predicted spectrum would be very different from the one obtained by an average of single-molecule absorption spectra.
9][210][211] We believe both these effects are important.However, which one dominates will likely depend on the experimental conditions.Both should be carefully considered when accurately modelling polaritonic properties, and further studies in this direction are needed.

Equivalent expressions for polaritonic properties
In this section, we derive equations of motion for polaritonic response functions.Inserting the expansion Eq. ( 64) into Ehrenfest theorem Eq. ( 42), we obtain a set of equations as in standard molecular response theory 84 Considering Eq. ( 98) and the position-momentum operator identity which holds for the dipole Hamiltonian in the length representation Eq. ( 35), we obtain the relation which states that the frequency-dependent dipole polarizability can be evaluated regarding the position or the conjugate momentum.This ensures the equivalence of expressions in the dipole and velocity form of transition dipole moments where |n is an excited state and ω n is the excitation energy.We note that for the velocity representation, in Eq. ( 30), the analogous relation is which states, exactly as Eq. ( 102), an equivalence between the transition dipole moment and the transition kinetic momentum of the electrons.Indeed, Eqs. ( 100) and ( 103) are connected by the transformation in Eq. (33).We recall that in second quantization, Eqs. ( 100) and ( 103) are fulfilled only in the limit of a complete basis set. 185 example of a relation introduced by the electromagnetic degrees of freedom is provided by the commutation relation between the Hamiltonian and the annihilation operators in the length form which through Eq. ( 98) leads to for any operator B. We obtain from the residues of Eq. ( 105) the relation which allows the computation of photonic transition moments in terms of matter transition moments.This reflects the entanglement between electronic and photonic degrees of freedom.
The analogous relation for the creation operators is Moreover, from the relation we obtain If we introduce the conjugate field displacement q α and momentum p α coordinates we can write Eq. ( 109) as which is the equivalent relation of Eq. ( 102) for the photonic conjugate momenta.Similarly, by employing the Hamiltonian in the velocity form in Eq. ( 30) or through the transformation in Eq. ( 33), we obtain equivalent relations in the velocity representation.Notice that relation Eq. ( 109) also holds for the dipole Hamiltonian in the velocity form, but the physical meaning of the operators ( 111) and ( 112) is changed since they refer to a different representation.This will be discussed extensively in the following sections.
Analogous relations can be obtained from the higher-order equations of motion in Eqs.(99), for instance for the transition dipole moments among excited states we have Notice that relations similar to Eq. ( 114) also hold for the photonic conjugate momenta:

Response functions in cavity QED
This section provides a discussion of response functions in cavity QED.Although far from exhaustive, we provide the reader with several examples of matter-photon and photonphoton response functions.We focus on the peculiarities introduced by the photon dressing, and we explicitly discuss the result of using different mathematical representations of the Hamiltonian.
The perturbation is described by a single frequency component of Eq. ( 44) where this also includes the case of static perturbations by imposing ω = η = 0.

External electromagnetic fields
Spectroscopic techniques probe the system by means of an external electromagnetic field, whose electric and magnetic components are coupled to the motion of particles.The external probe is treated as classical (non-quantized) fields, described by its own vector A e (r,t) and scalar φ e (r,t) potentials.It is not necessary to describe the internal and external fields through the same gauge.The Hamiltonian of the system in the Coulomb gauge and Born-Oppenheimer approximation reads 162 where α labels the photon modes and i refers to electrons.Expanding the first term, we get the polaritonic Hamiltonian in Eq. ( 17) and the interaction term V t with the external fields: While the last three terms of Eq. ( 118) are also found in molecular response theory in the semi-classical approximation, 192,194,195,206 the first term is a purely QED contribution.This term cancels the field contribution from the p terms and ensures the coupling is only to the matter subsystem.When we perform the dipole approximation for the cavity fields and apply the transformation in Eq. ( 33), we obtain the length Hamiltonian H l .The first term in Eq. ( 118) is cancelled, leaving only the familiar interaction terms.By a suitable expansion of the external potential, the interaction term is obtained in a multipolar fashion as for the standard semiclassical theory [192][193][194]206 where we defined respectively the electric dipole d, magnetic dipole m, electric quadrupole Θ αβ and magnetic susceptibility χ αβ operators (in atomic units): Note that since we employ the dipole approximation in the length form, these multipolar operators refer exclusively to the matter subsystem and have the same physical meaning as in standard response theory.On the contrary, in the velocity form the operators has a mixed matter-photon character because of the field component in the conjugate momentum p (see Eq. ( 21)).
If we only retain the dipole interaction term in (119), the external electric field E e (0) computed at the molecular position is a constant that can be factorized out of the response functions.The dipole-dipole response function d i ; d j ω can then be interpreted as the negative molecular polarizability at frequency ω where i and j refer to cartesian components of the dipole operator d and k labels excited polaritonic states.This function describes the first-order time evolution of the molecular dipole moment subject to an external electric field: The poles occur when the frequency ω matches the energy difference between the ground and excited states, and the cor- provide information on the transition dipole moments between the ground and the excited states.While these expressions are the same as in standard molecular response theory, 84 the states involved are here polaritonic and the predicted properties will consequently be modified.
If we replace one of the electric dipole operators with the magnetic dipole m we obtain the electric dipole-magnetic dipole response function which describes the electric dipole response when a magnetic field is applied From the residues of the frequency-dependent electric dipolemagnetic dipole polarizability in Eq. ( 130), we obtain the rotational strength of optically active molecules lim As previously discussed, the operator m can only in the length representation be interpreted as the molecular magnetic dipole operator, while it has a different physical meaning in the velocity Hamiltonian.Notice that in polaritonic systems, the optical activity can appear both from molecular and field chirality.As discussed in Sec. ( 2), several groups have recently developed chiral cavities, 115,[147][148][149][150][151][152][153][154][155][156][157][158][159][160][161] which can host only one-handedness of light polarization within themselves.This opens up the possibility of engineering chiral properties through the chirality of the photon field, 36,37,147,[212][213][214] which is transferred to molecular states via light-matter strong coupling.6][217][218] On the other hand, chiral molecules interacting with a chiral field should give rise to "polaritonic diastereoisomers", possibly leading to novel approaches to chiral discrimination.
As the states of the system here belong to the radiation-matter Hilbert space in Eq. ( 29), we can study how the perturbation of the matter subsystem leads to modifications of the photonic properties.If we focus on the perturbation in Eq. ( 127), in the velocity representation, the response function describes the time evolution of the photon number in the αmode The residues computed at the polaritonic excitations lim are now connected to dipole and photon number transition moments between the ground and the excited state.Note that to correctly model the coupling to an external field in the velocity representation, we should also include the first interaction term in Eq. ( 118).In the length form, the electric field mode creation (annihilation) operators are described by Eq. ( 37).The response function equivalent to Eq. ( 135) is then given by The response function Eq. ( 135) represents the variation of the number of displacement field modes, in agreement with the transformation in Eq. (33).
In an analogous way, we can consider the time evolution of the field displacement coordinates Eqs. ( 111) and (112).Their time development is obtained by means of the response functions with residues lim that provide information on the transition dipole and photon moments.

Photonic environment perturbation
Classical charge currents J e (r,t) are a source of electromagnetic fields, which means that they can directly perturb the cavity photon field.The interaction term for the Pauli-Fierz Hamiltonian in Eq. ( 17) coupled to a classical external current reads 162 Following Refs. 56,94,187, resolving the external currents into modes which act directly on the α mode of the cavity, in the dipole approximation and the length form, the interaction Hamiltonian can be written as where p α is the photon conjugate momentum of mode α, and j α is connected to the α-mode of the time-derivative of the external current J e . 56,94,187The interaction Hamiltonian Eq. ( 147) generates mode excitations into the system.Due to the coupling of light and matter, this perturbation provides an indirect way to probe the matter and photon-matter correlation properties.
In the length form, the photonic response function describes the time evolution of the displacement field photon number.Interestingly, the more general response function which describes the time evolution of the β -photon number when the system is coupled to the current j α (t) is different from zero even when α = β .This happens because of the coupling to matter degrees of freedom, which provide an indirect interaction between different photon modes.This result opens up the possibility of photon transfer between modes with different frequencies or polarization.In the limit of zero coupling, the response function Eq. ( 151) would differ from zero only when α = β , as photons do not directly interact.Analogously, the response function describes the time evolution of the β -cavity mode momentum when perturbing the α-cavity mode As discussed in the other sections, the transformation in Eq. ( 33) links the response functions for the physical observables in the velocity and length form.For the creation and annihilation operators of the electric photon states, we then need to employ the transformation in Eq. (37).
We can also describe the time evolution of matter degrees of freedom when the system is perturbed through Eq. ( 147).The response function is analogous to Eq. ( 141) as it describes the time evolution of the dipole moment when the photon degrees of freedom are perturbed by j α (t) The response functions in Eqs. ( 135), (140), and ( 154) are examples of response functions which describe the cross-talk between photon and matter degrees of freedom.In fact, in the limit of λ = 0, the eigenstates of the system would be given by the direct product between bare molecular states and photons states.As a consequence, Eqs. ( 135), ( 140), (141), and ( 154) would be identically zero.This means that a perturbation acting on the molecular Hilbert space cannot induce any time evolution in the (decoupled) cavity-radiation space, and vice versa.

Static perturbations and energy derivatives
This section provides several examples of static perturbations on the polaritonic system and their connection to the energy derivatives.a. External electric and magnetic fields.A static external electric field is described by the scalar potential φ e that can be expanded in a Taylor series around the origin where we employed Einstein's summation convention.If we retain only the lowest expansion term, the interaction operator in the Hamiltonian is the familiar electric dipole term Note that the transformation in Eq. ( 33) commutes with this operator, which implies it is the same in both the length and the velocity representation.The static electric dipole polarizability α 0 then reads (158) An analogous interaction term is introduced when the system is immersed in a static uniform external magnetic field, which can be described by the following vector potential 206 In the length form, the interaction Hamiltonian reads where and since we are employing the length form, these quantities refer only to the matter subsystem.The static molecular magnetizability ξ 0 is computed as the second derivative of the energy with respect to the magnetic field The first term in Eq. ( 162) is called paramagnetic contribution, while the second term is called diamagnetic contribution and usually dominates for closed-shell molecules.Magnetizabilities have been extensively studied in the framework of molecular response theory, [219][220][221][222][223][224][225] but the effect of photon-dressing, to the best of our knowledge, is yet to be explored.The evaluation of magnetic properties for approximate wave functions is affected by an origin dependence on the choice of the gauge origin of the external field.In molecular response theory, the origin independence can be recovered by using gauge invariant atomic orbitals (GIAO). 87,226.Spin interactions.Nuclear magnetic resonance (NMR) and electronic paramagnetic resonance (EPR) provide important information on the molecular structure and find wide applications in chemistry.The electron spin is connected to magnetic moment of the electron where g e is the electronic g factor, µ B is the Bohr magneton and s i is the spin operator of electron i.In the same way, the nuclear spin is connected to the nuclear magnetic moment where γ N is the magnetogyric ratio and I N is the spin operator of nucleus N. The spin degrees of freedom modify the Pauli-Fierz Hamiltonian in Eq. ( 30), introducing spin-spin, spin-orbit and spin-external field interaction terms.The cavity environment introduces a twofold modification on NMR and EPR spectra: it modifies the ground state wave function, which is now dressed by the photons, and introduces further terms in the Hamiltonian due to the interaction of the spins with the cavity field.The effects of these additional QED terms are still to be explored and require going beyond the dipole approximation.
We now focus on the NMR properties of closed-shell molecules and disregard the electronic spin contribution to the magnetic field.The total vector potential A tot (r i ) has contributions from the cavity field, the external field, and the nuclear spins and the total magnetic induction is Therefore, the following interaction terms need to be included in the Hamiltonian in Eq. ( 30): where the first term mediates the indirect coupling of nuclear spins, the second term is the nuclear Zeeman interaction, and the last term in Eq. ( 167) is the direct couplings between nuclear dipole magnetic moments.Moreover, using Eq. ( 165) in the velocity Hamiltonian leads to additional interaction terms.The NMR properties are then connected to the following energy derivatives 87 where E (1,1) N is related to the nuclear shielding tensor σ N which describes the electron shielding effects on nucleus N.
At the same time, M,N is connected to the direct D M,N and indirect (electron-mediated) K M,N spin-spin coupling tensors As for magnetizabilities, the computation of NMR properties is affected by origin dependencies, which are usually eliminated using GIAO orbitals.The modifications of the ground state density induced by the QED environment change how the nuclei are shielded from the electrons, affecting the NMR properties.Moreover, the novel interaction terms between spins and the cavity magnetic field can possibly affect spinspin coupling, resulting in modifications to the shape of NMR signals.
c. QED environment perturbations.A modification of the QED environment is reflected in the coupling strength and the frequency of the cavity modes.If we consider a singlemode and perturb the coupling strength the Hamiltonian reads =H l +V (1) +V (2) , where H l is the unperturbed dipole Hamiltonian in the length representation Eq. ( 31), and we have defined the first and second-order perturbators as We therefore obtain expressions for the first and second energy derivatives dV (1) d∆λ ; dV (1) d∆λ ω=0 .
These derivatives measure how sensitive the system is to a change in the coupling strength.In particular, Eq. ( 178) can be interpreted as the static ground state polarizability with respect to the cavity field fluctuations.We note that because of the dipole self-energy in H l , at very large coupling strength, these expressions can be approximated as As we increase λ, from Eq. ( 179) the energy increases since 0|dd|0 is positive definite, but the Hessian in Eq. ( 181) is negative.Therefore, the system is endowed with a stable ground state.Riso et al. showed that for infinitely large coupling strength, the polaritonic ground state energy will converge to the SC-QED-HF energy. 64On the other hand, if the dipole self-energy is disregarded, the remaining terms are While the Hessian is still negative, there is no guarantee from Eq. ( 182) that the energy is increasing, and the ground state energy could decrease indefinitely.The dipole self-energy has been shown to be fundamental for the stability of the ground state. 186 modification of the cavity environment can effectively change the photon frequency From the expansion for ∆ω α ω α , the perturbed Hamiltonian reads The first and second-order perturbation operators are therefore The energy first and second energy derivatives read dV (1) d∆ω α ; dV (1) d∆ω α ω=0 .
These derivatives measure how sensitive the system is to a change in the photon frequency.
In the following section, we develop response theory for approximate wave functions, focusing on QED-HF and QED-CC response theory.

AB INITIO APPROXIMATIONS
In this section, we derive the linear and quadratic response equations for the ab initio QED-HF and QED-CC methods.

QED-HF
The QED-HF is the generalization of the HF model to molecules in a QED environment. 60,61The system is de-scribed through the dipole Hamiltonian where the dipole self-energy term can be included in the molecular Hamiltonian by a proper modification of the one and two-electron integrals where we have defined The QED-HF wave function ansatz is a factorized state where |HF is a single Slater determinant, |0 is the photon vacuum, and c n are expansion coefficients for the photon number states.The QED-HF state is obtained by minimization of the mean value of the Hamiltonian with respect to electronic and photonic parameters.Haugland et al. 60,61 showed that the photonic parameters for the groundstate QED-HF wave function define a coherent state through the following transformation By transforming the dipole Hamiltonian using Eq. ( 201) we get We note that the transformed Hamiltonian Eq. ( 202) is now origin invariant also for charged systems.Minimizing the photon-averaged coherent state transformed dipole Hamiltonian in Eq. ( 202), with respect to the orbital coefficients, we obtain the QED-Fock matrix where F e pq is the standard electronic Fock matrix and the indices a and i refer respectively to virtual and occupied orbitals.The optimization condition that defines the QED-HF molecular orbitals is Brillouin's theorem, F ia = 0, as in standard HF theory.However, the eigenvalues of the Fock matrix F , usu-ally interpreted as orbital energies, are now origin dependent for charged systems.Nevertheless, the total QED-HF energy is origin invariant since the occupied-virtual blocks F ia are origin independent.

Time-dependent QED-HF
Once the optimized QED-HF reference state |R has been obtained, the time development due to an external perturbation V t is parametrized as exp − iΛ(t) |R = exp − iΛ(t) |HF ⊗U QED-HF |0 , (204)   where Λ(t) is the Hermitian operator The κ ai parameters in Eq. ( 205) are orbital rotation parameters, γ α describe the evolution of the QED-HF coherent state, and we compactly defined the orbital rotation operator as Notice that electronic and photonic operators commute, implying that the exponential in Eq. ( 204) can be split into the product of an orbital rotation exponential and a timedependent coherent state for the field Equivalently, we can parameterize the time dependence as where we moved the coherent state transformation to the left.Since Eq. ( 209) and Eq. ( 208) differ by an unimportant phase factor, these parameterizations lead to the same time evolution.By using the Ehrenfest theorem and developing the equations in orders of the perturbation, we obtain the zero, first and second-order equations 1) , R| ∂ ∂t Λ (2) , Ω (0) + ∂ ∂t Λ (1) , Ω (1) + i Λ (1) , ∂ ∂t Λ (1) , Ω (0) |R = R|[V t , Ω (1) ]|R + R|[H, Ω (2) ]|R + i R| Λ (2) , H, Ω (0) |R + i R| Λ (1) , V t , Ω (0) |R + i R| Λ (1) , H, Ω (1) |R − 1 2 R| Λ (1) , Λ (1) , [H, where we also accounted for the possibility that Ω depends on the external perturbation Ω ≡ Ω(V t ).The apexes (0)-( 2) refer to the expansion order and As long as Ω is a one-electron operator or a purely photonic 1/2 operator, the zero-order condition is satisfied due to the QED-HF optimization condition.The first-order equation is non-trivial, and its solution gives us the time evolution of the parameters to first-order in the perturbation.Following the derivation of the response equations in molecular response theory, 84,227 we make use of the operators Making this choice, we ensure the response equations are identical to the ones derived from the time-dependent variational principle. 84,228,229It is convenient to move from the time domain to the frequency domain by performing a Fourier transformation of the parameters such that the operators are the Fourier components of Eqs. ( 213) and ( 214) a. Linear response equations.From Eq. ( 211), performing a Fourier transform and by using Eq. ( 218) and Eq. ( 215), after some algebra, we obtain the standard matrix equation where the vectors X and Y collect the Fourier transformed parameters The right hand side of Eq. ( 222) is the generalized gradient The matrix is the generalized Hessian matrix while is the generalized metric matrix.Therefore, Eq. ( 222) is a generalization of the Casida equations of TDHF in molecular response theory. 230,231Indeed, the explicit expressions of A and B are The electronic blocks have the same definition as the standard Casida matrices, A and B, in the TDHF theory.Nevertheless, we point out that these blocks differ from the bare TDHF matrices for two reasons.First, the QED-HF orbitals and orbital energies differ from the bare molecular ones.Second, the two-electron integrals now contain the dipole self-energy contributions, as shown in Eq. (196).This contribution can be explicitly separated 95 following Eqs.(196) and ( 197) where the ∆ and ∆ matrices now include the dipole selfenergy contribution to the Hamiltonian, while A HF el and B HF el only include the standard electronic integrals.The TD-QED-HF matrices have additional dimensions due to the cavity modes.The coupling between the molecular and the photonic parameters is given by the projection of the transition dipole moments onto the coupling strength vector λ in the off-diagonal terms of Eqs. ( 229) and (230).Moreover, the A matrix has a non-zero photonic block that contains the frequencies ω α of the cavity modes.As expected, in Eq. ( 227) there is no direct interaction between the photon modes, which are coupled indirectly through matter degrees of freedom.We see that these equations are an extension of the familiar TDHF response equations since in the zero coupling limit we recover TDHF solutions and the one photon lines of the empty cavity.We can then write the response function from the time evolution of an operator Ω Inverting equation Eq. ( 222), we identify the linear response function where g Ω has the same structure as g with the operator Ω replacing the perturbation V ω in Eq. (224).We note that the response function fulfils the symmetry relation in Eq. (83).From the generalized eigenvalue equation we obtain the eigenvectors (x i , y i ) T and the spectral decomposition We identify the excitation energies of the system as the eigenvalues of Eq. ( 235).The transition moments are obtained from the residues of the response function, identified as Although the eigenvalue problem in Eq. ( 235) is non Hermitian, it is possible to show that if the computed reference state is close to the ground state, Eq. ( 235) has solutions with real and positive eigenvalues. 84We note that in Eq. ( 222), setting ω = 0, we obtain the static coupled-perturbed QED-HF equations.
Several approximations can be proposed for these response equations.A hierarchy of approximations has been discussed by Yang et al. 95 for their TDDFT-Pauli-Fierz (TDDFT-PF) model, which defined a eigenvalue problem analogous to Eq. (222).They define the Tamm-Dancoff approximation (TDA) by neglecting B el .However, contrary to standard electronic TDA, this is not equivalent to a CI-singles approach.For this reason, we suggest it is more natural to define the TDA approximation by neglecting the whole B matrix, since this is equivalent to a CI problem with singly excited determinants with zero photons and the ground state determinant with one photon.Neglecting B and the dipole self-energy contribution in ∆, they define the TDDFT-Jaynes-Cummings approximation, similar to a JC calculation with all the single excited determinants.
b. Quadratic response equations The quadratic response equations are obtained in a similar fashion as for the linear response. 84,85,232By employing Eq. ( 212) and the set of operators (215), we obtain the second-order response equation where P(ω 1 , ω 2 ) sums the permutations of ω 1 and ω 2 and we used the compact notation To compute the quadratic response function A; B,C ω 1 ,ω 2 , together with the E [2] and S [2] matrices defined in Eqs. ( 225) and ( 226) for the linear response, we need the additional matrices where X = A, B or C. The second-order variation of a time-independent observable A is [Λ (1) , [Λ (1) , where we made the expression explicitly symmetric in the frequencies.From Eq. (241), we can identify the quadratic response function, which has poles where the frequencies or their sum match an excitation energy 233 For the residues and the response function, we then need the following vectors: [2] − (ω 1 + ω 2 )S [2] ] −1 N b (ω 1 ) = [E [2] − ω 1 S [2] ] −1 g B N c (ω 2 ) = [E [2] − ω 2 S [2] ] −1 g C .
As for the linear response equations, these vectors and matri-ces have additional dimensions due to the photonic parameters.

On the definition of the photonic character of excited states
From the QED-Casida equations, it is possible to define a relative electronic/photonic contribution to the polaritonic excitation. 95,99Given a normalized eigenvector of Eq. ( 235), we define the "photonic character" χ n as the sum over the cavity modes α of the squares of the photonic response parameters γ α 95,99 The electronic character of the excitation ρ n is defined as the sum of the orbital rotation κ ai , such that χ n + ρ n = 1 However, such a definition of the excitation character is not straightforwardly connected to photon number nor to electronic excitations in a molecular sense.In fact, for the QED-HF equations, we used the dipole Hamiltonian in the length form.As a consequence of this choice, as pointed out at the end of section (4.2), the photon creation and annihilation operators cannot be simply identified with b † and b.These operators refer instead to the displacement field and therefore include matter contributions from the polarization in Eq. (38).Moreover, as the coupling strength increases, the ground state might gain significant contributions from photon states.
The "photon character" of the excitation becomes, therefore, a slippery concept that requires careful examination also of the ground state.A clarification of this photon character ambiguity could be obtained by implementing the quadratic response equations, from which we obtain the expectation values of observables for excited states.Accordingly, we may compute the following expectation values In the length representation, Eq. ( 245) refers to the mean occupation number of the displacement field number states, while Eq. ( 246) to the mean electric field number states.We can then compare these calculations to the QED-HF ground state expectation values to obtain an estimate of the cavity-field contribution to the excitations.
Employing these less straightforward but more precise definitions might clarify some results presented in the literature.When using the definition in Eq. ( 243), Yang et al. 95 found that the photon contribution of the lower polariton increases with the coupling strength, while the opposite is found for the upper polariton. 95At the same time, the intensity of the lower (upper) polariton increases (decreases) with the coupling strength, which seems in contradiction with the zero oscillator strength associated with the photon states.It is argued that this apparent discrepancy results from the intrusion of higher-energy electronic states, which strengthen the intensity of the lower polariton, overcompensating for the increased photon character of the excitation. 95Nevertheless, in this representation, b † involves both the photons and the matter polarization, which suggests that this interpretation might be revised.
We also note that since QED-HF introduces a coherent state transformation, it would be interesting to investigate the following expectation values where Eq. (249) refers to the occupation of the displacement field coherent states, Eq. (250) refers to the displacement field number states, and Eq.(251) refers to the occupation of the electric field number states.

Equivalent transition moments in TD-QED-HF
Following the same procedure as for standard HF response theory, 84,227 we can show that Eqs. ( 102), (110), and (113), derived for exact wave functions, also hold in the QED-HF response framework.In Tab.(I) we report the TD-QED-HF transition moments for the lower polariton of an ethylene molecule using different basis sets.We report the oscillator strengths in the length f L and velocity f V forms and the transition photon displacement coordinate and momentum We numerically demonstrate the validity of Eqs. ( 102), (110), and (113).Notice that since Eq. ( 102) is fulfilled only in the complete basis set limit, Eqs. ( 252) and (253) will converge only for large basis sets, as seen from the table.At the same time, Eq. ( 113) holds independently of the basis set size and the photon space truncation.Notice that although the equivalence expressed by Eqs. ( 109) and ( 110) do not depend on the quality of the basis set, the computed values of the transition moments change with the basis, as expected from the connection between matter and photon moments expressed by (110).
In particular, the use of diffuse functions strongly affects the computed results.Equivalence relations such as Eqs.( 102) and ( 113) are well known in molecular response theory.How-ever, they have not been explicitly investigated for QED systems.Furthermore, these relations do not hold for the Tamm-Dancoff approximation.The formal equivalence between the transition dipole and velocity momenta can be exploited in the computation of electronic circular dichroism (ECD) and optical rotation.From standard molecular response theory, the ECD spectrum is known to be proportional to the rotational strength The dipole formulation in Eq. ( 256) is origin dependent for a finite basis set, which can lead to unphysical results.If we are interested in the chirality effects promoted by a chiral cavity, such a formulation might be misleading.However, the equivalent velocity form is always origin independent, and Eqs. ( 256) and (257) are equivalent in the limit of a complete basis.Oscillator strength in dipole f L and velocity f V form and transition photon displacement coordinate 0|q α |n and momentum i 0|p α |n for the lower polariton of an ethylene molecule, for different basis sets at TD-QED-HF level.The cavity frequency is set to 0.2695 a.u., and the coupling strength to 0.01 a.u. with field polarization along the transition dipole moment of the first non-dark excitation of ethylene.Eq. ( 102) holds only in the complete basis limit, and convergence between the velocity and length form of the oscillator strength strongly depends on the quality of the basis set.On the other hand, Eqs. ( 109) and ( 110) hold independently of the truncation of the electronic or photonic space.

QED-HF static response equations
For static perturbations, the response equations for the optimized QED-HF state |R can be derived using the same parametrization shown in Eq. (204), which now is timeindependent Here Λ has the same form showed in Eq. ( 205) and may be written compactly as where Θ collects both photon and electron parameters and ξ the respective operators: Following the perturbation expansion of the energy in Sec. ( 5), the first-order energy derivative is which expresses the Hellmann-Feynman theorem for the QED-HF state.In Eq. ( 260), H (1) refers to the first-order interactions.The second-order energy derivative is E QED-HF = R|H (2) |R + R|[Λ (1) , H (1) ]|R , where H (2) is the second-order interaction and Λ (1) is the first-order correction to the parameterization.Since QED-HF state is a variational theory, 60 the optimization conditions are equivalent to those in standard HF 87 R|[ξ, where Eq. ( 262) corresponds to the QED-HF Brillouin's theorem , 60 and Eq. ( 263) to the first-order static response equations.The left-hand side of Eq. ( 263) includes the generalized Hessian matrix, whose blocks are given in Eq. ( 227)-Eq.( 228), and the first-order parameters Θ (1) .The righthand side contains the first-order Hamiltonian, which includes the first-order interaction terms.We note that Eq. ( 263) is equivalent to the time-dependent response equation (222) with ω = 0.

QED-CC
The reference wave function used in QED-CC is QED-HF, where we employ the coherent-state transformed Hamiltonian in Eq. ( 202). 60The QED-CC state is defined as where T is the cluster operator The electronic cluster T e is the standard electronic excitation operator 185 T e = ∑ µ t µ τ µ (266) where |µ is an excited HF determinant.The cluster operator T p includes pure photonic excitations where n is a vector of integers n α referring to the α-mode.
Finally, the interaction operator T int includes simultaneous excitations of matter and cavity modes where fir instance The Schrödinger equation is solved by projection onto the space S spanned by the electronic determinants and the number states 60 In Eq. ( 272), the state |µ, n is a simultaneous excitation to the µ electronic excited determinant and the n photon state The ground state QED-CC equations are therefore µ, n|e −T He T |HF, 0 = 0 (274) HF, 0|H|QED-CC = E QED-CC . (275) The QED-CC dual state is defined as in standard CC theory where tµn are the Lagrangian multipliers. 60,185A collective index t µn can also be defined for the cluster operator such that we can write T more compactly as where As in standard CC theory, QED-CC is based on a hierarchy of approximations where the cluster operators and the projection space are truncated. 60,185Since the space S is larger than in standard CC theory, as it also includes the photonic excita-tions, the QED-CC Jacobian A has additional dimensions In addition to the electronic block A e,e similar to electronic CC, there is also a photonic block A p,p and blocks involving the electronic-photonic parameters From the Jacobian A and the η vector we obtain the Lagrangian multipliers and the equation of motion (EOM)-QED-CC formalism for properties and excited states. 60
We then obtain the equations γk . (307) The zero-order equation in Eq. ( 305) is the ground state multiplier equation, and the vectors η are defined as As for QED-HF, the response matrices and vectors have additional dimensions connected to the electromagnetic degrees of freedom.The X µn and Y νm vectors include, together with the standard CC electronic excitation response parameters, the purely photonic and the simultaneous electronic-photonic response parameters of Eqs. ( 268) and (269).As for the Jacobian matrix A, the F matrix has additional photon and electron-photon blocks Moreover, the perturbation V ω can act on both the electronic and the photon degrees of freedom, as in Eq. ( 147).
we can identify the QED-CC response functions where µn (ω 1 ) µ, n|A|QED-CC + X µn (ω 1 ) Λ|[A, τ µn ]|QED-CC (316) µn (ω 1 , ω 2 ) µ, n|A|QED-CC + X As for QED-HF, the static response equations can be obtained by setting the frequency of the external to zero.The poles and the residues of the response functions can be obtained assuming that A can be diagonalized where ω n is the n-th (real) excitation energy.We further introduce the notation From these equations, we can compute the polaritonic properties. 234Since we employ the QED-HF coherent-state transformed Hamiltonian in Eq. ( 202), the operators and the perturbations in the response functions must be described in the same representation.

QED-CC electronic and photonic excitation character
Since QED-CC is a highly correlated method, defined for the QED-HF coherent-state transformed Hamiltonian in Eq. ( 202), the definition of the electronic or photonic character of the excitation is not trivial.In the EOM framework, Haugland et.al 60 defined the electronic weight w el of the state by means of projection operators where Λ k and R k are the k-th left and right state, P is the projection operator onto the space S of Eq. ( 272) and P el is the projector onto the states of S with zero photons.In the response framework, a similar definition can be obtained by considering the electronic components of the X µn and Y µn vectors, for instance (327) Notice, however, that such a definition does not consider the contribution of the simultaneous electron-photon excitations.
A more consistent definition of the photonic character could be obtained by considering the mean values of the field number operators from the quadratic response function, and comparing them to the ground state values, as discussed for QED-HF

CONCLUDING REMARKS
In this paper, we proposed a systematic discussion of polaritonic response theory based on the well-established molecular response theory routinely employed in quantum chemistry.The fundamental definitions and features of the response functions are still valid, but the explicit treatment of the electromagnetic degrees of freedom allows for novel perspectives.Additional equivalence relations between matter and photonic observables are introduced, and novel ways to probe the system are discussed.Particular care is needed when using different mathematical representations of the operators, as this can lead to misinterpretations of the computed results.We also provided QED-HF and QED-CC response equations that resemble the standard electronic response theory, providing the reader with a general framework for ab inito QED response theory.While significant progress in the theoretical description of polaritonic systems has been made, we emphasize several challenges that future research will have to address.These include a description of light-matter interaction beyond the electric dipole approximation, the issue of disorder in optical devices, the role of collective effects on polaritonic properties and chemical reactions, and chiral polaritonics.

ACKNOWLEDGMENTS
We acknowledge Tor S. Haugland for insightful discussions.
., A.B., and H.K. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (grant agreement No. 101020016).R.R.R and H.K. acknowledge funding from the Research Council of Norway through FRINATEK Project No. 275506.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).We acknowledge computing resources through UNINETT Sigma2-the National Infrastructure for High Performance

TABLE I .
e A e,ep A e,p A ep,e A ep,ep A ep,p A p,e A p,ep A p,p  .