Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics

Experiments at the interface of quantum-optics and chemistry have revealed that strong coupling between light and matter can substantially modify chemical and physical properties of molecules and solids. While the theoretical description of such situations is usually based on non-relativistic quantum electrodynamics, which contains quadratic light-matter coupling terms, it is commonplace to disregard these terms and restrict to purely bilinear couplings. In this work we clarify the physical origin and the substantial impact of the most common quadratic terms, the diamagnetic and self-polarization terms, and highlight why neglecting them can lead to rather unphysical results. Specifically we demonstrate its relevance by showing that neglecting it leads to the loss of gauge invariance, basis-set dependence, disintegration (loss of bound states) of any system in the basis set-limit, unphysical radiation of the ground state and an artificial dependence on the static dipole. Besides providing important guidance for modeling strongly coupled light-matter systems, the presented results do also indicate under which conditions those effects might become accessible.

D riven by substantial experimental progress in the field of cavity-modified chemistry, 1−11 theoretical methods at the border between quantum-chemical ab initio methods and optics have become the focus of many recent investigations. 12−58 The high complexity of a molecular system, which can undergo, e.g., chemical reactions or quantum phase transitions, coupled strongly to photons makes the use of some sort of approximation strategy necessary. A common approach is to use approximation strategies designed for atomic two-level-like systems in highquality optical cavities 59−61 and to apply them to the quite different situation of molecular systems. However, under the generalized conditions of cavity-modified chemistry, contributions in the theoretical description that are usually disregarded, e.g., quadratic terms describing the coupling between light and matter, can become important 15,62,63 and might even dominate the physical properties. [13][14][15][16]64,65 While the existence of these quadratic terms is well-known, 66−73 their origin, interpretation, and consequences are less clear, and when to include them has become the subject of recent intense discussions. 15,40,62,74−82 Here we will elucidate these terms for the most relevant setting of cavity-modified chemistry (i.e., in the Coulomb gauge and in the long-wavelength limit), clarify their origins, physical interpretations, and consequences, and show under which conditions and for which observables they become relevant. This will also highlight the domain of applicability of common approximations that disregard these quadratic terms and at the same time indicate under which conditions substantial influence can be expected, 12 accessible with ab initio techniques such as quantum-electrodynamic density functional theory (QEDFT). 12,16,22,25,83 Before we do so, let us briefly outline the theory that we will consider and collect a set of fundamental conditions that we deem to be important for a reasonable theoretical description.
Any theory that we employ to model coupled light−matter systems should obey certain fundamental constraints. Which ones these are often depends on the specific situation being considered. For instance, in the case of high-energy physics, an adherence to special relativity (physical laws should be Lorentzinvariant) is paramount, and hence, the use of Dirac's equation becomes necessary to capture the behavior of electrons. If we further want to ensure that all of the interactions among the electrons are local and our theory should stay invariant under local phase transformations, we find the Maxwell field coupled to Dirac's momentum operator in a linear (minimal) fashion.
However, the resulting theorywhich, if quantized, is called quantum electrodynamics (QED) and perfectly describes highenergy scattering eventshas many subtle issues. 84 For lowenergy physics, a simplified version, where instead of the relativistically invariant momentum the nonrelativistic momentum is employed, has been shown to be able to resolve many of these issues. 73 The resulting theory of nonrelativistic QED (also sometimes called molecular QED 69,72,85,86 ) is ideally suited to describe atoms, molecules, or solids interacting with the quantized light field. 87−89 However, the coupling between light and matter is only defined up to a phase, and we need to make a specific choice for this phase, i.e., we need to fix a gauge. Changing the gauge or performing a local unitary transformation should not modify physical observables but merely affect their representation in terms of canonical coordinates. While gauge independence is respected by nonrelativistic QED, this constraint is specifically challenging for dimensionally reduced, simplified models. 15,77,80,82 Besides gauge independence, nonrelativistic QED guarantees a set of further intuitive and essential conditions. For instance, it guarantees that the physical observables are independent of the chosen coordinate system (or, in more quantum-chemical terms, where a specific spatial basis is just one of many basis set choices, basis-setindependent), and it also guarantees the stability of matter, i.e., atoms and molecules are stable if coupled to the vacuum of the electromagnetic field. 73 A direct consequence of this fundamental condition is that the combined ground state of light and matter has a zero transverse electric field expectation value. If this were not the case, the system could emit photons and lower its energy. To summarize, a few basic constraints to which we want a theory of light−matter interactions to adhere are the following: all physical observables should be independent of the gauge choice and the choice of coordinate system (e.g., it would be unphysical that the properties of atoms and molecules would depend on the choice of the origin of the laboratory reference frame); the theory should support stable ground states (or else we could not define equilibrium properties and identify specific atoms and molecules); and the coupled light−matter ground state should have a zero transverse electric field (or else the system would radiate and cascade into lowerenergy states).
In From Microscopic to Macroscopic Maxwell's Equations, we will introduce nonrelativistic QED and some of its unitarily equivalent realizations and highlight the physical implications of the associated transformations and further approximations that lead to nonrelativistic QED in the long-wavelength limit. In Necessity and Implications of Quadratic Couplings in the Dipole Approximation, we will illustrate that the aforementioned fundamental physical conditions will not be retained when the quadratic components are disregarded. Finally, in the Summary we will discuss implications and perspectives. We provide further details in the Appendix. In Fundamental Coupling of Light with Matter and the Emergence of Diamagnetism, we present the basic approximations leading to nonrelativistic QED. In The Power−Zienau−Wooley Gauge Transformation and Transverse Basis and Distributions, we provide additional details that complement our discussion regarding the Power−Zienau−Wooley transformation and transverse basis functions. In Spectral Features of Operators, we discuss some implications that go hand in hand with approximating operators. Our discussions will be presented first using a field-theoretical convention in which the four vector potential A μ is given in volts and the four charge-current density j μ is given in coulombs per square meter per second and later in atomic units. By multiplying A μ by 1/c we find the standard convention in terms of volt seconds per meter.

