Avoiding gauge ambiguities in cavity quantum electrodynamics

Systems of interacting charges and fields are ubiquitous in physics. Recently, it has been shown that Hamiltonians derived using different gauges can yield different physical results when matter degrees of freedom are truncated to a few low-lying energy eigenstates. This effect is particularly prominent in the ultra-strong coupling regime. Such ambiguities arise because transformations reshuffle the partition between light and matter degrees of freedom and so level truncation is a gauge dependent approximation. To avoid this gauge ambiguity, we redefine the electromagnetic fields in terms of potentials for which the resulting canonical momenta and Hamiltonian are explicitly unchanged by the gauge choice of this theory. Instead the light/matter partition is assigned by the intuitive choice of separating an electric field between displacement and polarisation contributions. This approach is an attractive choice in typical cavity quantum electrodynamics situations.

Despite superficial differences the theory is still gauge invariant and so, irrespective of the choice of χ-field, H m will describe the same physics. From Eqn. (5b) it is clear that the transverse part of the vector potential is gauge invariant. This is quantised in the usual way by satisfying the commutation relation where f k = (2ν k V ) −1/2 is the coupling strength to the mode with frequency ν k = c|k| in volume V , kλ are polarisation unit vectors orthonormal to k, withâ kλ (â † kλ ) being a gauge-dependent photon annihilation (creation) operator [1][2][3][4]. In the following we will henceforth derive equations for the simple dipole system in the main text instead of the generic charge distribution. From Eqn. (SM 1) the magnetic field follows as The canonical momentum of the field is gauge dependent, but has the same expression in each gauge found through the commutator , and is given by

A. Coulomb gauge
The Coulomb gauge Hamiltonian is defined by choosing ∇χ = A and results in φ C = −ε 0 E and so Π C = −ε 0 E ⊥ where we have replaced superscript 'prime' with subscript 'C' to denote the gauge [5]. The photon field is purely transverse with A C = A ⊥ with conjugate field Π C = −ε 0 E ⊥ . This also means that the Coulomb gauge canonical momentum is On substituting these into the arbitrary-gauge Hamiltonian Eqn. (SM 1) for the simple system described in the main text, one finds that where U ext is an externally applied potential and U b = ε0 2 d 3 x [E (x)] 2 is the electron-hole electrostatic potential which contains the divergent self-energies and interaction between the particles. This potential follows from Gauss's law Eqn. (1c) and the lack of free charge, which leads to where V b = V e +V h is the sum of electron and hole potentials. Integration by parts, assuming fields vanish at infinities, and again using Gauss's law leads to where U e−e and U h−h are the divergent self energies of the electron and hole and U e−h is the electrostatic energy due to the mutual attraction. The A ⊥ , B and Π ⊥ C fields are quantised as in Eqns. (SM 1), (SM 2) and (SM 3) albeit with the arbitrary gauge ladder operators replaced with Coulomb-gauge ones. Using these, the terms inside the integral in Eqn. (SM 5) becomes the usual bosonic field energy term and so the Coulomb-gauge Hamiltonian becomes In the Coulomb gauge, the light/matter coupling is between the dipole canonical momentum and the gauge invariant transverse vector potential, of the form p C · A ⊥ . From Eqn. (7b), we see that the field canonical momentum is proportional to the transverse electric field, Π C = −ε 0 E ⊥ . We now expand in the basis of N + 1 energy eigenstates. The system Hamiltonian is formed from Finally, the Coulomb-gauge light matter interaction is written in this basis as follows where nn ≡ n − n , d nn = −e n,C |r | n ,C andσ y nn = −i (| n ,C n,C | − | n,C n ,C |). In the second line of Eqn. (SM 9) we have written p = mṙ, which is only true in Coulomb gauge in the absence of a vector potential, and is why the unperturbed Hamiltonian is used in the expectation valueṙ = i[H 0 , r]. Finally, the Coulomb-gauge Hamiltonian written in the energy eigenbasis is where the ladder operators are implicitly in Coulomb gauge. We have made the electric dipole approximation which allows us to evaluate the vector potential at the centre of the dipole at the origin at the coordinate system, i.e. A ⊥ ≡ A ⊥ (0). Additionally, we have ignored the vacuum energy of the photon field, and the divergent self energies of the electron and hole, U e−e and U h−h , and also the electron-hole interaction U e−h which is itself divergent in the electric dipole approximation.

