Uniqueness of the Phase Transition in Many-Dipole Cavity Quantum Electrodynamical Systems

TL;DR

This paper demonstrates a unique superradiant phase transition in many-dipole cavity QED systems, clarifying the conflicting no-go theorems by emphasizing gauge invariance and the definition of radiation.

Contribution

It introduces an arbitrary gauge approach to show the existence of a single phase transition and clarifies the role of gauge choice in superradiance in cavity QED.

Findings

A unique phase transition occurs in many-dipole cavity QED systems.

Gauge choice influences the inclusion of polarisation in the radiative subsystem.

The paradox of no-go and counter no-go theorems is resolved by different radiation definitions.

Abstract

The possibility of a superradiant phase transition in light-matter systems is the subject of much debate, due to numerous apparently conflicting no-go and counter no-go theorems. Using an arbitrary gauge approach we show that a unique phase transition does occur in archetypal many-dipole cavity QED systems, and that it manifests unambiguously via a macroscopic gauge-invariant polarisation. We find that the gauge choice controls the extent to which this polarisation is included as part of the radiative quantum subsystem and thereby determines the degree to which the abnormal phase is classed as superradiant. This resolves the long-standing paradox of no-go and counter no-go theorems for superradiance, which are shown to refer to different definitions of radiation.

Equations

$${\mathbf{p}}_{\mu \alpha }$$

$${\mathbf{\Pi }}_{\alpha }(\mathbf{x})$$

$$H^{\alpha ,2}=$$

$$-\frac{{\mathcal{C}}_{\alpha }}{N}{\left(J_{\alpha }^{+}+J_{\alpha }^{-}\right)}^2-i\frac{g_{\alpha }^{\prime }}{\sqrt{N}}(J_{\alpha }^{+}-J_{\alpha }^{-})(c_{\alpha }^{†}+c_{\alpha })$$

$$+i\frac{g_{\alpha }}{\sqrt{N}}(J_{\alpha }^{+}+J_{\alpha }^{-})(c_{\alpha }^{†}-c_{\alpha }),$$

$${{\tilde{\omega }}_m^{\alpha }}^2$$

$${\omega }_m\ge 4{\mathcal{C}}_{\alpha }=2\rho d^2(1-{\alpha }^2).$$

$$H_{\mathrm{th}}^{\alpha ,2,\mathrm{i}}=$$

$$2{E_{\alpha \pm }^{\mathrm{n}}}^2=$$

$$\displaystyle+16[{\tilde{\omega}}_{m}^{\alpha}{\tilde{g}}_{\alpha}^{\prime}+\omega_{\alpha}{\tilde{g}}_{\alpha}][{\tilde{\omega}}_{m}^{\alpha}{\tilde{g}}_{\alpha}+\omega_{\alpha}{\tilde{g}}_{\alpha}^{\prime}]\bigg{)},$$

$${\omega }_m\left({\omega }_m-2\rho d^2\right)\left({\omega }_{\alpha }^2-2{\omega }_m\rho d^2(1-\alpha {)}^2\right)\ge 0.$$

$${\omega }_m\ge 2\rho d^2.$$

$$2{E_{\alpha \pm }^{\mathrm{a}}}^2=$$

$$\displaystyle+16[\underaccent{\tilde}{\omega}_{m}^{\alpha}\underaccent{\tilde}{g}_{\alpha}^{\prime}+\omega_{\alpha}\underaccent{\tilde}{g}_{\alpha}^{\prime}][\underaccent{\tilde}{\omega}_{m}^{\alpha}\underaccent{\tilde}{g}_{\alpha}+\omega_{\alpha}\underaccent{\tilde}{g}_{\alpha}^{\prime}]\bigg{)},$$

$$\left([2\rho d^2{]}^2-{{\omega }_m}^2\right)\left({\omega }_{\alpha }^2-{{\omega }_m}^2(1-\alpha {)}^2\right)\ge 0.$$

$${\Pi }_{\alpha ,\mathrm{th}}^{\mathrm{a}}=\alpha \rho d\sqrt{1-{\tau }^2}=\frac{\alpha }{2d}\sqrt{{\left(2\rho d^2\right)}^2-{{\omega }_m}^2}.$$

$$P_{\mathrm{T},\mathrm{th}}^{\mathrm{a}}=-\rho d\sqrt{1-{\tau }^2}$$

$$\displaystyle{1\over N\omega}{d^{2}G_{s}^{\alpha,2}\over d\eta^{2}}\bigg{|}_{\rm th,a}=-\frac{\omega_{m}}{\omega}{d^{2}\over d\eta^{2}}{(1-\tau)^{2}\over 4\tau}$$

$$\mathbf{V}={\mathbf{V}}_{\mathrm{T}}+{\mathbf{V}}_{\mathrm{L}},$$

$$\nabla \cdot {\mathbf{V}}_{\mathrm{T}}(\mathbf{x})=0,$$

$$\nabla \times {\mathbf{V}}_{\mathrm{L}}(\mathbf{x})=0.$$

$$\int d^{3}x\mathbf{V}(\mathbf{x})\cdot \nabla f(\mathbf{x})=-\int d^{3}xf(\mathbf{x})\nabla \cdot \mathbf{V}(\mathbf{x}).$$

$$\int d^{3}x{\mathbf{V}}_{\mathrm{T}}(\mathbf{x})\cdot {\mathbf{U}}_{\mathrm{L}}(\mathbf{x})=0$$

$$\rho (\mathbf{x})$$

$$\mathbf{J}(\mathbf{x})$$

$$\mathbf{P}(\mathbf{x})=-\int d^3{x}^{\prime }\mathbf{g}(\mathbf{x},{\mathbf{x}}^{\prime })\rho ({\mathbf{x}}^{\prime })=-\sum _{\mu =1}^N\int d^3{x}^{\prime }\mathbf{g}(\mathbf{x},{\mathbf{x}}^{\prime }){\rho }_{\mu }({\mathbf{x}}^{\prime })\equiv \sum _{\mu =1}^N{\mathbf{P}}_{\mu }(\mathbf{x})$$

$$\nabla \cdot \mathbf{g}(\mathbf{x},{\mathbf{x}}^{\prime })=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })$$

$${\mathbf{g}}_{\mathrm{L}}(\mathbf{x},{\mathbf{x}}^{\prime })=-\nabla \frac{1}{4\text{uppi}|\mathbf{x}-{\mathbf{x}}^{\prime }|}.$$

$$L=\frac{1}{2}\sum _{\mu =1}^N\left[m_{\mathrm{p}}{\dot{\mathbf{r}}}_{\mathrm{p}\mu }^2+m_{\mathrm{e}}{\dot{\mathbf{r}}}_{\mathrm{e}\mu }^2\right]-V_{\mathrm{ext}}-\int d^{3}x\left[j^a(\mathbf{x})A_a(\mathbf{x})+\frac{1}{4}{F}_{ab}(\mathbf{x})F^{ab}(\mathbf{x})\right]$$

$$\chi (\mathbf{x})=\int d^3{x}^{\prime }\mathbf{g}({\mathbf{x}}^{\prime },\mathbf{x})\cdot {\mathbf{A}}_{\mathrm{T}}({\mathbf{x}}^{\prime })\equiv \int d^3{x}^{\prime }{\mathbf{g}}_{\mathrm{T}}({\mathbf{x}}^{\prime },\mathbf{x})\cdot {\mathbf{A}}_{\mathrm{T}}({\mathbf{x}}^{\prime })$$

$$A_0$$

$$\mathbf{A}$$

$${\varphi }_{\mathrm{Coul}}(\mathbf{x})=\int d^3{x}^{\prime }\frac{\rho ({\mathbf{x}}^{\prime })}{4\pi |\mathbf{x}-{\mathbf{x}}^{\prime }|}.$$

$$L^{\prime }=L-\frac{d}{dt}\int d^{3}x\mathbf{P}(\mathbf{x})\cdot \mathbf{A}(\mathbf{x}).$$

$$L^{\prime }=L_0-\frac{d}{dt}\int d^{3}x{\mathbf{P}}_{\mathrm{T}}(\mathbf{x})\cdot {\mathbf{A}}_{\mathrm{T}}(\mathbf{x})$$

Full text