MAXWELL'S EQUATIONS
Classical electrodynamics is at the heart of QED. Consider for instance the inhomogeneous microscopic Maxwell equations: 85 This representation of light−matter coupling is by no means unique, and many different formulations, such as the Riemann− Silberstein 90,91 or macroscopic Maxwell's equations, have been developed over the years. To arrive at the latter, let us first rewrite j = j b + j f , where j b is a bound current and j f a free current, and define the bound charge current as where M is the magnetization and P is the polarization of the matter system. If we then define the displacement field D = ϵ 0 E + P and the magnetization field = − which takes the back-reaction of a given medium on the electromagnetic field into account by the constitutive equations. Clearly, the classical description of electromagnetic interactions can take different forms for which, without further simplifications, none is superior over the other. These forms deviate merely in their choice of canonical variables, the very same variables that will be quantized to reach QED. The electromagnetic field energy is of quadratic form, and therefore, substituting D = ϵ 0 E + P into eq 1 will naturally lead to quadratic self-polarization (P 2 ) and selfmagnetization (M 2 ) terms. While many equivalent ways of formulating Maxwell's equations exist, there will be accordingly also several (unitarily equivalent) forms of the resulting nonrelativistic QED Hamiltonian. Let us in the following see how this equivalence in QED is manifested. A relativistic quantization procedure with a subsequent nonrelativistic limit, as illustrated in Fundamental Coupling of Light with Matter and the Emergence of Diamagnetism, is indeed equivalent to introducing the covariant derivative for the electronic system and then quantizing the resulting gauge field. 12,25,73 This minimal coupling procedure makes the invariance under local phase transformations Ψ′ = e iθ(r) Ψ explicit and, with the Coulomb gauge condition ∇·Â(r) = 0, fixes the local phase θ(r) uniquely. The momentum of each particle is shifted according to −iℏ∇ → −iℏ∇ − (q/c)Â(r), where q is the charge of the particle and the quantized vector potential is where we have defined the transverse polarization vectors ϵ n,λ for mode and polarization (n, λ), 92 for harmonic oscillators with the allowed frequencies ω n = c|k n |. Then, with ĵ0(r) = cn(r) = c∑ i=1 N e +N n q i δ(r − r i ), the nonrelativistic minimally coupled Hamiltonian (including N e electrons and N n nuclei) in the Coulomb gauge reads as Ä Each charged particle then evolves under the influence of the kinetic energy operator [−iℏ∇ − (q/c)Â(r)] 2 and at the same time experiences the instantaneous longitudinal field (the Coulomb potential Ĥi nt,∥ ) created by all of the other charged particles. The nonrelativistic limit of the minimal coupling procedure therefore naturally leads to the appearance of a quadratic term (see also Fundamental Coupling of Light with Matter and the Emergence of Diamagnetism). This quadratic term provides the diamagnetic shift 93 of the bare modes and introduces a lowest allowed frequency, 94 which then removes the infrared divergence. 73 It is therefore not a drawback to have this term. 74 In contrast, it affects, for instance, optical spectroscopy, 95,96 and it is responsible for diamagnetism 97  Let us next for convenience assume linear polarization (i.e., ϵ n,λ = ϵ n,λ * ) and, to connect to the more common formulation, switch to atomic units, such that ϵ 0 = 1/4π, c = 1/α 0 , and the elementary charge (e) is equal to 1. There is a well-known procedure to connect to a Hamiltonian corresponding to the macroscopic where the jth nucleus has the effective positive charge Z j . This implies that all of the physical charges contribute to the bound current, such that j f = 0, ∇·D = ε 0 ∇·E + ∇·P = 0, and therefore D = D ⊥ . However, since the vector potential operator is purely transverse, it also only couples to the transverse part of the polarization operator, which can be expressed in terms of the transverse δ distribution 92 or we use the mode expansion of the vector potential directly. To do so, for notational simplicity we first introduce the abbreviation α ≡ (n, λ), and then we use the fact that the vector-valued functions S n,λ (r) in eq 2 form a basis for the transverse square-integrable vector fields (see Transverse Basis and Distributions). With S α (r) = ϵ α ·S α (r), we find that The resulting Hamiltonian 29,69,72 has the advantage that it can be conveniently expanded in multipoles of the interaction. We note, however, that the validity of such an expansion depends critically on whether it is considered as a perturbation of the wave function or to affect the operator itself and subsequently its self-consistent solution (see Spectral Features of Operators).
The mode expansion provides a consistent regularization such that terms like P̂⊥ 2 (r) are well-defined, a necessity when multiplying δ distributions (see Transverse Basis and Distributions). This avoids the usual auxiliary assumption that some of these terms, which contribute to the polarization self-energy, are only taken into account perturbatively. 79 It furthermore highlights how the condition of transversality of the Maxwell field also affects matter-only operators like the polarization. So far, the only restriction we have employed is that we have considered nonrelativistic particles. However, this simplification is usually not yet enough to allow for practical calculations. In the following we do not want to consider this more general case but rather assume only dipole interactions. This approximation is very accurate provided that the dominating modes of the photon field have wavelengths that are large compared with the extent of the matter subsystem. In the multipole form of the nonrelativistic QED Hamiltonian this means that we discard the integration over s in our transformation and the polarization operator. 69 The Hamiltonian we then find is the same as the one that we get if we approximate Â(r) ≈ Â(r matter ) for the bilinear and quadratic coupling terms. This does not restrict the form of the cavity modes itself but merely their spatial extension in relation to the matter subsystem. In practice, where, e.g., an ensemble of molecules interacts with the cavity mode, this simplification can become questionable. Such an ensemble might extend over macroscopic scales such that individual molecules will experience different couplings. For simplicity, in the following we take r matter = 0 such that S α (0) is real and we can straightforwardly perform the unitary PZW transformation (also called the length gauge transformation) Ĥ= ÛĤÛ † witĥ We accompany this by a canonical transformation that swaps the photon coordinates and momenta (−i∂ q α → −ω α p α , q α → −iω α −1 ∂ p α ) while preserving the commutation relations. 22  with the bare electron mass m e . The nuclear−electron interaction is given by Furthermore, the photonic contribution for M p modes is then given by Ä which incorporates the total dipole moment of electrons and nuclei, i.e., R = −∑ j=1 N e r j + ∑ j=1 N n Z j R j . The resulting bilinear coupling or R itself might be occasionally defined with the opposite sign as a consequence of the inversion symmetry of eq 5. Even if we break the inversion symmetry of the matter Hamiltonian as, e.g., in Coordinate System and Dipole Dependence without Self-Polarization, the photonic symmetry p α ↔ −p α will retain the trivial connection between the two Hamiltonians and their respective observables. In eq 6, M p is a finite but arbitrarily large number of photon modes that are the most relevant modes but in principle run from the fundamental mode of our arbitrarily large but for simplicity finite quantization volume up to a maximum sensible frequency, e.g., an ultraviolet cutoff at the rest-mass energy of the electrons (an extended discussion can be found in Transverse Basis and Distributions). The operator given by eq 6 contains the bilinear matter−photon coupling and the quadratic dipole self-energy term The fundamental light−matter coupling to mode α is then denoted by which depends on the form of the mode functions and the chosen reference point for our matter subsystem. 25,100 This can lead to an increase in the fundamental coupling to a specific mode and is an inherent feature of the physical setup, e.g., the form and nature of the cavity. In the following we will treat λ α and the corresponding frequencies ω α as parameters that we can adopt freely to match different physical situations, motivated by the recent experimental progress to subwavelength effective cavity volumes. 101,102 This also highlights that the self-energy term depending on λ α is influenced directly by the properties of the cavity, i.e, it obtains increasing relevance with decreasing effective mode volume = α S L (0) 1/ 3 and increasing number of participating modes M p .
Importantly, since the PZW gauge is equivalent to Maxwell's equation in matter as introduced earlier, we now work in terms of the purely transverse displacement field 62,72,79 and the transverse polarization operator By construction, the electric field operator in the PZW gauge, which no longer represents the conjugate momentum, becomes The combination of the PZW and canonical momentum transformations changed our canonical operators B̂∝−i∂ p α , D̂⊥ ∝ p α and consequentially the representation of our original creation and annihilation operators to We might yet again define a new harmonic oscillator algebra solely in terms of our new canonical operators to reach a potentially more familiar representation in terms of different aP ZW and aP ZW † . However, we have to consider then that the expression of physical observables in terms of creation and annihilation operators is not invariant under the PZW transformation, i.e., aP ZW ≠ a. Special care has to be taken in how we interpret observables and design possible approximations, as otherwise unphysical consequences arise, as highlighted explicitly in Radiating Eigenstates without Self-Polarization. We also see that in accordance with the Maxwell's equation in matter, by working with D̂⊥ we implicitly take into account the back-reaction (polarization) of matter on the electromagnetic field. The physical field is found with the constitutive relation of eq 8. Finally we note that the PZW transformation has removed the explicit diamagnetic contribution of the current, and the physical current is now equivalent to the paramagnetic current. However, the diamagnetic term has not vanished but is contained in the introduced phase of the coupled light−matter wave function. 62 Let us clearly state a warning at this point. Unitary equivalence or gauge invariance, which implies that we obtain the same predictions in the Coulomb and PZW gauges, is only fulfilled when the full Hilbert space (full basis set) is considered. Any approximation in the molecular or photonic space will violate this equivalence and therefore result in deviating predictions. 72,77,80,82,103 We remain with three different strategies. We can acknowledge this failure and focus on reasonable domains in which the predictions remain in reasonably close agreement, e.g., focus on resonant interactions. 15,72 Alternatively, we can adjust the PZW transformation to compensate for the reduced space accordingly 80 or consider as much of the space as possible, which leads us into the realm of first-principles cavity QED. This breakdown of gauge invariance can have very ACS Photonics pubs.acs.org/journal/apchd5 Article fundamental consequences, as the long-standing debate of the (non)existence of a Dicke super-radiant phase shows. 63,104,105 Disregarding quadratic contributions (Â2/P⊥ 2 ) will consequentially merely allow perturbatively similar results to be obtained. Finally, there is one subtle question left. If we consider many photon modes, they give rise to radiative losses, that is, they constitute the photon bath of the matter subsystem into which the excited states can dissipate their energy. 50 Vacuum fluctuations give rise to effects like spontaneous emission, that is, turning the discrete eigenstates of the closed system described by a Schodinger equation into resonances with finite line width. 50 Selecting furthermore only one or a very limited set of modes α will restrict retardation effects and can lead to unphysical superluminal transfer appearing in the (deep) ultrastrong coupling regime. 106 In this work, however, we are not interested in lifetimes but in equilibrium states of the coupled light−matter system. In this case, instead of keeping many modes we can subsume the vacuum photon bath by renormalizing the bare masses m e and M j of the charged particles, i.e., we use the usual physical masses such as m e = 1 in atomic units and keep only a few important modes that are enhanced with respect to the free-space vacuum.
We have seen that already several approximations have to be employed to arrive at the above Hamiltonian, which represents the usual starting point for most considerations in cavity QED and cavity-modified chemistry. Each approximation restricts its applicability, but the basic physical constraints, i.e., gauge and coordinate system (basis set) independence, existence of a ground state, and radiationless eigenstates, are as of yet conserved. It is now the subject of the following sections to emphasize that ignoring the transverse self-polarization term 1 / 2 ∑ α=1 M p (λ α ·R) 2 and/or the diamagnetic termÂ (0) q c 2 2 2 will inevitably break some of those fundamental constraints. This by no means implies that perturbative treatments that ignore these terms, either by restricting the Hilbert space of the matter subsystem or by perturbation theory on top of free matter observables, might not provide accurate predictions. 15,77,107 However, it shows that care has to be taken when the quadratic terms are disregarded.