B. Multipolar gauge
The derivation of the multipolar gauge Hamiltonian is more involved and the reader is pointed to [5,6] for the details. It can be summarised by the choice ∇χ = A ⊥ − Φ 0 /e. This leads to φ m = P and therefore Π m = −D = −D ⊥ where the lack of longitudinal component follows from Gauss' law because there are no free charges. Subscript 'm' labels quantities in the multipolar gauge. By substitution into Eqn. (9), one finds that the Hamiltonian in the multipolar gauge is The EM fields are quantised in the same way as in Coulomb-gauge, despite Π ⊥ m and Π ⊥ C corresponding to different physical fields. In the multipolar gauge, the electrostatic interaction is governed by the longitudinal component of the polarisation 2 , which is clear from Eqn. (12a) and leads to the same expression given in Eqn. (SM 6). In the multipolar gauge, the photon field is much more complicated, Π m = −D ⊥ = −E ⊥ − P ⊥ , and contains both dipole and transverse field degrees of freedom. The light/matter coupling is between the mechanical dipole position and the canonical momentum of the field, d 3 x P ⊥ (x) · D ⊥ (x). Now expanding the matter into the energy basis, the interaction becomes where we used the dipole operator d = −er and in the last line we have made the electric dipole approximation. Within the EDA, the dipolar self-energy term can be rewritten using the transverse delta function as [4] (SM 13) Applying these, expanding the kinetic energy terms of the eigenbasis as in Eqn. (SM 8) and ignoring the magnetic contribution in Φ, we find the standard multipolar Hamiltonian where the ladder operators are implicitly in multipolar gauge. The similarities between the multipolar Hamiltonian derived here in the A-field approach and the gauge-invariant Hamiltonian we have derived in the new C-field approach in Eqn. (26) should be noted (discussed further in the main text). Furthermore, the expansion of the dipole operator into the N + 1 energy eigenstates is where we have assumed that the binding potential gives rise to eigenstates for which the diagonal elements of the dipole transition matrix are zero and d nn = −e n,m | r | n ,m . Note that the coupling strength between levels n and n scales like d nn ν in the multipolar gauge and d nn nn in Coulomb gauge. The increasing coupling strength with energy level separation in Coulomb gauge is related to the breakdown of gauge invariance for finite dipole level truncation, as reported by Refs. [7,8].

