Resolution of Gauge Ambiguities in Molecular Cavity Quantum Electrodynamics

Coupling molecules to the quantized radiation field inside an optical cavity creates a set of new photon-matter hybrid excitations, so-called polaritons [1, 2, 3, 4]. The rich dynamic interplay among these electronic, photonic, and nuclear degrees of freedom (DOF) has enabled a new paradigm for achieving unique chemical reactivities [5, 6, 7, 8, 9]. The non-relativistic quantum electrodynamics (QED) Hamiltonian that describes such quantum light-matter interactions should obey the gauge principle, i.e., giving rise to the same physical results (physical observables) upon a gauge transformation [10, 11]. While the QED Hamiltonian under both the Coulomb gauge and the dipole gauge (length gauge) indeed obeys this principle, these Hamiltonians under a finite electronic state truncation (the few-level approximation) are known to give different results for physical observable [12, 13, 14, 15, 16, 17, 18, 19], which is commonly referred to as the gauge ambiguity [14, 16, 17, 18]. In this letter, we provide a conjecture for the cause of this ambiguity, and demonstrate that it is possible to resolve it, hence providing identical physical results from either the Coulomb or the dipole gauge Hamiltonian upon the same level of the electronic state truncation.

We begin by defining the matter Hamiltonian and the corresponding total dipole operator as follows

$$\hat{H}_{\mathrm{M}}=\hat{\bf T}+\hat{V}({\hat{\bf x}})=\sum_{j}\frac{1}{2m_{j}}\hat{\bf p}_{j}^{2}+\hat{V}({\hat{\bf x}});~{}\hat{\bm{\mu}}=\sum_{j}{z}_{j}\hat{\bf x}_{j},$$

(1)

where $j$ is the index of the $j_{\mathrm{th}}$ charged particle (including all electrons and nuclei), with the corresponding mass $m_{j}$ and charge ${z}_{j}$. In addition, $\hat{\bf x}\equiv\{\hat{\bf x}_{j}\}=\{\hat{\bf R},\hat{\bf r}\}$ with $\hat{\bf R}$ and $\hat{\bf r}$ represent the nuclear and electronic coordinates, respectively, $\hat{\bf p}\equiv\{\hat{\bf p}_{\bf R},\hat{\bf p}_{\bf r}\}\equiv\{\hat{\bf p}_{j}\}$ is the mechanical momentum operator as well as the canonical momentum operator, such that $\hat{\bf p}_{j}=-i\hbar{\bm{\nabla}}_{j}$. Further, $\hat{\bf T}=\hat{\bf T}_{\bf R}+\hat{\bf T}_{\bf r}$ is the kinetic energy operator, where $\hat{\bf T}_{\bf R}$ and $\hat{\bf T}_{\bf r}$ represent the kinetic energy operator for nuclei and for electrons, respectively, and $\hat{V}({\hat{\bf x}})$ is the potential operator that describes the Coulombic interactions among electrons and nuclei.

The cavity photon field Hamiltonian under the single mode assumption is expressed as

$$\hat{H}_{\mathrm{ph}}=\hbar\omega_{\mathrm{c}}\big{(}\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\big{)}=\frac{1}{2}\big{(}\hat{p}_{\mathrm{c}}^{2}+\omega_{\mathrm{c}}^{2}\hat{q}_{\mathrm{c}}^{2}\big{)},$$

(2)

where $\omega_{\mathrm{c}}$ is the frequency of the mode in the cavity, $\hat{a}^{\dagger}$ and $\hat{a}$ are the photonic creation and annihilation operators, $\hat{q}_{\mathrm{c}}=\sqrt{\hbar/2\omega_{\mathrm{c}}}(\hat{a}^{\dagger}+\hat{a})$ and $\hat{p}_{\mathrm{c}}=i\sqrt{\hbar\omega_{\mathrm{c}}/2}(\hat{a}^{\dagger}-\hat{a})$ are the photonic coordinate and momentum operators, respectively. Choosing the Coulomb gauge, ${\bm{\nabla}}\cdot\hat{\bf A}=0$, the vector potential becomes purely transverse $\hat{\bf A}=\hat{{\bf A}}_{\perp}$. Under the long-wavelength approximation,

$\displaystyle\hat{\bf A}={\bf A}_{0}\big{(}\hat{a}+\hat{a}^{\dagger}\big{)}={\bf A}_{0}\sqrt{{2\omega_{\mathrm{c}}}/{\hbar}}~{}\hat{q}_{\mathrm{c}}$

(3)

where ${\bf A}_{0}=\sqrt{\hbar/2\omega_{\mathrm{c}}\varepsilon_{0}\mathcal{V}}~{}\hat{\bf e}$, with $\mathcal{V}$ as the quantization volume inside the cavity, $\varepsilon_{0}$ as the permittivity, and $\hat{\bf e}$ is the unit vector of the field polarization.

The minimal coupling QED Hamiltonian in the Coulomb gauge (the “$\mathrm{p\cdot A}$” form) is expressed as

$$\hat{H}_{\mathrm{C}}=\sum_{j}\frac{1}{2m_{j}}(\hat{\bf p}_{j}-{z}_{j}\hat{\bf A})^{2}+\hat{V}(\hat{\bf x})+\hat{H}_{\mathrm{ph}},$$

(4)

where $\hat{\bf p}_{j}=-i\hbar{\bm{\nabla}}_{j}$ is the canonical momentum operator. Upon a gauge transformation $\hat{U}_{\chi}=\exp[\frac{i}{\hbar}\sum_{j}z_{j}{\bm{\chi}}(\hat{\bf x}_{j})]$, $\hat{H}_{\mathrm{\chi}}=\hat{U}_{\chi}\hat{H}_{\mathrm{C}}\hat{U}_{\chi}^{\dagger}$ remains gauge-invariant as

$$\hat{H}_{\mathrm{\chi}}=\sum_{j}\frac{1}{2m_{j}}(\hat{\bf p}_{j}-{z}_{j}\hat{\bf A}_{\chi}(\hat{\bf x}_{j}))^{2}+\hat{V}(\hat{\bf x})+\hat{H}^{\chi}_{\mathrm{ph}},$$

(5)