COUPLINGS IN THE DIPOLE APPROXIMATION
Let us now consider concretely what happens if we discard the quadratic terms and which further physical constraints we violate. The example will be a simple molecular system, a slightly asymmetric one-dimensional Shin−Metiu model, coupled strongly to a single cavity mode, as illustrated in Figure 1. We subsume the rest of the photon bath in our description approximately into the physical masses of the electron and nuclei. The Shin−Metiu model features one nucleus moving in between two pinned nuclei with a total of one electron and is a paradigmatic model for nonadiabatic electron−nucleus coupling that gives rise to many interesting chemical processes. The Hamiltonian of this model with moving nuclear charge Z = +1 and removed vacuum shift ω/2 is given by where R = −x + ZX is the total dipole moment and the electron− nuclear and nuclear−nuclear potentials are given by Furthermore, we couple the electron and nucleus to a single cavity mode with the frequency ω = 0.00231, which is resonant with the vibrational excitation. We achieve ultrastrong vibrational coupling with ω λ ω = = g/ / 2 0.40748 in atomic units, but by no means do the following results qualitatively depend on the coupling or frequency. The strength of the light−matter interaction solely determines how quickly given effects will be visible, and the selected values are close to those of previous publications in this field of research. It is important to realize that the wavelength associated with this frequency is 1.9724 × 10 5 Å = 19.724 μm and thus differs by about 4 orders of magnitude from the size of the computational box that is considered for the matter system (∼60 Å). Our example is thus safely within the validity of the long-wavelength approximation when considering, e.g., one-dimensional cavities.
No Bound Eigenstates without Self-Polarization. Let us start to investigate the most fundamental problem of discarding ACS Photonics pubs.acs.org/journal/apchd5 Article the quadratic term in the nonperturbative regime: the instability of the coupled system, i.e., the fact that electrons and nuclei fly apart if coupled even in the slightest to the photon field unless we restrict the Hilbert space. 15,62 To illustrate this, we increase the simulation box stepwise by increasing the number of basis functions (grid points), keeping the other parameters fixed, and we present first the light−matter correlated energies as well as the total dipole moment in Figure 2. We find that when the space of allowed wave functions increases (i.e., approaching the basis set limit), the minimum-energy solution without the selfpolarization term does not converge and minimizes the energy (dashed blue line) by increasing the total dipole moment (dashed red line). To put it differently, the system is torn apart, and electrons occupy one side of the simulation box while nuclei occupy the other. On the other hand, with the self-polarization term we see how we quickly approach the basis set limit, such that we have a basisset-independent result (red and blue solid lines). The complete disintegration of the system without the self-polarization term happens at a critical box size that is just marginally larger than a box size leading to converged results when the self-polarization is included. As the box size is increased further, the energies resemble more and more those of an inverted shifted harmonic oscillator, which supports only scattering states. 108 This illustrates that by a small (∼20%) variation of the simulation box we lose the physical character of the model and enter a nonphysical regime. How drastic this effect will be is given by the ratio of the quadratic potential energy 1 / 2 (λ·R) 2 to the energy needed to ionize the system from a given eigenstate, −ε i (assuming noninteracting electrons for simplicity), such that a pure bilinear treatment would be perturbatively reasonable only for −(λ·R) 2 /2ε i ≪ 1. In this sense, the common ratio of the coupling and excitation energies, g/ω, assuming resonance at ε 2 − ε 1 = ω, with a slight adjustment to can be seen as an estimate of how quickly the given eigenstate i will become unstable without the self-polarization component. The extension criterion shown in eq 11 can be motivated by the s i n g l e -p a r t i c l e S c h r öd i n g e r e q u a t i o n Ä , such that the long-range exponential decay of the state i e (e.g., the Bohr radius in the case of the hydrogen atom). The simulation box has to be large enough to at least fit the state i to an amount that we resolve an exponential decay ∼e −1 (which is far from numerical convergence in fact). This provides an estimate of the extension of the eigenstate of interest and its associated self-polarization energy (λ·R) 2 ≈ (λa i ) 2 = −λ 2 /2ε i . While this might provide an orientation for theoretical calculations when instabilities are to be expected without self-polarization, even before the system is torn apart we see that the eigenvalues and the total dipole moment differ noticeably as the basis set increases. Also, other observables change without the self-polarization term, e.g., the nonperturbative Rabi splitting. The observables with selfpolarization remain completely size-independent once a sufficient basis resolution is reached.
Let us illustrate how weakly bound states are affected with the help of a second numerical example. We select a simple onedimensional soft-Coulomb hydrogen atom but screen the nuclear charge Z that binds the electron with to Z = 1 / 20 . We couple this system rather weakly (g/ω = 0.006) with frequency ω = 0.01368 in resonance with its first excitation (when converged) and, as before, increase the size of the simulation box stepwise. Figure 3 illustrates that although the ground state is merely perturbatively affected up to 200 Å, the excited states immediately turn into, for this case, unphysical scattering states. As before, adding the selfpolarization term results in the expected spectrum, very much in contrast to the spectrum without the self-polarization component. The extension criterion λ 2 /4ε i 2 leads to 0.0011 for the ground state and to 0.1664 for the first excited state. While the ratio g/ω = 0.006 gives the impression of rather weak Figure 2. First four light−matter correlated eigenvalues with (blue, solid) and without (blue, dashed) the self-polarization contribution 1 / 2 (λ·R) 2 as well as the total dipole moment R = −x + ZX with (red, solid) and without (red, dashed) self-polarization. While observables start to be converged with a box size around 50 Å, without selfpolarization the system starts to disintegrate already for slightly larger box sizes, as highlighted by the inset. Figure 3. First four light−matter correlated eigenvalues with (blue, solid) and without (red, dashed) the self-polarization contribution for the Rydberg-type weakly bound hydrogen model with grid spacing Δx = 0.8 and 120 photon number states. Until the box reaches a large extent (∼200 Å), the ground state without self-polarization deviatess merely slightly from the correct one. However, the excited states are relatively weakly bound and experience unphysical behavior (and therefore so do the spectra and all of observables involving excited states) even before entering a converged regime. The inset magnifies the unphysical crossover from physically bound into scattering states. The disintegration effect is qualitatively independent of the frequency.