SUPPLEMENTARY MATERIAL 2: ADDING AN IONIC LATTICE IN THE NEW APPROACH
An ionic lattice can be added to the new approach consistently under the assumption that it does not contribute any macroscopic current. This is justified given that ions vibrate about a mean position. To do so, we note that the partitioning of ρ and J into free and bound charges was arbitrary. We could just as well make the partition ρ = ρ d + ρ i and J = J d + J i where subscripts 'd' and 'i' denote charges that are bound in dipoles, and charges that are ions, respectively. The lattice is described by ions of charges Q k at positions R k with charge density and J i = 0. Note that to conserve charge without a current density, the charge density must not be an explicit function of time.
A field P, analogous to the polarisation field P, is then defined to be sourced by the dipoles: ∇ · P = −ρ d and a field D analogous to the displacement field D to be sourced by the ions: ∇ · D = ρ i . If there are no ions, then the new C-field theory outlined in the main text is unchanged so long as P → P and D → D because ∇ · E = ρ is still satisfied for E = (D − P)/ε 0 . Including the ions, the Maxwell equations Eqns. (13) become where we have set J i = 0. Analogously to the main text, we write the fields in terms of vector potentials to satisfy Eqns. (SM 17). As before H = ∇C 0 +Ċ, however, the displacement field analogue requires a longitudinal component with the restrictions that ∇ · D = ρ i andḊ = 0. The latter restriction comes from Eqn. (SM 17a) and the lack of a macroscopic free current. Note that as before D ⊥ (x) = ∇ × C(x). These restrictions are satisfied by choosing where V i is the electrostatic potential of the ionic lattice. The inclusion of an ionic lattice to the Lagrangian Eqn. (??) brings an additional kinetic term along with the replacements P → P and In order for the lattice to have been added consistently we must be able to derive the Lorentz force acting on the ions by using the Euler-Lagrange equation with ion position R k . The velocity derivative clearly results in ∂L/∂Ṙ k = M kṘk , though the position derivative is more difficult. To start we rewrite Eqn. (SM 19) as We can then write the ion position derivative of the Lagrangian as where we have enforced that P and D ⊥ are not functions of R k and that ∂/∂R k ×D = 0 because D is a longitudinal field with respect to R k . The j-th component of the integrand is therefore where the sum over i is implied. The longitudinal and transverse delta functions are [9] where δ ij is a Kronecker-delta and δ ij (x) ≡ δ ij δ(x). Therefore, which gives the Lorentz force The lack of magnetic force and transverse electric force originates from assuming that the ions do not generate any current. Therefore, the ions are not affected by magnetic forces and the electrostatic force is conservative which means it has zero curl and so no transverse component. SinceḊ = 0, the mechanical degree of freedom D does not have a canonical momentum. The ion position does, however, have a canonical momentum given by which brings an additional kinetic energy term to the total HamiltonianH. This Hamiltonian is then found to be (using As described in the main text, we can quantise the C ⊥ -field as in Eqn. (25). The radiation fields are quantised exactly as without the ionic lattice. That is, for the C ⊥ -field quantised as in Eqn. (25), D ⊥ takes the form and the magnetic field is just the canonical momentum. Therefore, using the canonical commutation relation where the second equality follows from ∇ · B = 0, we find that After quantising the radiation fields in this way, it is clear that the first line of Eqn. (SM 29) will lead to the Hamiltonian describing the dipole interacting with the light field that is given in the main text in Eqn. (26), so long as the magnetic effects in Φ are ignored and the electric dipole approximation is made to simplify P(x) ≈ −erδ(x).
Quantising the second line of Eqn. (SM 29) leads to a description of the influence (on the dipole) of vibrations in the ionic lattice from their mean positions. From Eqn. (SM 19) we can see that the total energy of the phonon field is described by the terms where b qµ (b † qµ ) is the annihilation (creation) operator of a phonon of wavenumber q and polarisation µ. Details of the derivation of Eqn. (SM 32) can be found in, for example, Refs. [1][2][3]10]. The final term in the Hamiltonian, d 3 x D · P, describes the interaction between the phonons and the dipole. After linearising this interaction by assuming that the ions are only displaced from equilibrium by small amounts, one can then show that in the EDA [1-3, 10] which is the usual phonon interaction encountered throughout the literature where the coupling strength of dipole level n to phonon mode with wavenumber q and polarisation µ is g nqµ . We can finally write the full Hamiltonian describing the dynamics of a dipole in an EM field interacting with a vibrating ionic lattice (with V = V 0 = 0): In this section, we will derive the equations of motion and the Hamiltonian of the C-field theory from the Lagrangian given in the main text by Eqn. (18). For ease of reading, this Lagrangian is where D = ∇ × C and H = ∇C 0 +Ċ, which automatically satisfy Gauss's law ∇ · D = 0 and Maxwell-Ampere's law ∇ × H =Ḋ. Note that the equations of motion for the Lagrangian transformed by Eqns. (17) must be identical because the transformation preserves the equations of motion by construction.

