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 the chemical and physical properties of molecules and solids. While the theoretical description of such situations is usually based on nonrelativistic quantum electrodynamics, which contains quadratic light-matter coupling terms, it is commonplace to disregard these terms and restrict the treatment 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 their relevance by showing that neglecting these terms 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 of strongly coupled light-matter systems, the presented results also indicate conditions under which those effects might become accessible.

Driven 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 high-quality 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−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 Lorentz-invariant) 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 high-energy scattering events—has many subtle issues.[84] For low-energy 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-set-independent), 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 lower-energy 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.

From Microscopic to Macroscopic 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 = jb + jf, where jb is a bound current and jf a free current, and define the bound charge current aswhere M is the magnetization and P is the polarization of the matter system. If we then define the displacement field D = ϵ0E + P and the magnetization field , we end up withwhich 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 energyis of quadratic form, and therefore, substituting D = ϵ0E + P into eq will naturally lead to quadratic self-polarization (P2) and self-magnetization (M2) 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 Ψ′ = eiθ(Ψ 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 iswhere we have defined the transverse polarization vectors ϵ for mode and polarization (n, λ),[92] and the creation and annihilation operators can be expressed in terms of displacement coordinates and their conjugate momenta for harmonic oscillators with the allowed frequencies ω = c|k|.

Then, with ĵ0(r) = cn̂(r) = c∑qδ(r – r), the nonrelativistic minimally coupled Hamiltonian (including Ne electrons and Nn nuclei) in the Coulomb gauge reads asEach 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 Ĥint,∥) 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] and hence implies very important physical processes, such as the famous Meissner effect. Recent theoretical[14,94,98] and experimental[76] studies focused on ultrastrongly coupled light–matter dynamics as well as the prediction of enhanced electron–phonon coupling[64,65,99] highlight the non-negligible influence of the collective diamagnetic shift.

Let us next for convenience assume linear polarization (i.e., ϵ = ϵ*) 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 Maxwell’s equation by employing the unitary Power–Zienau–Woolley (PZW) transformation (see also The Power–Zienau–Wooley Gauge Transformation):for all charged particles contributing to the polarization,where the jth nucleus has the effective positive charge Z. This implies that all of the physical charges contribute to the bound current, such that jf = 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(r) in eq form a basis for the transverse square-integrable vector fields (see Transverse Basis and Distributions). With Sα(r) = ϵα·Sα(r), we find thatThe 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) ≈ Â(rmatter) 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 rmatter = 0 such that Sα(0) is real and we can straightforwardly perform the unitary PZW transformation (also called the length gauge transformation) Ĥ = ÛĤÛ† with . We accompany this by a canonical transformation that swaps the photon coordinates and momenta (−i∂ → −ωαpα, qα → −iωα–1∂) while preserving the commutation relations.[22] The nonrelativistic QED Hamiltonian then reads asIn eq , the nuclear Hamiltonian iswith the bare nuclear masses M. The electronic Hamiltonian iswith the bare electron mass me. The nuclear–electron interaction is given byFurthermore, the photonic contribution for Mp modes is then given bywhich incorporates the total dipole moment of electrons and nuclei, i.e., R = −∑r + ∑ZR. The resulting bilinear coupling or R itself might be occasionally defined with the opposite sign as a consequence of the inversion symmetry of eq . 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 , Mp 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 contains the bilinear matter–photon coupling and the quadratic dipole self-energy term .

The fundamental light–matter coupling to mode α is then denoted bywhich 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 and increasing number of participating modes Mp.

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 operatorBy construction, the electric field operator in the PZW gauge, which no longer represents the conjugate momentum, becomesThe combination of the PZW and canonical momentum transformations changed our canonical operators B̂ ∝–i∂, D̂⊥ ∝ pα and consequentially the representation of our original creation and annihilation operators toWe might yet again define a new harmonic oscillator algebra solely in terms of our new canonical operatorsto reach a potentially more familiar representation in terms of different âPZW and âPZW†. 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., âPZW ≠ â. 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 . 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 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 Schödinger 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 me and M of the charged particles, i.e., we use the usual physical masses such as me = 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(λα·R)2 and/or the diamagnetic term 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.

Necessity and Implications of Quadratic 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 . 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 bywhere R = −x + ZX is the total dipole moment and the electron–nuclear and nuclear–nuclear potentials are given bywhere erf(z) is the error function. For the following calculations, we consider the parameter values Z+ = 1, Z– = 1.05, M = 1836me, L = 18.8973, Rc = 2.8346, and Rf = 3.7795 with electronic and nuclear spacings of Δ = 0.4 and Δ = 0.04 between the equidistant grid points and 40 photon number states. 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 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 × 105 Å = 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.