ACS Photonics
pubs.acs.org/journal/apchd5 Article coupling, the extension criterion provides a first indication that the excited states will be substantially affected by the selfpolarization component. While the bilinear interaction reduces the ground-state energy with increasing coupling, the self-polarization contribution counteracts this by an increase in energy and dominates for typical couplings the bilinear contribution, i.e., even the sign, and thus the qualitative behavior, of the energetic shift within the cavity can be altered depending on the presence/absence of the self-polarization. 13,15 This qualitative change is also represented in spatial observables, i.e., the self-polarization term favors a reduced polarizability and thus focuses charge density in domains where charge is already present. 14−16 The bilinear coupling, which furthermore scales with the frequency, is typically weaker affecting the ground state and features the contrary tendency, and their competition determines the qualitative distribution of charges inside the cavity. 15,109 The resulting consequences can, e.g., include a reduced equilibrium bond length 13,14 with an earlier onset of static correlation 14 that could be steered on demand by controlling the polarization of the field and therefore implies interesting opportunities for chemical considerations and electronic devices.
Let us briefly inspect how the same system behaves in the Coulomb gauge instead. Figure 4 illustrates the same correlated eigenvalues with increasing box size as in Figure 3. The Coulomb and PZW gauges lead to accurate agreement, with a numerical difference in energy of less than 10 −7 eV for a box size of more than 40 Å, when both quadratic components are included. Omitting the diamagnetic term leads to a slight negative shift in the correlated eigenvalues and illustrates the relevance of the diamagnetic component, i.e., shifting the photonic excitations (also see Fundamental Coupling of Light with Matter and the Emergence of Diamagnetism). In the Coulomb gauge, couplings between higher excited states rescale lower matrix elements, demanding a well-converged set of electronic eigenstates, and therefore, the bilinear component accounts for the major effect of the self-polarization. 15,77,82 We should recall, however, that the diamagnetic contribution scales with the amount of polarizable material and attains increasing importance as the frequency of the field decreases. In this sense, when the same investigation is performed with a 10 times smaller frequency, ω = 0.00137 = 0.0372 eV, the lowest excited states are photon replicas with an energy spacing of 0.0372 eV between the ground and first excited states with the diamagnetic contribution and 0.0254 eV without, highlighting a substantial deviation.
As a side remark, we note that although the validity of the dipole approximation for high frequencies is questionable, the quadratic self-polarization term guarantees that the highfrequency photons are essentially decoupled from the matter subsystem. If only a purely linear coupling is assumed, then the UV behavior is completely wrong, as photons with arbitrarily high energies still interact with the matter subsystem. 100 Radiating Eigenstates without Self-Polarization. Let us look at another unphysical feature that appears when the selfpolarization is neglected. In the case of the simple Rabi or Dicke model, where the particle is assumed to be perfectly localized (assuming effectively classical particles), the polarization is zero, and we can associate the expectation values of the modes in the PZW gauge with the electric field. However, if we consider an ab initio treatment, this is no longer the case, and we need to use the correct definition of the electric field given in eq 8. In Figure 5 we show the expectation values of the displacement and the electric field as we increase the number of basis functions. This is again the same as allowing the electrons and nuclei to extend over an ever-increasing spatial region, which is equivalent to exploring the full Hilbert space. Also in these observables we see that the system with the self-polarization term included leads to results that are independent of the simulation box size after roughly 45−50 Å (blue and red solid lines). When the box is extended only slightly to ∼60 Å, the system without the self-polarization disintegrates (blue and red dashed lines). By construction, the system with the self-polarization term always obeys the basic constraint of zero electric field, while the system with only bilinear coupling leads for large extensions to an eigenstate with nonzero electric field. This cannot be a physical ground state.
Realizing the connection between the observable field and canonical momentum, let us turn our attention to the number of photons in the ground state. The photon number operator (the . First four light−matter correlated eigenvalues in the Coulomb gauge with (blue, solid) and without (red, dashed) the diamagnetic Â2(r) contribution for the Rydberg-type weakly bound hydrogen model (same parameters as Figure 3). Disregarding the quadratic (diamagnetic) term omits the diamagnetic shift that leads to accurate agreement between the two gauges. The inset magnifies the same region as in Figure 3. Recall that the diamagnetic contribution attains increasing impact as the frequency decreases and the amount of polarizable matter available increases. Figure 5. Displacement fields (blue) and electric fields (red) for the Shin-Metiu model with (solid) and without (dashed) the selfpolarization contribution. While the electric field is independent of the molecular convergence, and therefore even in a restricted subspace the system does not radiate and there is a well-defined equilibrium solution, the displacement field depends on the convergence of the molecular system.

ACS Photonics
pubs.acs.org/journal/apchd5 Article electromagnetic field occupation) is defined in the Coulomb gauge as In the PZW gauge, these annihilation and creation operators are given by eq 9, and as a consequence, the number operator N̂in that gauge is Ä As a result of the change of the conjugate momentum from the electric field to the displacement field, we see that the selfpolarization enters the definition of the photon number operator when we work in the PZW gauge. Without surprise, this leads to different occupations as if we would naively use Ä and we illustrate this difference in Figure 6. The alleged occupation N′ (blue) is higher than the physical occupation N (red) caused by the permanent dipole. Only for two-level models such as the Rabi model do the two definitions agree. 15 Upon comparison of Figures 5 and 6, it is instructive to observe that the behaviors of the displacement field D ⊥ and the naive mode occupation N′ are qualitatively very similar, obtaining relevant nonzero values only after a sufficiently large numerical box size is reached. In contrast, the electric field E remains system-size-independent, and the physical mode occupation N adjusts merely quantitatively to the simulation box. Not surprisingly, ignoring the self-interaction contributions in general leads to different results for different gauge choices.
Coordinate System and Dipole Dependence without Self-Polarization. The Hamiltonian of eq 5 and its variants guarantee that all of the physical observables in equilibrium are independent of the chosen coordinate system. This is obvious if we have a charge -neutral system, where the Hamiltonian of eq 5 is completely translationally invariant. This constraint is physically very reasonable because without a spatial dependence (i.e., the manifestation of the long-wavelength approximation), the electromagnetic field cannot break the translational symmetry of the bare molecular system. If the system is not charge-neutral, e.g., when we consider only electrons in an external binding potential, we no longer have trivial translational invariance. To see this, consider a shift of the origin of the coordinate system along the polarization of the field such that the total dipole moment operator R̂is also shifted. It should be noted that this also changes the origin of the cavity, as the longwavelength approximation enforces that all molecules see the same field (of the now also shifted reference point) of the cavity. However, because of the zero-electric-field condition of a physical ground state, we explicitly know that the relation between the (shifted) dipole moment expectation value ⟨R⟩ of the matter subsystem and the expectation values of the displacement fields is ⟨ ⟩ = ·⟨ ⟩ λ α ω α α p R . 15,25,62 If we then further re-express the light−matter coupling with fluctuation quantities ΔR = R − ⟨R⟩ such that eq 6 becomes Ä we find that at equilibrium the shifts cancel and the only remaining contribution in the Hamiltonian is given by the fluctuations around the mean values. Indeed, the equilibrium wave function in the new coordinate system is just the original wave function translated in space and the photon subsystem coherently shifted. As a consequence, the light−matter-coupled system is invariant under shifts of the origin in equilibrium, and no physical observable has a dependence on the permanent dipole. The fact that the equilibrium properties of light−matter-coupled systems do not depend on a possible permanent dipole is merely a consequence of how particles couple to the transverse electromagnetic field: only currents can interact with photons. A permanent dipole only shifts the photonic displacement field, which is not a physical observable, and the permanent dipoles of molecules contribute only when the combined system is moved out of equilibrium.
Only when the self-polarization term is neglected can an unphysical dependence on the permanent dipole in equilibrium arise. To illustrate that even for small systems and shifts this can have a large influence, we consider the Shin−Metiu model from before, but we slightly charge the complete system by using Z = +1.05e. We then perform a small shift x → x + μ in the coordinate system, solve the corresponding Shin−Metiu model, determine the ground-state electron density n e μ (x), and then translate back to obtain n e μ (x − μ) and compare it with the original (unshifted) solution n e (x). As expected, when the self-polarization is included, we just recover the old density, and n e μ (x − μ) − n e (x) ≡ 0. In contrast, in Figure 7 we show the differences without the self-polarization and find an ever-increasing difference for larger μ (increasing permanent dipole moment). The behavior of the Shin−Metiu model without the self-polarization is clearly unphysical, since observables should not depend on the coordinate system. Any approximate method tailored to perform self-consistent calculations should respect the above coordinate-system independence by retaining the balance between bilinear and quadratic contributions. Consider for example the performance of the nonvariational Krieger−Li−Iafrate approximation for ab initio QEDFT presented in ref 16, which breaks this balance.
It should be noticed that quadratic components also necessarily appear in other situations, i.e., they are indeed a Figure 6. Naive (N′, blue) and physical (N, red) photon occupations in the PZW gauge with (solid) and without (dashed) the self-polarization contribution during the self-consistent solution for the ground state of the Shin-Metiu model.