Superradiance was originally described by Dicke Dicke (1954), and since then it has received a great deal of attention (see e.g. Kirton et al. (2019) for a recent introduction). A superradiant phase of a light-matter system is one in which a macroscopic number of photons arises due to the interaction between many dipoles. The possibility of such a phase transition within the Dicke model was first recognised some time ago Hepp and Lieb (1973a); Wang and Hioe (1973). Later, seminal contributions were made in the connection with quantum chaos Holstein and Primakoff (1940); Emary and Brandes (2003a, b). The topic now includes extended Dicke models Carmichael et al. (1973); Hioe (1973); Pimentel and Zimerman (1975); Emeljanov and Klimontovich (1976); Sung and Bowden (1979), driven and open systems, semi-classical descriptions Grimsmo and Parkins (2013); Klinder et al. (2015); Gegg et al. (2018); Kirton and Keeling (2018); Peng et al. (2019), and artificial QED systems Yamanoi (1979); Haug and Koch (1994); Lee and Johnson (2004); Nataf and Ciuti (2010); Viehmann et al. (2011); Todorov and Sirtori (2012); Leib and Hartmann (2014); Bamba and Ogawa (2014); Bamba et al. (2016); Bamba and Imoto (2017); Jaako et al. (2016); De Bernardis et al. (2018a).

One of the most controversial aspects of theoretical studies has been the validity of so-called “no-go theorems”, which prohibit a superradiant phase and are proved in the Coulomb-gauge. The original no-go theorem Rzażewski et al. (1975) actually prohibits a phase transition of any kind, but neglects direct electrostatic interactions, whose presence are a defining feature of a correct Coulomb-gauge model. This theorem, and variants thereof, have been both refuted and confirmed in numerous subsequent works Kudenko et al. (1975); Rzażewski et al. (1976); Emeljanov and Klimontovich (1976); Knight et al. (1978); Bialynicki-Birula and Rza¸żewski (1979); Yamanoi (1979); Sung and Bowden (1979); Rzażewski and Wódkiewicz (1991); Keeling (2007); Nataf and Ciuti (2010); Vukics and Domokos (2012); Vukics et al. (2014); Bamba and Ogawa (2014); Tufarelli et al. (2015); Grießer et al. (2016); Andolina et al. (2019). It has been suggested that where natural atomic systems admit a no-go theorem certain artificial atomic systems do not Nataf and Ciuti (2010) (though see also Viehmann et al. (2011); Ciuti and Nataf (2012)). However, in the multipolar-gauge the superradiant phase transition also appears to be automatically recovered for conventional cavity QED systems Vukics and Domokos (2012); Vukics et al. (2014).

Further permutations of these results are available. For example, if explicit dipole-dipole interactions that are not naturally present are added into the multipolar-gauge description, then a no-go theorem re-emerges Bialynicki-Birula and Rza¸żewski (1979); Haug and Koch (1994); De Bernardis et al. (2018b); Stefano et al. (2019). A very recent contribution Andolina et al. (2019) argues without the two-level approximation that a superradiant phase is impossible, but this treatment considers only the radiative quantum subsystem and is again proved in the Coulomb gauge. If, rather than just the radiative subsystem, one also considers variations in the electrostatic interactions that are present within the Coulomb-gauge, then an apparently different ferroelectric phase transition is predicted. This, however, does not lead to superradiance Keeling (2007). Thus, despite numerous contributions spanning several decades, the occurrence and nature of the phase transition in generic many-emitter light-matter systems, and how this relates to the choice of gauge, are fundamental questions whose answers remain unclear, yet still highly relevant Mazza and Georges (2019); Andolina et al. (2019); Nataf et al. (2019); Guerci et al. (2020); Andolina et al. (2020).

Here we resolve these fundamental issues by proving that a unique physical phase transition does occur in generic many-dipole cavity QED systems and that the abnormal phase of the system is unambiguously signalled by a macroscopic average of the gauge-invariant transverse polarisation field ${\bf P}_{\rm T}$. This equals the longitudinal electric field ${\bf E}_{\rm L}$ except at the point-dipole positions themselves. Crucial to the resolution provided is the recognition that QED subsystems are gauge-relative, meaning that each gauge provides different gauge-invariant definitions of the light and matter subsystems. Whether the abnormal phase is characterised as ferroelectric or as superradiant depends on the extent to which ${\bf E}_{\rm L}={\bf P}_{\rm T}$ is included within the radiative quantum subsystem and this is controlled by the gauge choice. We thereby show that the different viewpoints provided by different gauges are not contradictory, but in fact equivalent, as required. In particular, correct no-go statements such as in Ref. Andolina et al. (2019) are reconciled with correct counter no-go statements such as in Ref. Vukics and Domokos (2012). Such results are found to be different ways of viewing the same phenomenon in terms of physically different quantum subsystems. By converting the apparent gauge non-invariance of the phase transition into a proof of gauge-invariance, our results resolve the associated long-standing controversies.

A related but separate point is that level truncation of material dipoles causes a breakdown of gauge-invariance De Bernardis et al. (2018b); Stokes and Nazir (2019a); Stefano et al. (2019); Roth et al. (2019). Using numerical results for finite numbers of dipoles we show that accurate two-level model predictions can be identified. It is reasonable to conclude that the same two-level truncation will be accurate in the thermodynamic limit. Thus, arbitrary-gauge QED is also capable of eliminating any further quantitative ambiguity resulting from the use of material two-level truncation.

We begin by deriving an arbitrary gauge Dicke Hamiltonian. We adopt a general formulation of QED in which the gauge is selected by a real parameter $\alpha$. We consider $N$ identical electric dipoles each described by a classical centre-of-mass position ${\bf R}_{\mu}$ and a dipole moment operator ${\hat{\bf d}}_{\mu}=-e{\bf r}_{\mu}$. The dipoles interact with a common electromagnetic field described by transverse-electric and magnetic fields ${\bf E}_{\rm T}$ and ${\bf B}$ respectively. We obtain the Hamiltonian for the system from first principles 111Please see the Supplemental Material, which includes Refs. Stokes and Nazir (2018) and Williamson (1936)., which can be written in the gauge-invariant form Stokes and Nazir (2019a) $H=E_{\rm matter}+E_{\rm field}$, where $E_{\rm matter}:=\sum_{\mu=1}^{N}{1\over 2}m{\dot{\bf r}}_{\mu}^{2}+V+V_{\rm dip}$ and $E_{\rm field}:={1\over 2}\int d^{3}x\,\left[{\bf E}_{\rm T}({\bf x})^{2}+{\bf B}({\bf x})^{2}\right]$. Here $V$ denotes the total intra-dipole potential, and $V_{\rm dip}$ denotes the inter-dipole electrostatic energy. The $\alpha$-dependent canonical momenta are found to be

[TABLE]

where ${\bf A}$ is the gauge-invariant transverse vector potential such that ${\bf E}_{\rm T}=-{\dot{\bf A}}$ and ${\bf P}_{{\rm T}\alpha}$ is the $\alpha$-gauge transverse polarisation given by ${\bf P}_{{\rm T}\alpha}({\bf x})=\alpha{\bf P}_{{\rm T}}({\bf x})$, with ${\bf P}_{\rm T}({\bf x})=\sum_{\mu=1}^{N}{\hat{\bf d}}_{\mu}\cdot\delta^{\rm T}({\bf x}-{\bf R}_{\mu})$. The canonical commutation relations are $[r_{\mu,i},p_{\nu,j}]=i\delta_{\mu\nu}\delta_{ij}$ and $[A_{{\rm T},i}({\bf x}),\Pi_{{\rm T},j}({\bf x}^{\prime})]=\delta_{ij}^{\rm T}({\bf x}-{\bf x}^{\prime})$. All other commutators between canonical operators vanish. The canonical momenta of different gauges are unitarily related via ${\bf X}_{\alpha}=R_{\alpha\alpha^{\prime}}{\bf X}_{\alpha^{\prime}}R_{\alpha^{\prime}\alpha}$, where ${\bf X}={\bf p},\,{\bf\Pi}$ and $R_{\alpha\alpha^{\prime}}=\exp\left(i(\alpha-\alpha^{\prime})\sum_{\mu=1}^{N}{\hat{\bf d}}_{\mu}\cdot{\bf A}\right)$.