Figure 1

Illustration of the Shin–Metiu model inside the cavity. The two outer nuclei with the distance L and charges Z± are fixed in position, while the central nucleus and electron can move freely within one dimension.

No Bound Eigenstates without Self-Polarization

Let us start to investigate the most fundamental problem of discarding 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 . We find that when the space of allowed wave functions increases (i.e., approaching the basis set limit), the minimum-energy solution without the self-polarization 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.

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 self-polarization the system starts to disintegrate already for slightly larger box sizes, as highlighted by the inset.

On the other hand, with the self-polarization term we see how we quickly approach the basis set limit, such that we have a basis-set-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, −ε (assuming noninteracting electrons for simplicity), such that a pure bilinear treatment would be perturbatively reasonable only for −(λ·R)2/2ε ≪ 1. In this sense, the common ratio of the coupling and excitation energies, g/ω, assuming resonance at ε2 – ε1 = ω, with a slight adjustment tocan 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 can be motivated by the single-particle Schrödinger equation in the limit r → ∞, v(r) → 0, ∇2 → ∂2, such that the long-range exponential decay of the state is defined by its characteristic extension (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)2 = −λ2/2ε. 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 self-polarization 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 one-dimensional 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 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 self-polarization term results in the expected spectrum, very much in contrast to the spectrum without the self-polarization component. The extension criterion λ2/4ε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 coupling, the extension criterion provides a first indication that the excited states will be substantially affected by the self-polarization component.

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.

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 illustrates the same correlated eigenvalues with increasing box size as in Figure . 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.

Figure 4

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 ). 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 . Recall that the diamagnetic contribution attains increasing impact as the frequency decreases and the amount of polarizable matter available increases.

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 high-frequency 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 self-polarization 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 . In Figure 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.

Figure 5

Displacement fields (blue) and electric fields (red) for the Shin-Metiu model with (solid) and without (dashed) the self-polarization 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.

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 electromagnetic field occupation) is defined in the Coulomb gauge asIn the PZW gauge, these annihilation and creation operators are given by eq , and as a consequence, the number operator N̂ in that gauge isAs a result of the change of the conjugate momentum from the electric field to the displacement field, we see that the self-polarization 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 useand we illustrate this difference in Figure . 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 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.

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.

Coordinate System and Dipole Dependence without Self-Polarization

The Hamiltonian of eq 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 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 long-wavelength 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 .[15,25,62] If we then further re-express the light–matter coupling with fluctuation quantities ΔR = R – ⟨R⟩ such that eq becomeswe 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 neμ(x), and then translate back to obtain neμ(x – μ) and compare it with the original (unshifted) solution ne(x).

As expected, when the self-polarization is included, we just recover the old density, and neμ(x – μ) – ne(x) ≡ 0. In contrast, in Figure 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.

Figure 7

Differences between the translated (by a shift μ) electron density neμ(x) and the electron density of the original system ne(x), i.e., neμ(x – μ) – ne(x), without the self-polarization term. If the self-polarization 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.

It should be noticed that quadratic components also necessarily appear in other situations, i.e., they are indeed a 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 ∇∇Ĥ, 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 ∫d3r P̂2 = ∫d3r 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 ∫d3r 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 , 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 Ê⊥ = 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 is directly applicable. Let us assume, however, that we are not interested in a self-consistent calculation, can safely disregard contact terms (i.e., ∫d3r 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 RAB, 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] which scale as . For smaller RAB, retardation might be omitted, and we would enter the realm of instantaneous attractive interactions captured by the London dispersion potential, which scale as . While those considerations are well-tested and allow for excellent perturbative results, they would not allow a self-consistent 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 RAB, modeled for example by the empirical term of the Lennard-Jones potential or the 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 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 long-wavelength 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 non-negligible 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-of-the-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 quantum-optical, chemical, and materials science aspects[121] and the subsequent merging of ab initio approaches with quantum-optical 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 ) 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 ), 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,62,75] Recent investigations of cavity-enhanced electron–phonon coupling[64,65] and its role for superconductivity[99] might indicate the substantial scientific impact of this realization. Exploring these situations where our theoretical descriptions begin to differ strongly thus holds promise to reveal further yet undiscovered effects.[12]