ACS Photonics
pubs.acs.org/journal/apchd5 Article quite general feature of nonrelativistic Hamiltonians. For example, if nuclear vibrations are approximated as phonon modes, the nonlinear Debye−Waller term, which is proportional to ∇ k ∇ k′ Ĥ, has to be added to the bilinear interaction. 110,111 This term originates from the quadratic elements in a Born− Huang expansion 15 with the very same physical effects as the quadratic components Â2 or P⊥ 2 , e.g., enforcing translational invariance and renormalizing the excitation energies. Collectivity, the Limit of the Dicke Model, and Plasmonic Systems. When considering a system of several molecules, we can separate their instantaneous interactions mediated by longitudinal and transverse polarization fields (in the PZW gauge) by ∫ d 3 r P̂2 = ∫ d 3 r P∥ 2 + P⊥ 2 . Here the first term on the right-hand side corresponds to the Coulomb interaction, and the second term corresponds to the self-polarization contribution. 85 We can further approximately distinguish between situations where the wave functions of the different constituents overlap strongly and situations where there is no strong overlap. The former situation, often termed intramolecular, demands careful consideration of the Coulomb and self-polarization contributions simultaneously, where the Coulomb contribution dominates in most situations. The latter situation of contact-free interactions between matter subsystems, termed intermolecular, leads to a perturbative (dipolar approximation) cancellation of instantaneous interactions such that ∫ d 3 r P̂A·P̂B ≈ 0, and we are left approximately with purely bilinear and retarded interactions between those separated matter subsystems. 112 This situation, however, does not allow longitudinal or transverse interactions to be neglected when calculations are performed locally for one of the subsystems. A consistent calculation considering, for example, molecular rearrangements during chemical reactions due to the influence of cavity-mediated strong light−matter coupling will thus also demand a consistent treatment of Coulomb and self-polarization contributions.
If we now enter into the realm where instantaneous contributions to the intermolecular interaction cancel, it is often instructive to assume that indeed the local (i.e., subsystem) eigenstates are not affected by intermolecular interactions and do not need to be updated during the process. In this case we can perform the pinned-dipole approximation, which implies that each subsystem is localized at a specific position and distinguishable. Starting from eq 5, we can then recover the Dicke model by absorbing the self-polarization contribution perturbatively by renormalizing the mass of the particles (similar to the perturbative treatment of the Lamb shift) such that the effective interaction reduces to the common bilinear coupling. 72 In the case of the pinned-dipole approximation, the bilinear coupling to the displacement field becomes equivalent to a coupling to the electric fields since the local polarization in E⊥ = 4π(D̂⊥ − P̂⊥) is zero by construction. To assume that the quantum subsystems are perfectly localized and distinguishable is in stark contrast to an ab initio quantum-mechanical description of molecules. Thus, applying the Dicke model to deduce the influence of strong coupling on the local molecular states calls for a very careful analysis of all of the applied approximations and their consequences. It furthermore permits physical features such as when charge distributions start to overlap, as is often the case in quantum-chemical calculations, leading to a dependence of the local observables on the surrounding (collective) ensemble. 16 The occurrence of quadratic (i.e., Debye−Waller) terms in the electron−nuclear coupling highlights how general quadratic components are in a nonrelativistic theory. More closely connected to our current situation is the coupling to modes of a plasmonic environment. In principle, if we describe the plasmonic environment as part of the full system, 113,114 the density oscillations of the plasmonic environment are captured in an ab initio description by the Coulomb interaction and the induced transverse photon field, and hence, eq 5 is directly applicable. Let us assume, however, that we are not interested in a self-consistent calculation, can safely disregard contact terms (i.e., ∫ d 3 r P̂A·P̂B = 0), and rather care about perturbative corrections. This consideration will lead to van der Waals (dipole−dipole)-type interactions with different scalings in terms of their intermolecular distance R AB , independent of the choice of Coulomb or PZW gauge. 72 The consideration of large distances subject to significant retardation is described by attractive Casimir−Polder interactions, 115,116  . While those considerations are well-tested and allow for excellent perturbative results, they would not allow a selfconsistent calculation, as these forces would merely result in a collapse of the wave function onto a singular point because of its unbalanced attraction. Assuming for instance a set of harmonic oscillators describing the plasmonic excitations coupled merely bilinearly to the system of interest would introduce divergent forces proportional to −∑ α λ α (λ α ·R). 15 The Coulomb potential, wave function overlap, and the Pauli principle give rise to repulsive components for small R AB , modeled for example by the empirical − R ( ) AB 12 term of the Lennard-Jones potential or the − (e ) R AB term of the Buckingham potential. It is therefore the higher-order components that ensure the stability of matter. A self-consistent treatment of molecules in a plasmonic cavity, which itself is modeled as, e.g., a simple harmonic oscillator, 40,117 thus needs to include higher-order couplings to describe a stable and physical system. Self-consistent calculations would therefore demand extending the quasistatic approximation 118−120 for plasmonic systems such that the plasmonic cluster responds to Figure 7. Differences between the translated (by a shift μ) electron density n e μ (x) and the electron density of the original system n e (x), i.e., n e μ (x − μ) − n e (x), without the self-polarization term. If the selfpolarization is included the difference is always zero, i.e., the equilibrium physics remains independent of the coordinate system and the permanent dipole moment. In the Shin−Metiu model, the charge of the moving nucleus was slightly increased to Z = +1.05, and electronic and nuclear box sizes of 59.27 and 5.93 Å, respectively, were chosen (i.e., before any scattering states appear). All of the other parameter values remained as before.