We now restrict our attention to a single cavity mode with volume $v$, frequency $\omega$ and unit polarisation vector ${\bm{\varepsilon}}$, described by bosonic operators $a_{\alpha},~{}a_{\alpha}^{\dagger}$ with $[a_{\alpha},a_{\alpha}^{\dagger}]=1$. The restriction is imposed consistently on all fields including the transverse delta-function $\delta^{\rm T}$. This eliminates the need to regularise ${\bf P}_{\rm T}$ Vukics et al. (2015), and ensures that the transverse commutation relation for the canonical fields is preserved. The fundamental kinematic relations given by Eqs. (54) and (55) are therefore also preserved. In order to obtain a Dicke Hamiltonian we next take the limit of closely spaced dipoles around the origin; ${\bf R}_{\mu}\approx{\bf 0}$ and we approximate the dipoles as two-level systems. Further details of all approximations used are given in Note (1). We introduce the collective operators $J_{\alpha}^{i}=\sum_{\mu=1}^{N}\sigma^{i}_{\mu\alpha}$, with $i=\pm,z$, where $\sigma^{\pm}_{\mu\alpha}$ are the raising and lowering operators of the $\mu$’th two-level dipole and $\sigma^{z}_{\mu\alpha}=[\sigma^{+}_{\mu\alpha},\sigma^{-}_{\mu\alpha}]/2$. We also introduce cavity bosonic operators $c_{\alpha}$ and $c_{\alpha}^{\dagger}$, which incorporate both the bare cavity energy and the ${\bf A}^{2}$-term that results when Eq. (54) is substituted into the energy $H$ [Supplementary Eqs. (72), (73)]. The resulting arbitrary-gauge Dicke Hamiltonian is

[TABLE]

where $\omega_{m}=\epsilon_{1}-\epsilon_{0}$, ${\omega_{\alpha}}^{2}=\omega^{2}+e^{2}(1-\alpha)^{2}\rho/m$, ${\cal C}_{\alpha}:=\rho d^{2}(1-\alpha^{2})/2$, $g_{\alpha}^{\prime}:=(1-\alpha)\omega_{m}d\sqrt{\rho/(2\omega_{\alpha})}$, and $g_{\alpha}:=\alpha d\sqrt{\rho\omega_{\alpha}/2}$, with $d:={\bm{\varepsilon}}\cdot{\bf d}$. Here $\rho=N/v$ remains finite in the thermodynamic limit $N\to\infty$, $v\to\infty$. Although the non-truncated Hamiltonian $H$ is unique, we now have a continuous infinity of Dicke Hamiltonians $H^{\alpha,2}$ such that $H^{\alpha,2}$ and $H^{\alpha^{\prime},2}$ are not equal when $\alpha\neq\alpha^{\prime}$ De Bernardis et al. (2018b); Stokes and Nazir (2019a); Stefano et al. (2019); Roth et al. (2019). This breaking of gauge-invariance will turn out not to be a barrier to eliminating all ambiguities regarding the occurrence and nature of a quantum phase transition.

To take the thermodynamic limit we use a Holstein-Primakoff map defined by $J_{\alpha}^{z}=b_{\alpha}^{\dagger}b_{\alpha}-{N/2}$, $J_{\alpha}^{+}=b_{\alpha}^{\dagger}\sqrt{N-b_{\alpha}^{\dagger}b_{\alpha}}$, and $J_{\alpha}^{-}=(J_{\alpha}^{+})^{\dagger}$, where $[b_{\alpha},b_{\alpha}^{\dagger}]=1$ Holstein and Primakoff (1940); Emary and Brandes (2003a, b). The Hamiltonian obtained by substituting these expressions into Eq. (Uniqueness of the Phase Transition in Many-Dipole Cavity Quantum Electrodynamical Systems) is denoted $H_{\rm th}^{\alpha,2}$. We first consider the material part of $H_{\rm th}^{\alpha,2}$, which can be written $H_{{\rm th,}m}^{\alpha,2}={\tilde{\omega}}_{m}^{\alpha}l^{\dagger}_{\alpha}l_{\alpha}+{1\over 2}({\tilde{\omega}}_{m}^{\alpha}-\omega_{m})$ where $[l_{\alpha},l_{\alpha}^{\dagger}]=1$ and

[TABLE]

The mode operators $l_{\alpha},~{}l_{\alpha}^{\dagger}$ are related to $b_{\alpha}$ and $b_{\alpha}^{\dagger}$ by a local Bogoliubov transformation that incorporates the contribution $V_{\rm dip}$ (see Eqs. (79) and (80) in Note (1)). This results in the renormalised frequency ${\tilde{\omega}}_{m}^{\alpha}$ in Eq. (4). Reality of ${\tilde{\omega}}_{m}^{\alpha}$ requires that

[TABLE]

When the electrostatic interaction strength ${\cal C}_{\alpha}$ is large enough this inequality may be violated signalling a phase transition. We refer to this transition as ferroelectric, because it is completely independent of the radiative mode. Inequality (5) generalises the result of Keeling obtained when $\alpha=0$ (Coulomb gauge) Keeling (2007). Violation of inequality (5) cannot occur in the multipolar gauge $\alpha=1$, which does not therefore admit a purely ferroelectric phase. In what follows this finding will be reconciled with our claim that a unique phase transition is predicted within all gauges. We show further that only in the Coulomb gauge does the phase transition appear purely ferroelectric.

We now consider the thermodynamic-limit of the total Hamiltonian, which is Note (1)

[TABLE]

where the superscript $\cdot^{\rm i}$ is either ${\rm i}={\rm n}$ for normal-phase, or ${\rm i}={\rm a}$ for abnormal-phase. The polariton operators $f^{\rm i}_{\alpha},\,c^{\rm i}_{\alpha}$ are bosonic satisfying $[f^{\rm i}_{\alpha},{f^{\rm i}_{\alpha}}^{\dagger}]=1=[c^{\rm i}_{\alpha},{c^{\rm i}_{\alpha}}^{\dagger}]$ with all other commutators vanishing. In the normal phase, ${\rm i}={\rm n}$, the zero-point constant in Eq. (6) is $C^{\rm n}=N\epsilon_{0}+\left(\rho d^{2}-\omega_{m}\right)/2$ and the polariton energies are

[TABLE]

where ${\tilde{g}}_{\alpha}=g_{\alpha}\sqrt{\omega_{m}/{\tilde{\omega}}_{m}^{\alpha}}$ and ${\tilde{g}}^{\prime}_{\alpha}=g_{\alpha}^{\prime}\sqrt{{\tilde{\omega}}_{m}^{\alpha}/\omega_{m}}$. The coupling strength at which the lower polariton energy $E^{\rm n}_{\alpha-}$ is no longer real signals the onset of the abnormal phase and the breakdown of $H_{\rm th}^{\alpha,2,{\rm n}}$. Reality of $E^{\rm n}_{\alpha-}$ requires that

[TABLE]

From the Thomas-Reiche-Kuhn (TRK) inequality ${e^{2}/m}\geq 2\omega_{m}d^{2}$ it follows that $\omega_{\alpha}^{2}\geq 2\omega_{m}\rho d^{2}(1-\alpha)^{2}$. Therefore, by inequality (8) $E_{\alpha-}$ is real if and only if

[TABLE]

This simple gauge-invariant result defines the normal phase. Inequality (9) is stronger than inequality (5), so ${\tilde{\omega}}_{m}^{\alpha}$ in Eq. (4) is also real when inequality (9) is satisfied.

The Hamiltonian $H_{\rm th}^{\alpha,2,\rm a}$ for the abnormal phase takes over from $H_{\rm th}^{\alpha,2,\rm n}$ when inequality (9) is violated. It is obtained as the thermodynamic limit of $H^{\alpha,2}$ written, via the Holstein-Primakoff map, in terms of displaced modes $f_{\alpha}$ and $c^{\prime}_{\alpha}$ such that $b_{\alpha}=f_{\alpha}-\sqrt{\beta_{\alpha}}$ and $c_{\alpha}=c^{\prime}_{\alpha}+i\sqrt{\gamma_{\alpha}}$, where $\beta_{\alpha}=\beta:=N(1-\tau)/2$ and $\gamma_{\alpha}={Ng_{\alpha}^{2}}\left(1-\tau^{2}\right)/\omega_{\alpha}^{2}$ are of order $N$, with $\tau={\omega_{m}/(2\rho d^{2})}$. Note that $\beta_{\alpha}=\beta$ is $\alpha$-independent indicating that the “material” mode is always displaced by the same macroscopic quantity. On the other hand, $\gamma_{\alpha}$ is $\alpha$-dependent, so the extent to which the “radiative” mode is displaced depends on the chosen definition of radiation. In particular, $\gamma_{0}=0$, so in the Coulomb gauge only the material mode is displaced. In the abnormal-phase, ${\rm i}={\rm a}$, the zero-point constant in Eq. (6) is $C^{\rm a}=N\left[\epsilon_{0}-{\omega_{m}}(1-\tau)^{2}\right]/(4\tau)-\rho d^{2}/2$ and the polariton energies are