where $\hat{H}^{\chi}_{\mathrm{ph}}=\hat{U}_{\chi}\hat{H}_{\mathrm{ph}}\hat{U}_{\chi}^{\dagger}$, $\hat{U}_{\chi}\hat{V}(\hat{\bf x})\hat{U}^{\dagger}_{\chi}=\hat{V}(\hat{\bf x})$ because $\hat{V}$ is a local potential operator for the matter, that is, only a function of $\hat{\bf x}$ and $\hat{\bf p}$-independent, and $\hat{\bf A}_{\chi}({\bf x}_{j})=\hat{\bf A}+{\bm{\nabla}}_{j}\chi(\hat{\bf x}_{j})$ is gauge transformed vector potential that provides the same physical field, because ${\bm{\nabla}}_{j}\times{\bm{\nabla}}_{j}\chi(\hat{\bf x}_{j})=0$. To prove Eqn. 5, we have used $e^{\hat{Y}}\hat{O}(\hat{X})e^{-\hat{Y}}=\hat{O}(e^{\hat{Y}}\hat{X}e^{-\hat{Y}})$ for a unitary operator, as well as $\hat{U}_{\chi}\hat{\bf p}_{j}\hat{U}_{\chi}^{\dagger}=\hat{\bf p}_{j}-z_{j}{\bm{\nabla}}_{j}{\bm{\chi}(\hat{\bf x}_{j})}$.

We further introduce the Power-Zienau-Woolley (PZW) gauge transformation operator [20, 21, 10] as

$$\hat{U}=\exp[-\frac{i}{\hbar}\hat{\bm{\mu}}\cdot\hat{\bf A}\big{]}=\exp[-\frac{i}{\hbar}\hat{\bm{\mu}}\cdot{\bf A}_{0}\big{(}\hat{a}+\hat{a}^{\dagger}\big{)}\big{]},$$

(6)

or $\hat{U}=\exp[-\frac{i}{\hbar}\sqrt{2\omega_{\mathrm{c}}/\hbar}\hat{\bm{\mu}}{\bf A}_{0}\hat{q}_{\mathrm{c}}\big{]}=\exp[-\frac{i}{\hbar}(\sum_{j}{z}_{j}\hat{\bf A}{\bf x}_{j})\big{]}$. Recall that a momentum boost operator $\hat{U}_{\mathrm{p}}=e^{-\frac{i}{\hbar}p_{0}\hat{q}}$ displaces $\hat{p}$ by the amount of $p_{0}$, such that $\hat{U}_{\mathrm{p}}\hat{O}(\hat{p})\hat{U}_{\mathrm{p}}^{\dagger}=\hat{O}(\hat{p}+p_{0})$. Hence, $\hat{U}$ is a boost operator for both the photonic momentum $\hat{p}_{\mathrm{c}}$ by the amount of $\sqrt{2\omega_{\mathrm{c}}/\hbar}\hat{\bm{\mu}}{\bf A_{0}}$, as well as for the matter momentum $\hat{\bf p}_{j}$ by the amount of ${z}_{j}\hat{\bf A}$. The PZW gauge operator (Eqn. 6) is a special case of $\hat{U}_{\chi}$, such that ${\bm{\chi}}=-\hat{\bf x}_{j}\cdot\hat{\bf A}$. Using $\hat{U}^{\dagger}$ to boost the matter momentum, one can show that

$$\hat{H}_{\mathrm{C}}=\hat{U}^{\dagger}\hat{H}_{\mathrm{M}}\hat{U}+\hat{H}_{\mathrm{ph}},$$

(7)

hence $\hat{H}_{\mathrm{C}}$ can be obtained [18] by a momentum boost with the amount of $-{z}_{j}\hat{\bf A}$ for $\hat{\bf p}_{j}$, then adding $\hat{H}_{\mathrm{ph}}$.

The QED Hamiltonian under the dipole gauge (the “$\mathrm{d\cdot E}$” form [20, 21]) can be obtained by performing the PZW transformation on $\hat{H}_{\mathrm{C}}$ as follows

	$\displaystyle\hat{H}_{\mathrm{D}}$	$\displaystyle=\hat{U}\hat{H}_{\mathrm{C}}\hat{U}^{\dagger}=\hat{U}\hat{U}^{\dagger}\hat{H}_{\mathrm{M}}\hat{U}\hat{U}^{\dagger}+\hat{U}\hat{H}_{\mathrm{ph}}\hat{U}^{\dagger}$		(8)
		$\displaystyle=\hat{H}_{\mathrm{M}}+\hbar\omega_{\mathrm{c}}(\hat{a}^{\dagger}\hat{a}+\frac{1}{2})+i\omega_{\mathrm{c}}\hat{\bm{\mu}}{\bf A}_{0}(\hat{a}^{\dagger}-\hat{a})+\frac{\omega_{\mathrm{c}}}{\hbar}(\hat{\bm{\mu}}{\bf A}_{0})^{2},$

where we have used Eqn. 7 to express $\hat{H}_{\mathrm{C}}$, and the last three terms of the above equation are the results of $\hat{U}\hat{H}_{\mathrm{ph}}\hat{U}^{\dagger}$ (see SI). The above result is also apparent by examining Eqn. 5, where the gauge DOF $-z_{j}{\bm{\nabla}}{\bm{\chi}}(\hat{\bf x}_{j})=z_{j}\hat{\bf A}$ exactly cancels the vector potential, leaving no explicit gauge-dependent quantity in $\hat{H}_{\mathrm{D}}$. Using $\hat{q}_{\mathrm{c}}$ and $\hat{p}_{\mathrm{c}}$, one can instead show that $\hat{H}_{\mathrm{D}}=\hat{H}_{\mathrm{M}}+\frac{1}{2}\omega_{\mathrm{c}}^{2}\hat{q}_{\mathrm{c}}^{2}+\frac{1}{2}(\hat{p}_{\mathrm{c}}+\sqrt{2\omega_{\mathrm{c}}/\hbar}\hat{\bm{\mu}}{\bf A}_{0})^{2}$, because the PZW operator boosts the photonic momentum $\hat{p}_{\mathrm{c}}$ by $\sqrt{2\omega_{\mathrm{c}}/\hbar}\hat{\bm{\mu}}{\bf A}_{0}$.