ACS Photonics
pubs.acs.org/journal/apchd5 Article the coupled molecule. This is precisely the physical origin of the quadratic terms in QED, which allow the photonic or matter system to respond to the coupling by adjusting their excitation energies. For instance, the Â2 part can be subsumed into adjusted mode frequencies and further defines a minimal frequency (i.e., cures the infrared divergence), while the P⊥ 2 term renormalizes the energies of the material, all within the longwavelength approximation. The very same effects should be present for a plasmonic cavity when consistently quantized. Such effects are already observed when ab initio calculations are performed with solely the longitudinal Coulomb interaction. 114 For small clusters, and therefore a small effective volume and high coupling strength, the modification of the response and volume due to the presence of the coupled molecule is nonnegligible and modifies the plasmonic modes of the cluster. A purely bilinear coupling dictates entirely different physics (see Spectral Features of Operators for a detailed discussion), violates all of the aforementioned basic constraints, and leaves such a simplification as inherently perturbative. While state-ofthe-art models might provide insightful perturbative results, the development of corrected models should obtain additional interest, and ab initio calculations could prove beneficial to foster this effort.

■ SUMMARY
It is the very nature of physics that our descriptions are necessarily approximate and that every theory has its limitations and drawbacks. Moreover, even when we have seemingly very accurate theories like QED, we need to reduce their complexity by employing further approximations and assumptions to render them practical. For QED this was historically done by employing perturbation theory and/or restrictions to a minimal set of dynamical variables. Because of their clarity and elegance, the resulting simplified versions of QED are a very good starting point for further investigations, provide for good reasons a common language for a variety of subjects, and have provided tremendous insight over decades of research. Nevertheless, we need to be aware of the conditions under which these simplifications are valid and what their consequences are.
With the recent experimental advances in combining quantumoptical, chemical, and materials science aspects 121 and the subsequent merging of ab initio approaches with quantumoptical methods, it has become important to scrutinize these common assumptions. 12 In this work, we have elucidated and illustrated the consequences of discarding quadratic terms that arise naturally in nonrelativistic QED. Omitting them breaks gauge invariance, introduces a dependence on the coordinate system (or basis set), leads to radiating ground states, introduces an artificial dependence on the total dipole moment, and in the basis set limit leads to disintegration of the complete system. However, many of these effects can be mitigated if one works perturbatively or restricts the parameter space. This is in accordance with many years of successful application of such approximations but also highlights their limits of applicability. However, estimates of their applicability, such as the extension criterion (eq 11) discussed in No Bound Eigenstates without Self-Polarization, nowadays become relevant for practical calculations. Certainly when strong coupling between light and matter modifies the local matter subsystem, as is suggested by recent experimental results, 3,11,99,122−124 the quadratic terms can become important and determine the physical properties.
When looking beyond the simple Rabi splitting of spectral lines, which is the accepted indicator of the onset of strong coupling, other observables that contain further information about the matter subsystem should be able to highlight the necessity to modify the common Dicke or Rabi models, e.g., as demonstrated in ref 76. By consideration of photonic and matter observables at the same time, the dipole-approximated bilinear coupling can be further scrutinized, the influence of quadratic coupling terms can be revealed, and effects that are due to spatially inhomogeneous fields (beyond the long-wavelength approximation) can be observed. Furthermore, when the light− matter coupling causes bilinear, self-polarization, and Coulomb interactions to act on comparable energy scales, nonperturbative effects can be expected. At this point it is important to realize that this statement also holds spatially, i.e., that while a coupling might be considered small for certain bond lengths/extensions of the molecular system, at other parts or on other scales it might become substantial. The extension criterion (eq 11), which weighs the divergent forces (which increase with increasing spatial extension) against the ionization energy, is motivated precisely with this spirit in mind, providing the g/ω complementing parameter estimate. Consider for example the binding curve of a molecule. Probing the dissociative regime at large distances will change the ratio of the aforementioned contributions until van der Waals-type interactions containing retardation effects have to be considered, a problem that is also of chemical interest. 125 Not just the equilibrium distance of molecules will change, but especially their behavior in the stretched configuration will be affected, 13,14 a feature essential to describe chemical reactions. For relatively large systems, which are yet still small compared with the relevant wavelengths of the photon field, stronger effects would be expected. In the simple models presented here, we could have used a smaller coupling strength yet a spatially more extended system, and we would have found similar effects. Dynamics that probes the long-range part of potential energy surfaces should be affected more strongly, and this is especially true compared with dynamics due to classical external laser fields in the dipole approximation ignoring the self-polarization term. While our focus remained on the single-molecule limit, an additional essential scale of the system is represented by the number of charge carriers, amplifying the dipole moment, polarizability, and therefore the self-polarization contribution. Extended systems (e.g., solids and liquids) and molecular ensembles with charge contact (e.g., biomolecules) are therefore expected to experience quite sizable influences by quadratic components, i.e., perturbatively seen renormalizing photonic or matter excitations due to the collective light−matter interaction. 15 Let us here briefly show how QED and its nonrelativistic limit can be set up starting from classical electrodynamics. In vector potential form, the microscopic description of the electromagnetic fields is given by ACS Photonics pubs.acs.org/journal/apchd5 Article = − ∂ electronic structure. The resulting quantum-optical models are designed to predict specific photonic observables and are consequentially limited in their predictability for the matter subsystem. 14,15 Recent interest in strong light−matter interaction is calling now for a consistent treatment of those limits that have historically been perceived as complementary. Under certain assumptions, the diamagnetic term can indeed be absorbed by a redefinition of the frequencies and polarizations of the field modes. 15,75,93,94 These redefinitions depend on the matter subsystem (more specifically the number of charged particles) and lead to a diamagnetic shift of the photon field, which can be observed experimentally. 98,127,128 Since the difference between the bare and diamagnetically dressed photonic quantities goes as N V / , where N is the total number of particles and V is the quantization volume, it is often argued 93 that one can use the bare quantities for finite systems. This is not entirely correct. The same argument would predict that the coupling between light and matter (see eq 7) would be zero. The reason for nonzero coupling lies in the fact that as the quantization volume becomes larger (approaching free space), the number of modes in any frequency interval approaches infinity as well. Thus, while the coupling to an individual mode indeed becomes zero, the coupling to the continuum of modes is nonzero. Therefore, when we keep individual modes in our theoretical description, we effectively treat a small but finite frequency interval of modes. This frequency interval can be related to the effective mode volume, since it gives the spacing between the effective modes. Consequently, we also need to dress the photon operators diamagnetically.