[TABLE]

where {\underaccent{\tilde}{\omega}_{m}^{\alpha}}^{2}={{\omega_{m}}^{2}}\left[1-(1-\alpha^{2})\tau^{2}\right]/\tau^{2}, while \underaccent{\tilde}{g}^{\prime}_{\alpha}=g_{\alpha}^{\prime}\sqrt{\tau\underaccent{\tilde}{\omega}_{m}^{\alpha}/\omega_{m}}, and \underaccent{\tilde}{g}_{\alpha}=g_{\alpha}\sqrt{\tau\omega_{m}/\underaccent{\tilde}{\omega}_{m}^{\alpha}}. The material frequency \underaccent{\tilde}{\omega}_{m}^{\alpha} is real provided $(2\rho d^{2})^{2}\geq{\omega_{m}}^{2}(1-\alpha^{2})$ and the lower polariton energy $E_{\alpha-}^{\rm a}$ is real provided

[TABLE]

In the abnormal phase we have $2\rho d^{2}\geq\omega_{m}$ implying that $\omega_{\alpha}^{2}-{\omega_{m}}^{2}(1-\alpha)^{2}\geq 0$ and therefore that $E_{\alpha-}^{\rm a}$ is real. At the critical coupling point where $2\rho d^{2}=\omega_{m}$ the Hamiltonians $H_{\rm th}^{\alpha,2,\rm n}$ and $H_{\rm th}^{\alpha,2,\rm a}$ coincide. We have therefore obtained a description of the thermodynamic limit for all coupling strengths. The polariton energies constitute different two-level approximated results in each different gauge $\alpha$. This is shown in Supplementary Figure 1. However, every gauge’s approximate (Dicke) model predicts exactly one ground state phase transition occurring when $\omega_{m}=2\rho d^{2}$ and the ground state is unique within the non-truncated theory. Therefore, inequality (9) should be interpreted as predicting a unique phase transition. However, the nature of the phase transition appears to be different depending on the value of $\alpha$. In the Coulomb gauge for example, it is necessarily purely ferroelectric, whereas this is impossible in the multipolar-gauge.

The radiative classification of a unique phase transition will naturally depend on the definition of radiation and the latter is controlled by the gauge. Evidently, the subsystem gauge-relativity of QED Stokes and Nazir (2019b), is strongly exemplified by the phase transition phenomenon. To understand the physical meaning of “matter” and “radiation” in the gauge $\alpha$ we note that the total multipolar polarisation of $N$ dipoles is ${\bf P}({\bf x})=\sum_{\mu=1}^{N}{\bf d}_{\mu}\delta({\bf x}-{\bf R}_{\mu})$. Since $\nabla\cdot{\bf E}=\rho=-\nabla\cdot{\bf P}$ it follows that for ${\bf x}\neq{\bf R}_{\mu}$ we have ${\bf P}_{\rm T}={\bf E}_{\rm L}$ and therefore ${\bf\Pi}_{\alpha}=-{\bf E}_{\rm T}-\alpha{\bf E}_{\rm L}$ [cf. Eq. (55)]. Similarly, the material momentum ${\bf p}_{\alpha}$ of a dipole is given by Eq. (54) in which $q{\bf A}({\bf 0})$ is the electric dipole approximation (EDA) of $q{\bf A}({\bf r})={\bf P}_{\rm long}=\int d^{3}x{\bf E}_{\rm L\bf r}({\bf x})\times{\bf B}({\bf x})$, which is the momentum of the longitudinal field generated by $q$ at ${\bf r}$ with $\nabla\cdot{\bf E}_{\rm L\bf r}({\bf x})=q\delta({\bf x}-{\bf r})$.

As an example, one may consider the Coulomb-gauge in which “matter” is fully dressed by ${\bf E}_{\rm L}$, i.e., ${\bf p}_{0}=m{\dot{\bf r}}+{\bf P}_{\rm long}$, so “matter” as defined by ${\bf p}_{0}$ is not fully localised. Correspondingly, “radiation” is defined using the field ${\bf\Pi}_{0}=-{\bf E}_{\rm T}$ alone. In the multipolar-gauge (within the EDA) matter is completely bare, i.e., ${\bf p}_{1}=m{\dot{\bf r}}$, and therefore fully localised. “Radiation” is correspondingly defined for ${\bf x}\neq{\bf 0}$ by the local (causal) total field ${\bf\Pi}_{1}=-{\bf D}_{\rm T}=-{\bf E}_{\rm T}-{\bf E}_{\rm L}=-{\bf E}$. More generally, $\alpha$ controls how the longitudinal electric degrees of freedom are shared out, thereby controlling the balance between localisation and electrostatic dressing in defining the quantum subsystem called “matter”. “Radiation” is then defined using the canonical degrees of freedom left over.

There are noteworthy gauges in between $\alpha=0$ and $\alpha=1$, such as gauges relative to which ground state “virtual photons” are highly suppressed and for which the corresponding two-level model can sometimes offer a more accurate representation of the ground state than conventional quantum Rabi models Stokes and Nazir (2019a) (see also Note (1)). What differs between gauges is the spacetime locatisation properties of “material sources” and their dressing by virtual photons. In general the most operationally relevant definitions of the subsystems may depend on the available measurements, including their time- and length-scales. As a result, general statements about measurable photon condensation (superradiance), that are independent of experimental context, cannot be made. What can be demonstrated and is demonstrated below, is that there are no internal theoretical inconsistencies and no fundamental paradoxes. Previous no-go and counter no-go theorems refer to different definitions of radiation and so are not contradictory. They are in fact equivalent.

We now calculate the ground-state momentum ${\bf\Pi}_{\alpha}$ of radiation defined relative to gauge $\alpha$. This directly demonstrates strict equivalence of all gauges and reveals an unambiguous macroscopic manifestation of the abnormal phase. We allow the two-level truncation to be performed in an arbitrary gauge $\alpha^{\prime}$. The $\alpha^{\prime}$-gauge two-level approximation of an operator $o_{\alpha}=o_{\alpha}({\bf p}_{\alpha},{\bf\Pi}_{\alpha})$, denoted $o^{\alpha^{\prime},2}_{\alpha}$, is found by expressing $o_{\alpha}$ in terms of $\alpha^{\prime}$-gauge canonical operators followed by two-level truncation. For $\Pi_{\alpha}={\bm{\varepsilon}}\cdot{\bf\Pi}_{\alpha}$ we have $\Pi_{\alpha}^{\alpha^{\prime},2}=\Pi_{\alpha^{\prime}}-{d}(\alpha-\alpha^{\prime})(J_{\alpha^{\prime}}^{+}+J_{\alpha^{\prime}}^{-})/v$. We will see that in the thermodynamic limit the ground state value of $\Pi_{\alpha}^{\alpha^{\prime},2}$ is actually independent of $\alpha^{\prime}$, i.e. the prediction is gauge-invariant, so we return to the simpler notation $\Pi_{\rm\alpha,th}$. Using the Holstein-Primakoff representation, we find that $\Pi_{\alpha,\rm th}$ vanishes in the normal phase and in the abnormal phase is proportional to the identity. The calculation in Note (1) yields the simple result

[TABLE]

The factor of $\alpha$ in Eq. (12) is highly significant. It demonstrates that the degree of superradiance in the abnormal phase is proportional to $\alpha$, with the minimum value of zero occurring only in the Coulomb-gauge.