The gauge invariance is explicitly enforced between $\hat{H}_{\mathrm{C}}$ (Eqn. 4) and $\hat{H}_{\mathrm{D}}$ (Eqn. 8) through the unitary PZW Gauge transformation (Eqn. 6). However, this gauge invariance will explicitly breakdown when a truncation of electronic states is applied to both Hamiltonians [17, 18]. Consider a finite subset of electronic states $\{|\alpha\rangle\}$, where the projection operator $\hat{\mathcal{P}}=\sum_{\alpha}|\alpha\rangle\langle\alpha|$ defines the truncation of the full electronic Hilbert space $\hat{\mathds{1}}_{\mathrm{r}}=\hat{\mathcal{P}}+\hat{\mathcal{Q}}$ to the corresponding subspace $\hat{\mathcal{P}}$. The truncation reduces the size of the Hilbert space from originally $\hat{\mathds{1}}_{\mathrm{r}}\otimes\hat{\mathds{1}}_{\mathrm{R}}\otimes\hat{\mathds{1}}_{\mathrm{ph}}$ to now $\hat{\mathcal{P}}\otimes\hat{\mathds{1}}_{\mathrm{R}}\otimes\hat{\mathds{1}}_{\mathrm{ph}}$, where $\hat{\mathds{1}}_{\mathrm{R}}$ and $\hat{\mathds{1}}_{\mathrm{ph}}$ represent the identity operator of the nuclear and the photonic DOF, respectively. The truncated matter Hamiltonian is

$$\hat{\mathcal{H}}_{\mathrm{M}}=\hat{\mathcal{P}}\hat{H}_{\mathrm{M}}\hat{\mathcal{P}}=\hat{\mathcal{P}}\hat{\bf T}\hat{\mathcal{P}}+\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}.$$

(9)

Throughout this letter, we use calligraphic symbols (such as $\hat{\mathcal{H}}_{\mathrm{M}}$) to indicate operators in the truncated Hilbert space. Truncating the momentum operator and dipole operator as $\hat{\mathcal{P}}\hat{\bf p}_{j}\hat{\mathcal{P}}$ and $\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}$, the QED Hamiltonians under the truncated subspace are commonly defined as

	$\displaystyle\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$	$\displaystyle=\hat{\mathcal{P}}\hat{H}_{\mathrm{C}}\hat{\mathcal{P}}=\hat{\mathcal{P}}\hat{U}^{\dagger}\hat{H}_{\mathrm{M}}\hat{U}\hat{\mathcal{P}}+\hat{H}_{\mathrm{ph}}$		(10)
		$\displaystyle=\hat{\mathcal{H}}_{\mathrm{M}}+\hat{H}_{\mathrm{ph}}+\sum_{j}\big{(}-\frac{{z}_{j}}{m_{j}}\hat{\mathcal{P}}\hat{\bf p}_{j}\hat{\mathcal{P}}\hat{\bm{A}}+\frac{{z}_{j}^{2}\hat{\bm{A}}^{2}}{2m_{j}}\big{)}$
	$\displaystyle\hat{\mathcal{H}}_{\mathrm{D}}$	$\displaystyle=\hat{\mathcal{H}}_{\mathrm{M}}+\hat{H}_{\mathrm{ph}}+i\omega_{\mathrm{c}}\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}{\bf A}_{0}(\hat{a}^{\dagger}-\hat{a})+\frac{\omega_{\mathrm{c}}}{\hbar}(\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}{\bf A}_{0})^{2}.$		(11)

It is a well-known fact that the above two Hamiltonians do not generate the same polariton eigenspectrum [17, 18, 19] under the ultra-strong coupling regime [22], hence explicitly breakdown the gauge invariance, leading to the ambiguities [14, 16] of using which Hamiltonian, $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$ or $\hat{\mathcal{H}}_{\mathrm{D}}$, to compute the physical quantities when applying to a state truncation through $\hat{\mathcal{P}}$. This is attributed [17, 23] to the fact that $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$ usually requires a larger subset of the matter states to converge or generate consistent results of $\hat{\mathcal{H}}_{\mathrm{D}}$ and apparently, under the complete basis limit, they should be gauge invariant. Further, this fundamentally different behavior of $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$ and $\hat{\mathcal{H}}_{\mathrm{D}}$ upon electronic states truncation is also attributed to the fundamental asymmetry of the $\hat{\bf p}$ and $\hat{\bm{\mu}}=\sum_{j}\hat{\bf x}_{j}$ operators, especially when the potential $\hat{V}(\hat{\bf x})$ is highly anharmonic [17].