A. No magnetic monopoles
The non-existence of magnetic monopoles means that B is sourceless and hence divergenceless: ∇ · B = 0. This equation of motion is found from the Euler-Lagrange equation Since L does not depend onĊ 0 we can immediately write that ∂L/∂Ċ 0 = 0. We can also show that where, to arrive at the third line, we have integrated by parts and, to arrive at the final line, we have used that ∂C 0 (x )/∂C 0 (x) = δ(x − x ). This completes the proof. Note that treating C 0 as an independent coordinate results in a singular Hessian and so a non-invertible Lagrangian. It is better to not treat C 0 as a coordinate, and instead impose ∇ · B = 0 as a primary constraint on the Lagrangian. This is identical to Gauss's law ∇ · E = ρ and A 0 in the conventional theory. We start with the derivative with respect to C, for which we find that and note that this is also the conjugate momentum of the field. The other derivative is This is slightly more complicated so we will take each term individually. The D 2 term is evaluated using D = ∇ × C and the identity ∂/∂W (∇ × W) 2 = 2∇ × ∇ × W for any vector field W. The D · P term is evaluated as For the first term of Eqn. (SM 40) we find that

C. Lorentz force
The proceeding derivation of the Lorentz force is very lengthy. It can be bypassed by noting that our Hamiltonian in Eqn. (23) is identical to the multipolar Hamiltonian in [5] when written in terms of the physical fields. Since the multipolar gauge Hamiltonian is equivalent to the Coulomb gauge Hamiltonian, and it is readily proven that the Coulomb gauge Hamiltonian reproduces the Lorentz force through Heisenberg's equation, so must the C-field Hamiltonian. This argument can only be made for the Lorentz force derivation because this does not depend on the choice of potentials to characterise the physical fields. Nevertheless, for concreteness here we will derive the Lorentz force explicitly from the C-field Lagrangian. The It can be shown that Therefore, which is also equal to the conjugate momentum p in Eqn. (21a) if V = V 0 = 0. We then take the total time derivative of this to find d dt The position derivative separates into electric and magnetic parts and the terms involving partial derivatives of D with respect to r are zero. To arrive at the last line we have used the identity a × (∇ r × b) =∇ b (a · b) − (a · ∇ r ) b where the gradient operator with a tilde is the Feynman subscript notation, i.e.∇ b (a · b) = a · ∂b ∂r x , a · ∂b ∂r y , a · ∂b ∂r z .
(SM 51) Likewise the magnetic part is and partial derivatives of H with respect to r are zero. The electric part evaluates to and the magnetic part is Before substituting into the Euler-Lagrange equation we note that by using Faraday's law we can rewrite the second term in Eqn. (SM 48) as where we have also written B (λr) = d 3 x B(x)δ(x − λr). After some algebra we can then prove that From here the Lorentz force is derived by noting that because of the non-existence of magnetic monopoles and that the following identities hold These identities can be proven using the Fourier representation of the delta function [5]. Using these we find that mr = −e [E(r) +ṙ × B(r)] as required.

D. Hamiltonian
The Hamiltonian is derived as is standard by where q i and p i are the generalised coordinates and canonical momenta. For this theory, the canonical momenta of the matter and C-field are given in Eqns. (21) with V = V 0 = 0. We find that Focusing first on the terms within the integral, we can rewrite B ·Ċ = B · H − B · ∇C 0 . The term B · ∇C 0 then vanishes after integration by parts and enforcing ∇ · B = 0. Since the canonical field momentum is B, we want to eliminate H in favour of B to make the quantisation process easier. Therefore, writing H + M = B/µ 0 we arrive at Using the definitions of θ and Φ 0 from the main text we are then able to show that Finally, after combining this with the terms outside the integral using mṙ = p + Φ 0 , we have (SM 68) Also, recall the definitions P = P + P V , M = M − M V , where we have written the implicit change in C (C 0 ) due to the change in light/matter partition explicitly as C → C (C 0 → C 0 ). It is clear that the field canonical momentum will not change from the V = V 0 = 0 case, because it is equal to the physical field B. Therefore we trivially find that For the matter canonical momentum we must evaluate To do so we recall some algebraic relations used in the V = V 0 = 0 case: We can expand the integral as Using Eqns. (SM 75) we can show that Then, using standard vector calculus identities we find that the remaining terms can be written as Combining these, we find that as in Eqn. (21a) in the main text.