To demonstrate equivalence between all gauges we calculate the $\alpha$-gauge transverse polarisation $P_{\rm T\alpha}=\alpha{\bm{\varepsilon}}\cdot{\bf P}_{\rm T}=\alpha(\Pi_{0}-\Pi_{1})$, which is such that $P_{\rm T\alpha}^{\alpha^{\prime},2}=\alpha{d\over v}(J_{\alpha^{\prime}}^{+}+J_{\alpha^{\prime}}^{-})$. This quantity is also $\alpha^{\prime}$-independent in the thermodynamic limit. In the normal phase $P_{\rm T\alpha,th}$ vanishes, whereas in the abnormal phase it is found by the same method that leads to Eq. (12) to be $P_{\rm T\alpha,th}^{\rm a}=-\alpha\rho d\sqrt{1-\tau^{2}}$. Eq. (12) then yields $\Pi_{\alpha,{\rm th}}^{\rm a}=-P_{\rm T\alpha,th}^{\rm a}$, which since $-E_{\rm T,th}^{\rm a}=\Pi_{0,{\rm th}}^{\rm a}=0$, is seen to be nothing but the fundamental kinematic relation (55). This establishes consistency between all gauges. The quantity $|P_{\rm T,th}^{\alpha^{\prime},2,{\rm a}}|=|P_{\rm T\alpha,th}^{\alpha^{\prime},2,{\rm a}}/\alpha|$ provides a gauge-invariant monotonic measure of the coupling-distance past the phase transition point. Thus, independent of the gauge the onset of the abnormal phase manifests in the form of a macroscopic value of the gauge-invariant field ${\bf P}_{\rm T}$;

[TABLE]

which is plotted in Supplementary Figure 2. Within the present simplified Dicke-type treatment the field ${\bf P}_{\rm T}$ is independent of spatial position ${\bf x}$, but at a more fundamental level ${\bf P}_{\rm T}$ coincides with the longitudinal electric field ${\bf E}_{\rm L}$ away from the dipole positions, i.e., for ${\bf x}\neq{\bf R}_{\mu}$. Whether one considers ${\bf E}_{\rm L}={\bf P}_{\rm T}$ to be “material” or “radiative” determines whether one calls the phase transition “purely ferroelectric” or “superradiant”, and this in turn is determined by the gauge choice as discussed earlier.

We finally consider a concrete example. We assume that each dipole has canonical operators pointing along ${\bm{\varepsilon}}$ and a double-well potential $V(\theta,\phi)=-\theta r^{2}/2+\phi r^{4}/4$ where $\theta$ and $\phi$ control the shape of the double-well. The Hamiltonian of each dipole is therefore $H_{m}^{\alpha}={{\cal E}\over 2}\left(-\partial_{\zeta}^{2}-\beta\zeta^{2}+{\zeta^{4}\over 2}\right)$ De Bernardis et al. (2018b) where we have defined $\zeta=r/r_{0}$ with $r_{0}=(1/[m\phi])^{1/6}$, along with ${\cal E}=1/(mr_{0}^{2})$ and $\beta=\theta mr_{0}^{4}$. We also define the gauge-invariant dimensionless coupling parameter $\eta=(e/\omega)\sqrt{\rho/m}$. The parameters $e,\,m$ and $\rho$ can now be eliminated in favour of ${\cal E},\beta$ and $\eta$.

To demonstrate that it is possible to obtain accurate two-level model predictions and to show a clear precursor to the phase transition, in Fig. 1 we consider the normalised second derivative of the (shifted) ground energy Emary and Brandes (2003b). In the abnormal phase of the thermodynamic limit this is given by

[TABLE]

where $G_{s}=G-\rho d^{2}/2$ is the ground energy $G$ shifted by the coupling-dependent term $\rho d^{2}/2$ in Eq. (Uniqueness of the Phase Transition in Many-Dipole Cavity Quantum Electrodynamical Systems). We choose $\beta=3.3$, which provides a highly anharmonic single-dipole spectrum such that $(\epsilon_{2}-\epsilon_{0})/\omega_{m}\approx 36$. The two-level truncation within the multipolar gauge $\alpha=1$ is subsequently found to be accurate in predicting low energy properties. This was first confirmed in the case of the Rabi model $N=1$ in Ref. De Bernardis et al. (2018b). The accuracy actually increases with $N$ and convergence of exact (gauge-invariant, no two-level truncation) and approximate predictions already occurs at $N=4$. The situation may change if the double-well is parameterised differently such that the multipolar truncation is not optimal Stokes and Nazir (2019a), see also additional analysis in Note (1).

The situation may also change if additional cavity modes are taken into account Mu oz et al. (2018); Roth et al. (2019). In particular, the multipolar-gauge coupling scales as $\sqrt{\omega}$ such that the single-mode approximation appears least favourable in this gauge, and has been shown to breakdown in the ultrastrong-coupling regime Mu oz et al. (2018). To incorporate some of the effects of non-resonant modes within a Dicke-type model a formal procedure of adiabatic elimination can be used and this also has the advantage of enabling an exploration of more diverse dipolar geometries De Bernardis et al. (2018a). Nevertheless, for our purpose of determining whether a physical phase transition can be supported by systems describable using a Dicke model and on understanding its relationship to the choice of gauge, the single-mode restriction is sufficient, because the qualitative behaviour of the thermodynamic limit of the single-mode Dicke model is known to carry over to the multi-mode case Hepp and Lieb (1973b); Pimentel and Zimerman (1975). The extension to general dipolar arrangements, and to more sophisticated cavity models warrants further study.

We have shown that a unique physical phase transition can occur in simple many-dipole cavity QED systems. We have resolved all ambiguities pertaining to the choice of gauge by determining both the origin and properties of the phase transition in terms of any gauge’s definitions of the quantum subsystems, and by demonstrating equivalence between all gauge choices. We have shown that the original “no-go theorem” Rzażewski et al. (1975) does not apply, and also that one need not look beyond ordinary cavity QED in order to find systems supporting a superradiant phase transition. A no-go theorem for ground state superradiance occurs for, and only for, the Coulomb-gauge definition of radiation. We have shown that although the two-level approximation ruins the gauge-invariance of the theory, unambiguous predictions can still be obtained. The framework developed here should be straightforwardly extendable to artificial solid-state and superconducting systems, as well as to driven and dissipative systems. This will elucidate both qualitatively and quantitatively the underlying causes and physical natures of thermodynamic phase transitions therein, and in each case, determine optimal approximate descriptions.

Arbitrary gauge quantisation of the matter-radiation system

Throughout this section we will frequently use the Helmholtz decomposition of a vector field ${\bf V}$ into transverse and longitudinal parts ${\bf V}_{\rm T}$ and ${\bf V}_{\rm L}$ such that for all ${\bf x}$

[TABLE]

We assume that all vector fields vanish at the boundaries $|{\bf x}|\to\infty$, which allows free use of integration by parts such as

[TABLE]

Recalling that $\nabla\times\nabla f({\bf x})={\bf 0}$ for any $f$ and for all ${\bf x}$, we have that for any longitudinal field ${\bf U}_{\rm L}$ there exists an $f$ such that ${\bf U}_{\rm L}=\nabla f$. It follows from Eq. (18) that

[TABLE]

for any vector fields ${\bf V}$ and ${\bf U}$. These formulae will be frequently used in what follows.

We consider $N$ hydrogen-like non-relativistic atoms comprised of positive charges $e_{\rm p}=e$ with mass $m_{\rm p}$ at positions ${\bf r}_{\rm p\mu}$, each of which is paired with a negative charge $e_{\rm e}=-e$ with mass $m_{\rm e}$ at ${\bf r}_{\rm e\mu}$. The four-current has components $(j^{a})=(\rho,{\bf J})$ with

[TABLE]

such that $\partial_{a}j^{a}={\dot{\rho}}({\bf x})+\nabla\cdot{\bf J}=0$. Since the system is globally neutral we can define the polarisation field ${\bf P}$ by the equation $-\nabla\cdot{\bf P}=\rho$, which can be solved to give

[TABLE]

where

[TABLE]

defines the green’s function ${\bf g}$ for the divergence operator. Since $\nabla\cdot{\bf g}({\bf x},{\bf x^{\prime}})\equiv\nabla\cdot{\bf g}_{\rm L}({\bf x},{\bf x}^{\prime})$, Eq. (22) only fixes ${\bf g}_{\rm L}={\bf g}-{\bf g}_{\rm T}$ uniquely as

[TABLE]

Any field ${\bf g}_{\rm T}$ with $\nabla\cdot{\bf g}_{\rm T}({\bf x},{\bf x}^{\prime})=0$, can be added to ${\bf g}_{\rm L}$ in Eq. (24) to obtain a ${\bf g}$ that satisfies Eq. (22). It follows that ${\bf P}_{\rm L}$ is fixed uniquely by Eqs. (22) and (24) while ${\bf P}_{\rm T}$ is arbitrary, being determined by Eq. (22) and ${\bf g}_{\rm T}$. To avoid any confusion we note that here we are using the notation ${\bf P}_{\rm T}$ to refer to the arbitrary transverse part of ${\bf P}$ given in Eq. (22). In the main text we reserve the notation ${\bf P}_{\rm T}$ for specifically the multipolar transverse polarisation, which is just one possible example of the transverse part of ${\bf P}$.