We conjecture this gauge-ambiguity is caused by the fact that $\hat{\mathcal{P}}\hat{U}^{\dagger}$ and $\hat{U}\hat{\mathcal{P}}$ operators used to construct $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$ (Eqn. 10) do not properly constrain all operators in the subspace $\hat{\mathcal{P}}$, such that some of them are entering into the subspace $\hat{\mathcal{Q}}=\hat{\mathds{1}}_{\mathrm{r}}-\hat{\mathcal{P}}$. Indeed, this is the case for $\hat{U}\hat{\mathcal{P}}=(\hat{\mathcal{P}}+\hat{\mathcal{Q}})\hat{U}\hat{\mathcal{P}}$. This is almost apparent by examining the diamagnetic term ${z}_{j}^{2}\hat{\bm{A}}^{2}/2m_{j}$ in $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$, which is effectively evaluated in the full space $\hat{\mathds{1}}_{\mathrm{r}}$ (based on the Thomas-Reiche-Kuhn sum rule), hence is not properly confined in $\hat{\mathcal{P}}$. This diamagnetic term overestimates what it should be in the subspace [17, 18], and by confining it inside $\hat{\mathcal{P}}$, the results can be significantly improved [18]. Note that using $\hat{\mathcal{P}}\hat{U}\hat{\mathcal{P}}=\sum_{n=1}^{\infty}\frac{1}{n!}(-\frac{i}{\hbar})^{n}\hat{\mathcal{P}}\hat{\bm{\mu}}^{n}\hat{\mathcal{P}}\hat{\bf A}^{n}$ does not resolve this problem either (as oppose to the original claim in Ref. 18), because for all non-linear terms in $\hat{\mathcal{P}}\hat{U}\hat{\mathcal{P}}$ (for $n\geq 2$), $\hat{\mathcal{P}}\hat{\bm{\mu}}^{n}\hat{\mathcal{P}}$ is not properly confined in the $\hat{\mathcal{P}}$ subspace [19]. For example, $\hat{\mathcal{P}}\hat{\bm{\mu}}^{2}\hat{\mathcal{P}}=\hat{\mathcal{P}}\hat{\bm{\mu}}(\hat{\mathcal{P}}+\hat{\mathcal{Q}})\hat{\bm{\mu}}\hat{\mathcal{P}}=(\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}})^{2}+\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{Q}}\hat{\bm{\mu}}\hat{\mathcal{P}}\neq(\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}})^{2}$, where we have used $\hat{\mathcal{P}}\hat{\mathcal{P}}=\hat{\mathcal{P}}$. Nevertheless, using $\hat{\mathcal{P}}\hat{U}\hat{\mathcal{P}}$ still provide slightly better numerical results [19] compared to using $\hat{U}\hat{\mathcal{P}}$, because the former has a fewer terms in $\hat{\mathcal{Q}}$, corroborating our conjecture. Further, this conjecture also provides the reason why $\hat{\mathcal{H}}_{\mathrm{D}}$ needs to be constructed in the form of Eqn. 11, because by doing so, all of the operators are properly confined in the $\hat{\mathcal{P}}$ subspace. Indeed, $\hat{\mathcal{H}}_{\mathrm{D}}$ gives very accurate results compared to the exact solution when a few level truncation is a valid approximation [17]. Note that $\hat{\mathcal{H}}_{\mathrm{D}}\neq\hat{\mathcal{P}}\hat{H}_{\mathrm{D}}\hat{\mathcal{P}}$, because $\hat{\mathcal{P}}\hat{H}_{\mathrm{D}}\hat{\mathcal{P}}$ also has an operator (the Dipole self-energy term) $\hat{\mathcal{P}}\hat{\bm{\mu}}^{2}\hat{\mathcal{P}}$ that is not properly confined inside $\hat{\mathcal{P}}$, and will indeed give much worse numerical results than $\hat{\mathcal{H}}_{\mathrm{D}}$. Numerical examples are provided in SI to support our conjecture.

Based on the above conjecture, the gauge ambiguity will be resolved by defining the following unitary operator

$$\hat{\mathcal{U}}=\exp[-\frac{i}{\hbar}\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}\cdot\hat{\bf A}\big{]}\equiv\exp[-\frac{i}{\hbar}\tilde{\bm{\mu}}(\hat{\bf x},\hat{\bf p})\cdot\hat{\bf A}\big{]},$$

(12)

such that all terms in $\hat{\mathcal{U}}=\sum_{n=0}^{\infty}\frac{1}{n!}(-\frac{i}{\hbar})^{n}(\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}})^{n}\hat{\bf A}^{n}$ are properly confined within the subspace $\hat{\mathcal{P}}$. Note that $\hat{\mathcal{U}}$ is defined in analogous to the PZW gauge operator (Eqn. 6) in the full space, and $\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}\equiv\tilde{\bm{\mu}}(\hat{\bf x},\hat{\bf p})$ in principle is a function of both $\hat{\bf x}$ and $\hat{\bf p}$, due to the finite level projection that ruins the locality of $\hat{\bf x}$ [18, 24]. Further, $\hat{\mathcal{U}}$ is a unitary transformation operator in the $\hat{\mathcal{P}}$ subspace and the identity operator in the subspace of $\hat{\mathds{1}}_{\mathrm{r}}-\hat{\mathcal{P}}$, such that we still have $\hat{\mathcal{U}}\hat{\mathcal{U}}^{\dagger}=\hat{\mathds{1}}_{\mathrm{r}}\otimes\hat{\mathds{1}}_{\mathrm{R}}\otimes\hat{\mathds{1}}_{\mathrm{ph}}=\hat{U}\hat{U}^{\dagger}$.

Using the above operator $\hat{\mathcal{U}}$, one can define the following Coulomb gauge Hamiltonian

$$\hat{\mathcal{H}}_{\mathrm{C}}=\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{M}}\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}},$$

(13)

in analogous to $\hat{H}_{\mathrm{C}}$ in Eqn. 7 in the full space. One can then formally show that $\hat{\mathcal{H}}_{\mathrm{C}}$ (Eqn. 13) and $\hat{\mathcal{H}}_{\mathrm{D}}$ (Eqn. 11) are related through $\hat{\mathcal{U}}$ (Eqn. 12) as follows

$$\hat{\mathcal{U}}\hat{\mathcal{H}}_{\mathrm{C}}\hat{\mathcal{U}}^{\dagger}=\hat{\mathcal{H}}_{\mathrm{M}}+\hat{\mathcal{U}}\hat{H}_{\mathrm{ph}}\hat{\mathcal{U}}^{\dagger}=\hat{\mathcal{H}}_{\mathrm{D}},$$

(14)

as well as $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{D}}\hat{\mathcal{U}}=\hat{\mathcal{H}}_{\mathrm{C}}$. Note that to establish the last equality in Eqn. 14 we have used the fact that $\hat{\mathcal{U}}\hat{H}_{\mathrm{ph}}\hat{\mathcal{U}}^{\dagger}=\hat{\mathcal{U}}(\frac{1}{2}\omega_{\mathrm{c}}^{2}\hat{q}_{\mathrm{c}}^{2}+\frac{1}{2}\hat{p}_{\mathrm{c}}^{2})\hat{\mathcal{U}}^{\dagger}=\frac{1}{2}\omega_{\mathrm{c}}^{2}\hat{q}_{\mathrm{c}}^{2}+\frac{1}{2}(\hat{p}_{\mathrm{c}}+\sqrt{2\omega_{\mathrm{c}}/\hbar}\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}{\bf A}_{0})^{2}$. Thus, we have formally demonstrated that the gauge ambiguities between the Coulomb and dipole gauge Hamiltonians can be resolved for an arbitrary matter-cavity hybrid system, under the same level of electronic state truncation.