B. Hamiltonian
As usual we write the Hamiltonian in terms of the coordinates and momenta as in Eqn. (SM 61) We then define X = Φ 0 + Φ and note the following three relations: • Using the definition˙ C + ∇ C 0 = µ −1 0 B − M and that ∇ · B = 0 leads to • The terms outside the volume integral can be written as Inserting these three relations into Eqn. (SM 82) leads to The third term in Eqn. (SM 86) is a result of not being able to invert the matter canonical momentum for an expression forṙ in terms of p. We therefore impose a constraint on the choice of V and V 0 such that: This extra constraint on the transformation follows from being in general a complicated function ofṙ. Therefore, it is not possible obtain forṙ in terms of p by inverting the canonical momentum Eqn. (21a). As an example of an allowed transformation, consider V = 1 2 r × x and V 0 = 0. For this transformation, we find that Φ = 1 2 d 3 x x × B which means that we can writeṙ as a function of p by inverting Eqn. (21a). Importantly, this also means that the last term of Eqn. (SM 86) vanishes, and we arrive at Eqn. (23) in the main text. We have used various estimates when the dipole size has not been reported. In particular:

SUPPLEMENTARY MATERIAL 5: TABLE OF EXPERIMENTAL DATA POINTS
• The dipole size x 10 is approximately a couple of Bohr radii a 0 for atomic species, whereas the effective Bohr radius a eff 1 µm is used for Rydberg atoms.
• For quantum dots and quantum wells, the physical extent of the device is a good estimate for the dipole size.
• The estimates become slightly more complicated for superconducting circuits, but an approximate value can be extracted in two different ways, depending on the values reported: First from the definition of the dipole moment x 10 |d|/e with |d| = g/E vac , (SM 89) ⇒ x 10 g/(eE vac ), (SM 90) where g is the interaction strength, and E vac is the electric field amplitude of the vacuum fluctuations. Electric field fluctuations can in turn be estimated by E vac (V vac /length of relevant region), with V vac being the fluctuations in the electric potential. Alternatively, it is possibly to define an analogous Bohr radius a eff through the resonance frequency, as in ω 10 = 2πc/λ 10 = 2πc α 4πa eff , (SM 91) ⇒ a eff αc 2ω 10 , (SM 92) where α is the fine-structure constant and where we have interpreted the transition wavelength λ 10 = 4πa eff /α as an effective Bohr wavelength.
• In dye-filled polariton cavities, the spatial extent of the dye molecule acts as a bound of the dipole size.
• In the case of electron cyclotron resonances, we can estimate the dipole size as where l 0 = /(eB) is the magnetic length for a given magnetic field strength B, and ν = ρ 2deg 2πl 2 0 is the filling factor in the system with ρ 2deg being the electron density. Note that in some cases, the filling factor is directly reported. Figure SM 1 is the same calculation as for Figure 2 in the main text, however, the splitting between the two lowest lying dipole levels is no longer kept resonant with the single photon mode ν. Instead, we fix the length of the infinite square well such that the lowest lying dipole transition, ω 10 = 1 − 0 = 1 eV is held constant. Outwith resonance, the light-matter coupling ratios with ω 10 and ν in either gauge is not the same, and so the contours along the diagonal-axis and x-axis are different. Just as in Figure 2, even without resonance the C-field calculations yield more accurate results for a given coupling strength than the Coulomb gauge, and also converge onto the exact answer quicker as the number of dipole levels in the calculation is increased. Note also that it is for large dipoles in the strong-coupling limit that the C-field calculations become inaccurate, i.e. in the scenario when we expect the description of the polarisation field in terms of dipole degrees of freedom to become non-trivial.

SUPPLEMENTARY MATERIAL 6: ACCURACY OF FEW LEVEL TRUNCATION WITHOUT RESONANCE
FIG. SM 1. As in Figure 2 in the main text, we plot the relative error in calculating the energy spacing between the two lowest eigenstates of the full Coulomb gauge and C-field Hamiltonians within the two and three dipole-level approximations, using an infinite square well potential as Uext. The difference is that in this plot the lowest lying dipole transition is kept fixed at 1 eV, and so, is not on resonance with the mode frequency which is varied along the y-axes. For all other details refer to Figure ??.