Assuming for generality an external potential $V_{\rm ext}$ acting on the charges the standard Lagrangian describing the system of all charges coupled to the Maxwell field is given by

[TABLE]

where $(A^{a})=(A_{0},{\bf A})$ are the components of the electromagnetic four-potential and $F_{ab}=\partial_{a}A_{b}-\partial_{b}A_{a}$. The field tensor components $F_{ab}$ are invariant under a gauge transformation $A_{a}\to A_{a}-\partial_{a}\chi$ where $\chi$ is an arbitrary function. We encode this gauge-freedom into the arbitrary function ${\bf g}_{\rm T}$ by defining

[TABLE]

and by subsequently defining the arbitrary potentials

[TABLE]

where ${\bf A}_{\rm T}$ is the gauge-invariant transverse vector potential and

[TABLE]

The choice ${\bf g}_{\rm T}={\bf 0}$ defines the Coulomb gauge wherein ${\bf A}={\bf A}_{\rm T}$ and $A_{0}=\phi_{\rm Coul}$. A different gauge can by specified by choosing a different ${\bf g}_{\rm T}$.

The Lagrangian $L$ in Eq. (25) is not gauge-invariant. We therefore define the equivalent, but gauge-invariant Lagrangian

[TABLE]

Using Eqs. (27) and (28) in conjunction with Eq. (26) the Lagrangian $L^{\prime}$ can be written

[TABLE]

where

[TABLE]

in which ${\bf E}_{\rm T}=-{\dot{\bf A}}_{\rm T}$ and ${\bf B}=\nabla\times{\bf A}_{\rm T}$. The remaining total time derivative in Eq. (31) depends on ${\bf P}_{\rm T}$ which according to Eq. (22) is uniquely determined through a choice of gauge ${\bf g}_{\rm T}$. It is straightforward to show that

[TABLE]

where now according to Eq. (26) it is the arbitrary function $\chi$ that is determined by the gauge ${\bf g}_{\rm T}$. Note that ${\bf P}_{\rm T}$ in Eq. (31) is a completely arbitrary transverse field and need not coincide with the usual multipolar transverse polarisation field. In writing Eq. (32) we have used the Gauss law $\nabla\cdot{\bf E}_{\rm L}=-\nabla^{2}\phi_{\rm Coul}=\rho$ and integration by parts to separate the transverse and electrostatic parts of the electromagnetic Lagrangian as

[TABLE]

Combining the electrostatic part of Eq (34) with the term $-\int d^{3}x\rho\phi_{\rm Coul}$ coming from the component $-\int d^{3}xA_{0}\rho$ of the interaction Lagrangian we obtain the final electrostatic interaction term $-\int d^{3}x\rho\phi_{\rm Coul}/2$ that appears in Eq. (32).

We now introduce relative and centre-of-mass coordinates for each atom $\mu$, which are defined by

[TABLE]

along with total and reduced masses defined by

[TABLE]

We partition the electrostatic interaction term into intra-atomic and inter-atomic contributions as

[TABLE]

where

[TABLE]

The first term on the right-hand-side in Eq. (40) gives the divergent Coulomb self-energies of the charges at ${\bf r}_{\rm e\mu}$ and ${\bf r}_{\rm p\mu}$. The second term gives the potential energy that binds the charge $-e$ at ${\bf r}_{\rm e\mu}$ to its nucleus $+e$ at ${\bf r}_{\rm p\mu}$. The second term on the right-hand-side of Eq. (39) gives the inter-atomic Coulomb interactions.

We now restrict our attention to gauges in which the arbitrary transverse polarisation takes the form of a weighted multipolar transverse polarisation with weight $\alpha$;

[TABLE]

where $\alpha$ is real and dimensionless and where ${\bf P}_{{\rm T}\mu}$ denotes the transverse polarisation associated with the $\mu$’th atom expressed in terms of the centre-of-mass position ${\bf R}_{\mu}$. With this restriction the gauge is now completely determined by selecting a value of $\alpha$. The derivation of the Dicke-model also requires the electric-dipole approximation (EDA), which for simplicity we implement at this stage rather than later on. In the EDA the charge density within the inter-atomic Coulomb interaction is approximated by the first non-zero term in the multipole expansion of the single atom density $\rho_{\mu}$ about the atomic centre-of-mass at ${\bf R}_{\mu}$, viz.,

[TABLE]

where ${\bf d}_{\mu}=-e{\bf r}_{\mu}$ is the dipole moment of the $\mu$’th dipole. We thereby obtain

[TABLE]

Similarly to Eq. (Arbitrary gauge quantisation of the matter-radiation system) the multipole expansion of the current yields to leading order

[TABLE]

The current-dependent interaction component of $L_{0}$ in Eq. (32) therefore becomes in the EDA

[TABLE]

The multipole expansion of the $\alpha$-dependent transverse polarisation field in Eq. (41) yields to leading order

[TABLE]

At this stage we neglect the nuclear motions ${\dot{\bf R}}_{\mu}$, such that when we arrive at the quantum theory ${\bf R}_{\mu}$ will simply denote the fixed classical position of the $\mu$’th dipole. The polarisation-dependent interaction component of $L^{\prime}$ in Eq. (31) therefore becomes

[TABLE]

Altogether Eqs. (43), (45) and (Arbitrary gauge quantisation of the matter-radiation system) yield the Lagrangian $L^{\prime}$ within the EDA as

[TABLE]

where we have absorbed the external potential $V_{\rm ext}$ into the definition of the total intra-atomic potential as

[TABLE]

The electrostatic energies $V_{\mu}$ and $V_{\rm dip}$ are defined in Eqs. (40) and (43) respectively, while the transverse electric and magnetic fields are given by ${\bf E}_{\rm T}=-{\dot{\bf A}}_{\rm T}$ and ${\bf B}=\nabla\times{\bf A}_{\rm T}$ respectively. Thus, the Lagrangian in Eq. (48) is fully specified in terms of the dynamical variable set $\{{\bf r}_{\mu},~{}{\dot{\bf r}}_{\mu},{\bf A}_{\rm T},~{}{\dot{\bf A}}_{\rm T}\}$ together with the fixed dipolar positions ${\bf R}_{\mu}$.

It is now possible to switch to the canonical formalism by defining the canonical momenta

[TABLE]

and to then quantise the theory by assuming the canonical commutation relations

[TABLE]

The centre-of-mass variable ${\bf R}_{\mu}$ is a classical position of the $\mu$’th dipole. The $\alpha$-dependent canonical momenta are found to be

[TABLE]

For any two values $\alpha$ and $\alpha^{\prime}$ of the gauge parameter the canonical operators are related by the unitary gauge-fixing transformation $R_{\alpha\alpha^{\prime}}$ as

[TABLE]

where

[TABLE]

The Hamiltonian is defined by

[TABLE]

Through substitution of Eqs. (54) and (55) into Eq. (59) the Hamiltonian written in terms of the manifestly gauge-invariant operators $\{{\bf r}_{\mu},{\dot{\bf r}}_{\mu},{\bf A}_{\rm T},{\dot{\bf A}}_{\rm T}\}$ is found to coincide with the total energy expressed as the sum of material and transverse-electromagnetic energies;

[TABLE]

While this expression is clearly $\alpha$-independent (gauge-invariant), when expressed in terms of canonical operators the Hamiltonian has an $\alpha$-dependent functional form given by

[TABLE]

In writing Eq. (Arbitrary gauge quantisation of the matter-radiation system) we have used

[TABLE]