To derive the detailed expression of $\hat{\mathcal{H}}_{\mathrm{C}}$, we notice that $\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}\equiv\tilde{\bm{\mu}}(\hat{\bf x},\hat{\bf p})$ is in principle a non-linear function of both $\hat{\bf x}$ and $\hat{\bf p}$, as oppose to $\hat{\bm{\mu}}$ (see Eqn. 1) which is a pure linear function of $\hat{\bf x}$ . Thus, $\hat{\mathcal{U}}^{\dagger}$ no longer just boosts the matter momentum by ${z}_{j}\hat{\bf A}$. Using the Baker-Campbell-Hausdorff (BCH) identity, we have

	$\displaystyle\hat{\mathcal{U}}^{\dagger}\hat{\bf p}_{j}\hat{\mathcal{U}}$	$\displaystyle=\hat{\bf p}_{j}+\frac{i}{\hbar}[\tilde{\bm{\mu}}\hat{\bf A},\hat{\bf p}_{j}]+\frac{1}{2}\big{(}\frac{i}{\hbar}\big{)}^{2}[\tilde{\bm{\mu}}\hat{\bf A},[\tilde{\bm{\mu}}\hat{\bf A},\hat{\bf p}_{j}]]+...$
		$\displaystyle=\hat{\bf p}_{j}-{\bm{\nabla}}_{j}{\tilde{\bm{\mu}}(\hat{\bf x},\hat{\bf p})}\hat{\bf A}+{\bf\tilde{P}}_{j},$		(15)

where ${\bf\tilde{P}}_{j}\equiv\frac{1}{2}\big{(}\frac{i}{\hbar}\big{)}^{2}[\tilde{\bm{\mu}}\hat{\bf A},[\tilde{\bm{\mu}}\hat{\bf A},\hat{\bf p}_{j}]]+...$ is the residual momentum that accounts for terms with more than one commutator in the BCH identity. Hence, under the projection, $\hat{\mathcal{U}}^{\dagger}$ boosts the matter momentum by the amount of $-{\bm{\nabla}}_{j}\tilde{\bm{\mu}}\hat{\bf A}+{\bf\tilde{P}}_{j}$. Similarly, $\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}$ is in principle a non-local potential [15, 18, 24], such that $\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}=\hat{\mathcal{V}}(\hat{\bf x},\hat{\bf p})$. Hence, $\hat{\mathcal{U}}^{\dagger}$ also displaces the matter coordinate as well as boost the matter momentum inside $\hat{\mathcal{V}}(\hat{\bf x},\hat{\bf p})$, whereas $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{V}}\hat{\mathcal{U}}$ can be formally derived through the BCH identity.

Using Eqn. Resolution of Gauge Ambiguities in Molecular Cavity Quantum Electrodynamics, as well as the fact that both $\hat{\mathcal{U}}^{\dagger}$ and $\hat{\mathcal{U}}$ commute with $\hat{\mathcal{P}}$, such that $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{\bf T}\hat{\mathcal{P}}\hat{\mathcal{U}}=\hat{\mathcal{P}}\hat{\mathcal{U}}^{\dagger}\hat{\bf T}\hat{\mathcal{U}}\hat{\mathcal{P}}$, we can derive $\hat{\mathcal{H}}_{\mathrm{C}}$ (in Eqn. 13) as follows

		$\displaystyle\hat{\mathcal{H}}_{\mathrm{C}}=\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{\bf T}\hat{\mathcal{P}}\hat{\mathcal{U}}+\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}}$		(16)
		$\displaystyle=\sum_{j}\frac{1}{2m_{j}}\hat{\mathcal{P}}\big{(}\hat{\bf p}_{j}-{\bm{\nabla}}_{j}\tilde{\bm{\mu}}\hat{\bf A}+{\bf\tilde{P}}_{j})^{2}\hat{\mathcal{P}}+\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{V}}(\hat{\bf x},\hat{\bf p})\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}},$

where the sum $j$ includes all charged particles (electrons and nuclei). Note that $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$ (Eqn. 10) as well as $\hat{H}_{\mathrm{C}}$ (Eqn. 4) only contain the vector potential $\hat{\bf A}$ up to the second order. This is no longer the case for $\hat{\mathcal{H}}_{\mathrm{C}}$ in Eqn. 16. In fact, both the ${\bf\tilde{P}}_{j}$ term and the $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{V}}(\hat{\bf x},\hat{\bf p})\hat{\mathcal{U}}$ term in principle contain infinite orders of $\hat{\bf A}$. It is also self-evident that $\hat{\mathcal{H}}_{\mathrm{C}}$ (Eqn. 16) will return to $\hat{H}_{\mathrm{C}}$ (Eqn. 4) under the complete basis limit, such that $\tilde{\bm{\mu}}\equiv\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}\to\hat{\bm{\mu}}$, thus ${\bm{\nabla}}_{j}\tilde{\bm{\mu}}\to{\bm{\nabla}}_{j}\hat{\bm{\mu}}={z}_{j}$, hence ${\bf\tilde{P}}_{j}\to 0$, as well as $\hat{\mathcal{U}}\to\hat{U}$, hence $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}\hat{\mathcal{U}}\to\hat{U}^{\dagger}\hat{V}(\hat{\bf x})\hat{U}=\hat{V}(\hat{\bf x})$. Unfortunately, $\hat{\mathcal{H}}_{\mathrm{C}}$ no longer remains a gauge-invariant form as indicated by Eqn. 5, except when approaching to the complete basis limit, and of course, $\hat{\mathcal{H}}_{\mathrm{C}}$ is invariant from $\hat{\mathcal{H}}_{\mathrm{D}}$ through the $\hat{\mathcal{U}}$ transformation as indicated in Eqn. 14. Thus, $\hat{\mathcal{H}}_{\mathrm{C}}$ resolves the discrepancy (ambiguity) between the coulomb and dipole gauge by under matter state truncation, with the price that $\hat{\mathcal{H}}_{\mathrm{C}}$ is no longer gauge invariant in general.