The Power−Zienau−Wooley Gauge Transformation
While here we performed the unitary PZW transformation after choosing the Coulomb gauge quantization, which leaves the vector potential operator invariant but leads to an adjusted conjugate photon field momentum and coupling, one can equivalently use the PZW (multipolar) gauge of the field to perform the quantization procedure. 68 This gauge is connected to the Coulomb gauge by an adjustment of the phase of each particle by θ(r) = −qα 0 ∫ 0 1 r·A(sr) ds. 68 This extra phase removes the explicit diamagnetic component from the physical current but also assigns a longitudinal component to the vector potential that can be associated with the Coulomb interaction. While the PZW gauge features purely transverse light−matter coupling, similar to the Coulomb gauge, it mixes light and matter degrees of freedom in accordance with the macroscopic Maxwell's equations. [66][67][68]79,85

Transverse Basis and Distributions
For an arbitrary cavity geometry, the vector-valued eigenfunctions can be very complicated, deviating from simple plane waves. Still, this basis S α (r) is assumed to obey the condition ∇·S α (r) = 0. In this case, we need to perform the mode expansion of the vector potential operator and the polarization operator with the corresponding modes. Selecting a basis defined the representation of the δ distribution and the according polarization. It is important to note that while there are many equivalent representations of the δ distribution, e.g., by using different basis sets S α (r), multiplications of distributions are not uniquely defined. 129,130 Indeed, the origin of the divergence in quantum field theories stems from the fact that (operator-valued) distributions are multiplied 131 and a regularization and renormalization procedure needs to be employed to give a finite answer. The usual method of regularization is equivalent to introducing a cutoff in the mode expansion α, and hence, by keeping this explicit instead of working with an unspecified representation of the δ distribution we avoid nonuniqueness problems. 79 We could straightaway also use the full infinite space, 73 but this would just make the notation unnecessarily complicated, as the above Hamiltonian converges in the norm-resolvent sense to the infinite-space Hamiltonian for M p → ∞. 132 Hence, the above Hamiltonian can be made equivalent to the full infinite-space Hamiltonian if we increase the quantization volume.