which follows from Eqs. (43) and (Arbitrary gauge quantisation of the matter-radiation system), together with the property $\delta^{\rm T}_{ij}({\bf x})=-\delta^{\rm L}_{ij}({\bf x})$ for ${\bf x}\neq{\bf 0}$. The Hamiltonian reduces to the Coulomb-gauge result if we choose $\alpha=0$. In this case direct electrostatic interactions are fully explicit in the form of the dipole-dipole interaction $(1-\alpha^{2})V_{\rm dip}=V_{\rm dip}$. The field degrees of freedom are defined in terms of the transverse vector potential and its velocity, the transverse electric field; ${\bf\Pi}_{{\rm T}0}={\dot{\bf A}}_{\rm T}=-{\bf E}_{\rm T}$. Another common choice of gauge is the multipolar gauge obtained by choosing $\alpha=1$. In this gauge electrostatic interactions are eliminated; $(1-\alpha^{2})V_{\rm dip}=0$, while the field degrees of freedom are defined in terms of the transverse vector potential and the retarded transverse displacement field; ${\bf\Pi}_{{\rm T}1}=-{\bf D}_{\rm T}=-{\bf E}_{\rm T}-{\bf P}_{{\rm T}1}$. Outside of the atoms the transverse displacement field coincides with the total electric field, which in the EDA means that ${\bf D}_{\rm T}({\bf x})={\bf E}({\bf x})$ for ${\bf x}\neq{\bf R}_{\mu}$. More generally, in the $\alpha$-gauge the Hamiltonian has a hybrid form. Coulomb-gauge matter-transverse field interaction terms are weighted by $(1-\alpha)$ while multipolar-gauge matter-transverse field interaction terms are weighted by $\alpha$. Electrostatic interaction terms are weighted by $1-\alpha^{2}$. In the case $N=2$ the Hamiltonian in Eq. (Arbitrary gauge quantisation of the matter-radiation system) coincides with that given in Ref. Stokes and Nazir (2018).

The above expressions are applicable for general field operators ${\bf A}_{\rm T}$ and ${\bf\Pi}_{{\rm T}\alpha}$. We now define the operator

[TABLE]

where ${\tilde{A}}_{{\rm T},\lambda}({\bf k})={\bm{\varepsilon}}_{\lambda}({\bf k})\cdot{\tilde{\bf A}}_{{\rm T}}({\bf k})$ and ${\tilde{\Pi}}_{{\rm T}\alpha,\lambda}({\bf k})={\bm{\varepsilon}}_{\lambda}({\bf k})\cdot{\tilde{\bf\Pi}}_{{\rm T}\alpha}({\bf k})$. Here tildes denote the Fourier transform and ${\bm{\varepsilon}}_{\lambda}({\bf k}),~{}\lambda=1,\,2$ are mutually orthogonal unit vectors both orthogonal to ${\bf k}$. From the transverse canonical commutation relation

[TABLE]

it follows that

[TABLE]

The operators $a_{\alpha,\lambda}({\bf k})$ and $a_{\alpha,\lambda}^{\dagger}({\bf k})$ are recognisable as annihilation and creation operators for a photon with momentum ${\bf k}$ and polarisation $\lambda$. In terms of these operators the canonical fields support the Fourier representations

[TABLE]

where $\omega=|{\bf k}|$ and $g:={1/\sqrt{2\omega(2\uppi)^{3}}}$.

If we assume an implicit cavity with volume $v$ that satisfies periodic boundary conditions, the continuous label ${\bf k}$ becomes discrete. The pair ${\bf k}\lambda$ then labels a radiation mode and the operators $a_{\alpha,\lambda}({\bf k})$ are labelled with discrete index $\cdot_{\bf k}$ as $a_{\alpha,{\bf k}\lambda}$, which satisfy

[TABLE]

As a less realistic, but simpler model for the cavity we may restrict our attention to a single fixed mode ${\bf k}\lambda$ and ignore all modes ${\bf k}^{\prime}\lambda^{\prime}\neq{\bf k}\lambda$. In this case the field operators become

[TABLE]

where $g=1/\sqrt{2\omega v}$, $\omega=|{\bf k}|$, ${\bm{\varepsilon}}\equiv{\bm{\varepsilon}}_{{\bf k}\lambda}$. and $a_{\alpha}\equiv a_{\alpha,{\bf k}\lambda}$ with $[a_{\alpha},a_{\alpha}^{\dagger}]=1$. Eqs. (68) and (69) imply that the cavity canonical operators now satisfy the commutation relation

[TABLE]

To preserve Eq. (64) we must discretise the ${\bf k}$-space representation of the transverse delta-function and subsequently perform the single-mode approximation as

[TABLE]

With this the $\alpha$-gauge transverse material polarisation becomes

[TABLE]

Similarly, $V_{\rm dip}$ becomes

[TABLE]

Altogether, within the single-mode theory the Hamiltonian in Eq. (Arbitrary gauge quantisation of the matter-radiation system) becomes

[TABLE]

where ${\bf A}_{\rm T}({\bf R}_{\mu})$, ${\bf\Pi}_{{\rm T}\alpha}({\bf R}_{\mu})$ and $V_{\rm dip}$ are given by Eqs. (68), (69) and (73) respectively.

Within the single-mode restriction the Heisenberg equation with the Hamiltonian in Eq. (74) yields

[TABLE]

as required according to Eq. (55). Because the single-mode restriction has been imposed on both the mode operators and the material transverse polarisation the fundamental kinematic relations of the Hamiltonian theory, namely Eqs. (54) and (55), are preserved. We therefore obtain a self-consistent theory describing $N$ dipoles and a single-mode of radiation. We note that as in the multi-mode theory of Eq. (Arbitrary gauge quantisation of the matter-radiation system) electrostatic inter-dipole interactions are explicit in Eq. (74) in all gauges other than the multipolar gauge $\alpha=1$.

Finally we consider the limit of closely spaced dipoles around the origin ${\bf 0}$ such that ${\bf R}_{\mu}\approx{\bf 0}$. In this case, according to Eqs. (72) and (73) we obtain

[TABLE]

The canonical fields ${\bf A}_{\rm T}({\bf R}_{\mu})$ and ${\bf\Pi}({\bf R}_{\mu})$ within the Hamiltonian are replaced by ${\bf A}:={\bf A}_{\rm T}({\bf 0})$ and ${\bf\Pi}_{\alpha}:={\bf\Pi}_{\rm T\alpha}({\bf 0})$ respectively, where for notational convenience we have dropped the transversality subscript $\cdot_{\rm T}$. We also adopt this convention in the main text.

Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians

Here we diagonalise a coupled oscillator Hamiltonian with the generic structure that we will repeatedly encounter. We define arbitrary oscillator operators $\{y,\,y^{\dagger},\,z,\,z^{\dagger}\}$ with $[y,y^{\dagger}]=1=[z,z^{\dagger}]$ and where all other commutators between elements of $\{y,\,y^{\dagger},\,z,\,z^{\dagger}\}$ vanish. The generic Hamiltonian we wish to diagonalise is

[TABLE]

where $C$ is a constant. As an example, the Hamiltonian in Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians), which describes the thermodynamic limit of the $\alpha$-gauge Dicke model in the normal phase, has the form of $H$ above.

To diagonalise $h$ we introduce Hermitian quadratures $q_{\mu}=(\mu^{\dagger}+\mu)/\sqrt{2}$ and $p_{\mu}=i(\mu^{\dagger}-\mu)/\sqrt{2}$ where $\mu=y,\,z$. Subsequently we define the tuple of quadratures ${\bf r}=(q_{y},q_{z},p_{y},p_{z})$, which is such that $[r_{j},r_{k}]=i\Omega_{jk}$ where

[TABLE]

is a matrix representation of the standard symplectic form on ${\mathbb{R}}^{2\times 2}$. The Hamiltonian in Eq. (78) can now be written

[TABLE]

where $\cdot^{\rm T}$ denotes transposition and

[TABLE]

is assumed to be positive-definite. By Williamson’s theorem Williamson (1936) there exists a symplectic matrix $\Lambda$ such that

[TABLE]

where $D$ is diagonal. Denoting the elements of $D$ by $\nu_{j},\,j=1,2$, the quantity $\pm i\nu_{j}$ is an eigenvalue of $\Omega M$. We therefore make use of the canonically transformed quadratures ${\bf r}^{\prime}=(q^{\prime}_{y},q^{\prime}_{z},p^{\prime}_{y},p^{\prime}_{z})=\Lambda^{-1}{\bf r}$, which because $\Lambda$ is symplectic also satisfy $[r^{\prime}_{j},r^{\prime}_{k}]=i\Omega_{jk}$. The Hamiltonian $h$ can now be written in terms of upper and lower polaritons as

[TABLE]