To further present an equivalent yet convenient $\hat{\mathcal{H}}_{\mathrm{C}}$ for molecular cavity QED, we use the electronic states associated with the electronic Hamiltonian $\hat{H}_{\mathrm{el}}=\hat{\bf T}_{\bf r}+\hat{V}=\hat{H}_{\mathrm{M}}-\hat{\bf T}_{\bf R}$. The adiabatic electronic states $|\alpha({\bf R})\rangle$ are one type of the most commonly obtained electronic states by solving $\hat{H}_{\mathrm{el}}|\alpha({\bf R})\rangle=(\hat{\bf T}_{\bf r}+\hat{V})|\alpha({\bf R})\rangle=E_{\alpha}({\bf R})|\alpha({\bf R})\rangle$, where $E_{\alpha}({\bf R})$ is the so-called the adiabatic potential energy surface. Using a finite set of $\{|\alpha({\bf R})\rangle\}$ to form $\hat{\mathcal{P}}=\sum_{\alpha}|\alpha({\bf R})\rangle\langle\alpha({\bf R})|$, the projected electronic Hamiltonian is $\hat{\mathcal{H}}_{\mathrm{el}}=\hat{\mathcal{P}}\hat{H}_{\mathrm{el}}\hat{\mathcal{P}}=\sum_{\alpha}E_{\alpha}({\bf R})|\alpha\rangle\langle\alpha|$. Alternatively, diabatic electronic states [25, 26, 27, 28] $\{|\varphi\rangle,|\phi\rangle\}$ can be obtained by the unitary transform [25, 26, 27, 28, 29] from the adiabatic states $|\alpha({\bf R})\rangle$ or by constructions through preserving basic properties [30]. As oppose to adiabatic states, the character of the diabatic states do not depend on ${\bf R}$, such that $\langle\varphi|{\bm{\nabla}}_{\bf R}|\phi\rangle=0$. Using $\hat{\mathcal{P}}=\sum_{\varphi}|\varphi\rangle\langle\varphi|$, the truncated electronic Hamiltonian is $\hat{\mathcal{H}}_{\mathrm{el}}=\hat{\mathcal{P}}\hat{H}_{\mathrm{el}}\hat{\mathcal{P}}=\sum_{\varphi}{\mathcal{V}}_{\varphi\varphi}({\bf R})|\varphi\rangle\langle\varphi|+\sum_{\varphi\neq\phi}{\mathcal{V}}_{\varphi\phi}({\bf R})|\varphi\rangle\langle\phi|$, where ${\mathcal{V}}_{\varphi\phi}({\bf R})=\langle\varphi|\hat{H}_{\mathrm{el}}|\phi\rangle$ is the diabatic matrix element of $\hat{H}_{\mathrm{el}}$ (not to be confused by the previously introduced $\hat{\mathcal{V}}(\hat{\bf x},\hat{\bf p})=\hat{\mathcal{P}}\hat{V}(\hat{\bf x})\hat{\mathcal{P}}$).

The central result of this letter is reached by splitting the matter Hamiltonian as $\hat{H}_{\mathrm{M}}=\hat{\bf T}_{\bf R}+\hat{H}_{\mathrm{el}}$, then through a similar derivation procedure to obtain the following

		$\displaystyle\hat{\mathcal{H}}_{\mathrm{C}}=\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{T}_{R}\hat{\mathcal{P}}\hat{\mathcal{U}}+\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{P}}\hat{H}_{\mathrm{el}}(\hat{\bf p}_{r},\hat{\bf x})\hat{\mathcal{P}}\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}}$		(17)
		$\displaystyle=\sum_{j\in R}\frac{1}{2m_{j}}\hat{\mathcal{P}}\big{(}\hat{\bf p}_{j}-{\bm{\nabla}}_{j}\tilde{\bm{\mu}}\hat{\bf A}+{\bf\tilde{P}}_{j}\big{)}^{2}\hat{\mathcal{P}}+\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}},$

where the sum over $j$ only includes nuclei. In the above expression, we did not specify the choice of $\hat{\mathcal{P}}$, which could be either adiabatic or diabatic. Under the limiting case when ${\bf A}_{0}=0$ or $\tilde{\bm{\mu}}\cdot\hat{\bf A}=0$, both the $-{\bm{\nabla}}_{j}\tilde{\bm{\mu}}\hat{\bf A}$ and ${\bf\tilde{P}}_{j}$ terms (see Eqn. Resolution of Gauge Ambiguities in Molecular Cavity Quantum Electrodynamics) become 0, and $\hat{\mathcal{U}}^{\dagger}=\hat{\mathcal{U}}\to\hat{\mathcal{P}}\otimes\hat{\mathds{1}}_{\mathrm{R}}\otimes\hat{\mathds{1}}_{\mathrm{ph}}$. Thus, under a such limit, $\hat{\mathcal{H}}_{\mathrm{C}}\to\hat{\mathcal{H}}_{\mathrm{M}}+\hat{H}_{\mathrm{ph}}$ hence the matter and the cavity becomes decoupled. When using adiabatic states for the truncation, $\hat{\mathcal{P}}\hat{\bf p}^{2}_{j}\hat{\mathcal{P}}=\sum_{\alpha,\beta}(\hat{\bf p}_{j}\delta_{\alpha\beta}-i\hbar{\bf d}^{j}_{\alpha\beta})^{2}|\alpha\rangle\langle\beta|$, where ${\bf d}^{j}_{\alpha\beta}\equiv\langle\alpha|{\bm{\nabla}}_{j}|\beta\rangle$ is the well known derivative couplings. Besides these adiabatic derivative couplings, the light-matter interaction also induced additional “derivative”-type couplings, $-{\bm{\nabla}}_{j}\tilde{\bm{\mu}}\hat{\bf A}$ and ${\bf\tilde{P}}_{j}$, regardless of the electronic representation used in constructing $\hat{\mathcal{P}}$. When using the Mulliken-Hush diabatic states [27, 31] which are the eigenstates of the $\tilde{\bm{\mu}}\equiv\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}$ operator, such that $\tilde{\bm{\mu}}=\sum_{\phi}\mu_{\phi\phi}|\phi\rangle\langle\phi|$, one can prove that $\tilde{\bf P}_{j}=0$ for all nuclei. This is because that ${\bm{\nabla}}_{j}\tilde{\bm{\mu}}=\sum_{\phi}{\bm{\nabla}}_{j}\mu_{\phi\phi}|\phi\rangle\langle\phi|$, thus both $\tilde{\bm{\mu}}\hat{\bf A}$ and $[\tilde{\bm{\mu}}\hat{\bf A},\hat{\bf p}_{j}]$ become purely diagonal matrices, hence all of the higher order commutators in Eqn. Resolution of Gauge Ambiguities in Molecular Cavity Quantum Electrodynamics become zero, resulting in $\tilde{\bf P}_{j}=0$ for $j\in{\bf R}$.