Spectral Features of Operators
Most arguments for performing a multipole expansion for a Hamiltonian are based on perturbation theory, where local properties derived from a fixed wave function only slightly depend on higher-order contributions of a perturbing operator. 72 However, this does not mean that such arguments still apply for nonperturbative considerations, i.e., when the operator itself is changed and we solve the resulting equation self-consistently. Indeed, if we consider the influence of such expansions on the Hamiltonian directly, the opposite is usually true: the highest-order terms determine all of the basic properties. An instructive example is a one-dimensional model atom with Ĥ= −∂ x 2 /2 + v(x), where v(x) is some binding potential centered at x = 0 with v(x) → 0 for x → ∞. Its spectrum as a self-adjoint operator in L ( ) 2  contains both bound states (eigenfunctions exponentially localized around x = 0) and scattering states (distributional eigenfunctions corresponding to the continuous spectrum). In such cases, a harmonic approximation for certain ground-state properties, v(x) ≈ v(0) + v′(0)x + v″(0)x 2 /2, is reasonable (assuming v″(0) > 0), and the perturbative influence of higher-order terms proportional to x n with n > 2 is minor. However, if we consider the actual Hamiltonian and treat higher-order terms proportional to x n nonperturbatively in L ( ) 2  , we see that either we have an operator that is unbounded from below (having no ground state) with a purely continuous spectrum (no eigenfunctions but only scattering states) for odd n 62 or bounded from below with only bound states for even n (again assuming that all v (n) (0) > 0). Thus, all of the basic properties are determined only by the highest order of n. We therefore find that an expansion of an operator becomes meaningful only if we also indicate whether we consider it perturbatively for a fixed wave function or nonperturbatively. In this work, we focus on the nonperturbative situation. Let us also mention that as an alternative to perturbation theory, a nonperturbative consideration but on a different Hilbert space of a restricted domain ϑ ⊂ 3  becomes possible, i.e., we consider the operators on L 2 (ϑ) with appropriate boundary conditions or on a restricted state space {Ψ Ψ} ⊂ L , ..., ( ) r 1 2  . As illustrated in Necessity and Implications of Quadratic Couplings in the Dipole Approximation, this will render the restricted domain or subset a relevant parameter for the theoretical prediction since the physical properties can then crucially depend on this parameter.