where $y^{\prime},\,z^{\prime}$ are bosonic operators defined in terms of the transformed quadratures ${\bf r}^{\prime}$. They satisfy $[y^{\prime},y^{\prime\dagger}]=1=[z^{\prime},z^{\prime\dagger}]$ while all other commutators between elements of $\{y^{\prime},\,y^{\prime\dagger},\,z^{\prime},\,z^{\prime\dagger}\}$ vanish. The energies $E_{\pm}$ are given by the elements of $D$ as $E_{+}^{2}=(2\nu_{1})^{2}$ and $E_{-}^{2}=(2\nu_{2})^{2}$. The $\nu_{j}$ are found from the matrix $\Omega M$, which is found from Eq. (87). Explicitly, the polariton energies $E_{\pm}$ are given by

[TABLE]

The $\alpha$-gauge Dicke model Hamiltonian is

[TABLE]

where to obtain this expression the transverse vector potential and its conjugate momentum at the dipolar positions ${\bf R}_{\mu}$ have been approximated by their values at the origin as described in the previous section, and they are given by

[TABLE]

The bosonic operators $c_{\alpha}$ include the contribution of the ${\bf A}^{2}$-term of the Hamiltonian implicitly. This is why the frequency appearing in the corresponding mode expansions above is the renormalised frequency $\omega_{\alpha}$ and why there is no explicit ${\bf A}^{2}$-term within the Hamiltonian in Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians). The $c_{\alpha}$ are related to the unrenormalised bosonic operators $\alpha_{\alpha}$ by a local Bogoliubov tranformation within the cavity Hilbert space.

Substituting

[TABLE]

into Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians) gives

[TABLE]

All terms that depend on the square-root functions of the mode operators $b_{\alpha},\,b_{\alpha}^{\dagger}$ have coefficients that remain finite in the thermodynamic limit. Therefore, expanding the square-roots as

[TABLE]

and ignoring terms which vanish in the thermodynamic limit ($N\to\infty$) constitutes making the replacement

[TABLE]

in Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians), which yields

[TABLE]

We now combine the terms that depend only on $b_{\alpha},~{}b_{\alpha}^{\dagger}$ via a Bogoliubov transformation such that

[TABLE]

where $[l_{\alpha},l_{\alpha}^{\dagger}]=1$ and

[TABLE]

The quadrature operators defined in terms of the $b_{\alpha}$ are given in terms of the $l_{\alpha}$ by

[TABLE]

Substituting these expressions into Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians) gives

[TABLE]

where we have combined all constant terms into

[TABLE]

Note that if we choose the material potential such that $\epsilon_{0}=0$, i.e., if the material energy zero-point is zero, then $C_{\alpha}$ is $N$-independent and remains finite in the limit $N\to\infty$. The Hamiltonian in Eq. (105) can be diagonalised using the method presented at the start of this section, which leads to the final result given in the main text.

To find the abnormal-phase Hamiltonian, the relevant starting Hamiltonian $H^{\alpha,2}$ is again that of Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians), together with the Holstein-Primakoff representation

[TABLE]

We then introduce the displaced operator $f_{\alpha}=b_{\alpha}-\sqrt{\beta_{\alpha}}$ where $\beta_{\alpha}$ is assumed to be of order $N$, such that the material part $H_{m}^{\alpha,2}$ of $H^{\alpha,2}$ can be written

[TABLE]

where for convenience we have defined

[TABLE]

Expanding $\sqrt{\xi_{\alpha}}$ in Eq. (Diagonalisation of generic bilinear coupled-oscillator Hamiltonian, normal-phase, and abnormal-phase Hamiltonians) and neglecting terms which vanish in the thermodynamic limit one obtains after lengthy manipulations

[TABLE]

Next we replace the radiation mode operators with displaced operators such that $c_{\alpha}=i(k_{\alpha}+\sqrt{\gamma_{\alpha}})$. The total Hamiltonian therefore reads

[TABLE]

We now collect all terms that are linear in the $f_{\alpha}$ or in the $k_{\alpha}$ and choose $\beta_{\alpha}$ and $\gamma_{\alpha}$ such that these terms vanish. The trivial case in which $\beta_{\alpha}=0=\gamma_{\alpha}$ yields the normal phase Hamiltonian. We will see that the non-trivial solutions

[TABLE]

where

[TABLE]

yield a Hamiltonian describing the abnormal phase.

Expanding $\sqrt{\xi_{\alpha}}$ and neglecting terms which vanish in the thermodynamic limit one now obtains after lengthy manipulations

[TABLE]

We can remove the term quadratic in $f_{\alpha}^{\dagger}+f_{\alpha}$ by defining new material mode operators $f_{\alpha}^{\prime},\,f_{\alpha}^{\prime\dagger}$ such that

[TABLE]

where

[TABLE]

Letting $\omega_{m}^{\prime}=\omega_{m}(1+\tau)/(2\tau)$ and using the relations

[TABLE]

the Hamiltonian can be written

[TABLE]

where $c^{\prime}_{\alpha}=ik_{\alpha}$ while

[TABLE]

and

[TABLE]

The Hamiltonian in Eq. (119) has the form in Eq. (78) and can therefore be diagonalised by the method presented at the start of this section, leading to the final result denoted $H_{2,{\rm th}}^{\alpha,{\rm a}}$, and given in the main text.

Polariton Energies

Here we plot the polariton energies $E_{\alpha}^{\pm}$ in Fig. 2 using the example of double-well dipoles as considered in the main text. The thermodynamic limit of the Dicke-model is seen to be gauge-dependent, i.e., the $E_{\alpha}^{\pm}$ are $\alpha$-dependent. This occurs due to the use of material level truncation. Despite this, it is clearly seen that the occurrence of a unique phase transition is obtained as a gauge-invariant prediction. Furthermore, we are able to determine its gauge-invariant manifestation, which is plotted in the following section.

Calculation of radiative canonical operator averages

We begin with the cavity canonical operators $A={\bm{\varepsilon}}\cdot{\bf A}$ and $\Pi_{\alpha}={\bm{\varepsilon}}\cdot{\bf\Pi}_{\alpha}$. The ground state average of $A$ is trivially zero. Similarly the ground state average of $\Pi_{\alpha}$ in the normal phase is zero. The ground state average of $\Pi_{\alpha}$ in the abnormal phase can be calculated using the Dicke-model of any gauge $\alpha^{\prime}$. We begin with the expression

[TABLE]

which is the $\alpha^{\prime}$-gauge’s two-level approximation of $\Pi_{\alpha}$. Using the Holstein-Primakoff representation and then defining the displaced operators $f_{\alpha},~{}c_{\alpha}^{\prime}$ by

[TABLE]

one obtains

[TABLE]

where $\beta,~{}\gamma_{\alpha}$ and $\xi_{\alpha}$ are defined in Eqs. (112), (113), and (109) respectively. Expanding $\sqrt{\xi_{\alpha^{\prime}}}$ and retaining terms which do not vanish in the thermodynamic limit yields

[TABLE]

as given in the main text. We plot $\Pi_{\alpha,\rm th}^{\alpha^{\prime},2,{\rm a}}$ in Fig. 3 for different $\alpha$ again using the example of double-well dipoles. As a special case, this includes the gauge-invariant macroscopic manifestation of the abnormal phase as quantified by $P_{{\rm T}\alpha,\rm th}^{\alpha^{\prime},2,{\rm a}}$

Further numerical results

Here we consider a less anharmonic single-dipole double-well potential, which has $(\epsilon_{2}-\epsilon_{0})/\omega_{m}\approx 3.2$. This results from choosing $\beta=1.5$ rather than $\beta=3.3$ as was chosen in the main text. In this case single-dipole two-level models are able to remain accurate in predicting the low energy properties of the system, but the multipolar gauge no longer provides the optimal two-level model. The optimal gauge for the two-level truncation is shifted towards the Coulomb gauge, such that the Jaynes-Cummings gauge two-level model is close to optimal Stokes and Nazir (2019a). This is shown in Fig. 4, which compares $G$ (Fig. 4a), $E-G$ (Fig. 4b), and $d^{2}G/d\eta^{2}$ (Fig. 4c) each obtained from the Coulomb-gauge, Jaynes-Cummings gauge, multipolar-gauge two-level models, and the exact (non-truncated) theory. For $N>1$, two-level models become less accurate in predicting even low energy properties when the coupling is sufficiently strong, as shown for the case $N=2$ in Fig. 5 and for the case $N=3$ in Fig. 6. The accuracy of multipolar-gauge two-level truncation appears to improve relative to the other gauges, but the ground energy is not well represented by any two-level model for sufficiently strong coupling.