$\hat{\mathcal{H}}_{\mathrm{C}}$ in Eqn. 17 provides a general and convenient expression for the Coulomb gauge molecular cavity QED Hamiltonian under the truncation of electronic states. To construct it, one can follow the general procedure: (i) obtain the matter Hamiltonian in this truncated subspace $\hat{\mathcal{H}}_{\mathrm{el}}=\hat{\mathcal{P}}\hat{H}_{\mathrm{el}}\hat{\mathcal{P}}$ and the dipole operator in this subspace $\tilde{\bm{\mu}}=\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}$ (which in principle include both the transition and permanent dipoles), (ii) constructing the operator $\hat{\mathcal{U}}=\exp[-\frac{i}{\hbar}\hat{\mathcal{P}}\hat{\bm{\mu}}\hat{\mathcal{P}}\hat{\bf A}\big{]}$, and obtain $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}$, and (iii) evaluating ${\bm{\nabla}}_{j}\tilde{\bm{\mu}}$, as well as $\tilde{\bf P}_{j}$, then couple them with $\hat{\bf p}_{R}$ and project whole term with $\hat{\mathcal{P}}$. This procedure should extend to the semi-classical treatment of light-matter interaction when considering a few level matter interacting with an intense classical laser field [23], where the substitution $\hat{\bf A}\to{\bf A}(t)$ is often made. In the Coulomb gauge, one need do this replacement in $\hat{\mathcal{H}}_{\mathrm{C}}$ instead in $\hat{\mathcal{H}}^{\prime}_{\mathrm{C}}$.

Next, we use the above general principle to derive analytical results for a model system. Without loosing generality, let us consider a molecular system within the diabatic states $\{|0\rangle,|1\rangle\}$, which represents a broad range of chemical systems [32, 33, 34]. To simplify our algebra, we will assume there is only one nuclear DOF with the coordinate $\hat{R}$ and momentum $\hat{p}_{R}$. We will further assume that $\hat{\bm{\mu}}$ is always aligned along the polarization direction of $\hat{\bf A}$ (which is $\hat{\bf e}$). Under the truncated space, $\hat{\mathcal{P}}=|0\rangle\langle 0|+|1\rangle\langle 1|$, the dipole operator is expressed as $\tilde{\mu}\equiv\hat{\mathcal{P}}\hat{\mu}\hat{\mathcal{P}}=\Delta\mu\hat{\sigma}_{z}+\bar{\mu}\hat{\mathcal{P}}+\mu_{10}\hat{\sigma}_{x}$, where $\Delta\mu=\frac{1}{2}(\mu_{00}-\mu_{11})$, $\bar{\mu}=\frac{1}{2}(\mu_{00}+\mu_{11})$, and $\mu_{\varphi\phi}(\hat{R})=\langle\varphi|\hat{\mu}|\phi\rangle$. Note that these transition and permanent dipoles are functions of $\hat{R}$. The electronic Hamiltonian in this truncated subspace is $\hat{\mathcal{H}}_{\mathrm{el}}=\hat{\mathcal{P}}\hat{H}_{\mathrm{el}}\hat{\mathcal{P}}=\varepsilon(\hat{R})\hat{\sigma}_{z}+\bar{\mathcal{V}}(\hat{R})\hat{\mathcal{P}}+{\mathcal{V}}_{10}(\hat{R})\hat{\sigma}_{x}$, where $\varepsilon(\hat{R})=\frac{1}{2}({\mathcal{V}}_{00}(\hat{R})-{\mathcal{V}}_{11}(\hat{R}))$, $\bar{\mathcal{V}}(\hat{R})=\frac{1}{2}({\mathcal{V}}_{00}(\hat{R})+{\mathcal{V}}_{11}(\hat{R}))$, and ${\mathcal{V}}_{\varphi\phi}(\hat{R})=\langle\varphi|\hat{H}_{\mathrm{el}}|\phi\rangle$ (i.e., they are $\hat{H}_{\mathrm{el}}$’s matrix elements). Using the above spin representation for $\tilde{\mu}$ and $\hat{H}_{\mathrm{el}}$, as well as the BCH identity, one can analytically show (see details in SI) that for $\hat{\mathcal{H}}_{\mathrm{C}}$ in Eqn. 17, we have

	$\displaystyle\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}=\hat{\mathcal{H}}_{\mathrm{el}}+\Big{(}\varepsilon(\hat{R})\sin\theta-{\mathcal{V}}_{10}(\hat{R})\cos\theta\Bigr{)}\Bigl{(}\sin{[\xi\hat{A}]}\hat{\sigma}_{y}$
	$\displaystyle~{}+\cos{\theta}\bigl{(}1-\cos[\xi\hat{A}]\bigr{)}\hat{\sigma}_{x}+\sin{\theta}\bigl{(}\cos[\xi\hat{A}]-1\bigr{)}\hat{\sigma}_{z}\Bigr{)},$		(18)

$\xi=\sqrt{(\mu_{00}-\mu_{11})^{2}+4\mu_{10}^{2}}$, $\tan\theta=2\mu_{01}/(\mu_{00}-\mu_{11})$, and the residual momentum is ${\tilde{P}}_{R}=\frac{1}{2}\big{(}\nabla_{R}\tan\theta\big{)}\cos^{2}\theta\big{[}\big{(}1-\cos[\xi\hat{A}]\big{)}\hat{\sigma}_{y}+\big{(}(\sin\theta)\hat{\sigma}_{z}-(\cos\theta)\hat{\sigma}_{x}\big{)}\big{(}\sin[\xi\hat{A}]-\xi\hat{A}\big{)}\big{]}$. Note that for using adiabatic states projection $\hat{\mathcal{P}}=\sum_{\alpha}|\alpha(R)\rangle\langle\alpha(R)|$, the $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}$ expression in Eqn. 18 remains the same form with $\hat{\mathcal{V}}_{01}=0$, and so does the form of $\hat{\mathcal{H}}_{\mathrm{C}}$, except for the detailed expression of $\tilde{P}$. The above result has two interesting limits. The first limit is the Rabi model (two-level atom interacts with a cavity) under the Coulomb gauge. For such an atomic-cavity system, there is no nuclear DOF $\{R,p_{R}\}$, and $\hat{\mathcal{H}}_{\mathrm{M}}=\hat{\mathcal{H}}_{\mathrm{el}}=\varepsilon\hat{\sigma}_{z}+\bar{E}\hat{\mathcal{P}}$, where $\varepsilon=\frac{1}{2}({E}_{0}-{E}_{1})=\frac{1}{2}\hbar\omega_{01}$, $\bar{E}=\frac{1}{2}({E}_{0}+{E}_{1})$, and $\hat{\mathcal{P}}\hat{\mu}\hat{\mathcal{P}}={\mu}_{10}(|0\rangle\langle 1|+|1\rangle\langle 0|)={\mu}_{10}\hat{\sigma}_{x}$. Further, by noticing that ${\mathcal{V}}_{10}\rightarrow 0$, $\xi\rightarrow 2\mu_{10}$, $\tan\theta\rightarrow\infty$, $\theta\to\pi/2$, and $\sin\theta\rightarrow 1$, $\cos\theta\rightarrow 0$, Equation 18 simplifies to $\hat{\mathcal{H}}_{\mathrm{C}}=\bar{E}\hat{\mathcal{P}}+\varepsilon\sin[2\mu_{10}\hat{A}]\hat{\sigma}_{y}+\varepsilon\cos[2\mu_{12}\hat{A}]\hat{\sigma}_{z}+\hbar\omega(\hat{a}^{\dagger}\hat{a}+\frac{1}{2})$, which is the Rabi model under the Coulomb gauge [18] that provides the consistent results of the dipole gauge Rabi model. The second limit is when these diabatic states are also Mulliken-Hush diabatic states, which means that $\mu_{10}=0$ and $\tan\theta(R)=0$, hence $\nabla_{R}\tan\theta(R)=0$, and $\tilde{P}_{R}=0$, agreeing with our previous analysis of this residual momentum.

Fig. 1 demonstrates the validity of our theory with a numerical example of a molecule couple to the cavity [1]. Here, we use the Shin-Metiu model [35] to represent the molecular system, which contains two fixed ions, one moving electron and proton (whose position is $R$), all interacting with each other through modified Coulombic potentials. The details of this model, as well as the procedure to obtain the strict diabatic states (which is not the MH diabatic representation) are provided in SI. Fig. 1a presents the diabatic potential and the matrix elements of $\tilde{\mu}$. Here, we focus on comparing the polaritonic potential energy surface $\mathcal{E}_{k}(R)$, which is defined as $\hat{H}^{\mathrm{D}}_{\mathrm{pl}}|\Phi_{k}(R)\rangle=\mathcal{E}_{k}(R)|\Phi_{k}(R)\rangle$ where $\hat{H}^{\mathrm{D}}_{\mathrm{pl}}=\hat{H}_{\mathrm{D}}-\hat{T}_{R}$ represents the polariton Hamiltonian under the dipole gauge. In the truncated electronic subspace, $\hat{\mathcal{H}}^{\mathrm{D}}_{\mathrm{pl}}=\hat{\mathcal{H}}_{\mathrm{el}}+\hat{\mathcal{U}}\hat{H}_{\mathrm{ph}}\hat{\mathcal{U}}^{\dagger}$, $\hat{\mathcal{H}}^{\mathrm{C}}_{\mathrm{pl}}=\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}^{\mathrm{D}}_{\mathrm{pl}}\hat{\mathcal{U}}=\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}+\hat{H}_{\mathrm{ph}}$, and $\hat{\mathcal{H}}^{\mathrm{C^{\prime}}}_{\mathrm{pl}}=\hat{\mathcal{P}}\hat{U}^{\dagger}\hat{H}^{\mathrm{D}}_{\mathrm{pl}}\hat{U}\hat{\mathcal{P}}$. Note that the analytical results of $\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}_{\mathrm{el}}\hat{\mathcal{U}}$ is expressed in Eqn. 18, whereas the details of other expressions are provided in SI. The matrix elements of these Hamiltonians are evaluated with the two electronic diabatic states and a large number of Fock states, and diagonalizing this matrix gives $\mathcal{E}_{k}(R)$. Fig. 1b presents $\mathcal{E}_{k}(R)$ with $A_{0}=0.2$ a.u. and $\hbar\omega_{\mathrm{c}}=3$ eV (such that the light and matter excitations are in resonance at $R=0$). While the $\hat{\mathcal{H}}^{\mathrm{D}}_{\mathrm{pl}}$ (solid) and $\hat{\mathcal{H}}^{\mathrm{C}}_{\mathrm{pl}}$ (dotted) give identical results throughout all range of $R$, $\hat{\mathcal{H}}^{\mathrm{C^{\prime}}}_{\mathrm{pl}}$ (dashed) gives inconsistent results and breakdown gauge invariance. Fig. 1c-d presents $\mathcal{E}_{k}-\mathcal{E}_{0}$ at $\hbar\omega_{\mathrm{c}}=3$ eV as a function of the field strength $A_{0}$, at $R=0$ (resonance condition) and $R=-1.74$ a.u. (detuned), respectively. Again, the results from the Coulomb gauge and dipole gauge agree with each other exactly through out the entire range of the field strength, whereas simple state truncation on the Coulomb gauge QED Hamiltonian breaks the gauge invariance, especially in the ultra-strong coupling regime [18, 22].

In conclusion, we lay out the fundamental theoretical framework for the molecular cavity QED by presenting the general procedure to obtain the Coulomb gauge Hamiltonian that provide the consistent results from the dipole gauge Hamiltonian, under the same level of electronic state truncation. Using a finite level of projection on $\hat{H}_{\mathrm{D}}$ often provides accurate results of the polariton eigenspecturm, whereas $\hat{H}^{\prime}_{\mathrm{C}}$ often introduce the gauge ambiguity (especially in the ultra-strong coupling regime). Instead, $\hat{\mathcal{H}}_{\mathrm{C}}$ in Eqn. 17 resolves such gauge ambiguity by proving consistent results as $\hat{\mathcal{H}}_{\mathrm{D}}$. Investigations based upon the Coulomb gauge [36] should consider using $\hat{\mathcal{H}}_{\mathrm{C}}$.