Gauge-Independent Emission Spectra and Quantum Correlations in the Ultrastrong Coupling Regime of Open System Cavity-QED

A quantum dipole interacting with an optical cavity is one of the key models in cavity quantum electrodynamics (cavity-QED). To treat this system theoretically, the typical approach is to truncate the dipole to two levels. However, it has been shown that in the ultrastrong-coupling regime, this truncation naively destroys gauge invariance. By truncating in a manner consistent with the gauge principle, we introduce master equations {for open systems} to compute gauge-invariant emission spectra, photon flux rates, and quantum correlation functions which show significant disagreement with previous results obtained using the standard quantum Rabi model. Explicit examples are shown using both the dipole gauge and the Coulomb gauge.


I. INTRODUCTION
The intricate interactions between light and matter allow one to observe drastically different behavior depending on the relative strength of the light-matter coupling. In the weak-coupling regime, the losses in the system exceed the light-matter coupling strength, and energy in the system is primarily lost before it has the chance to coherently transfer between the matter and the light. Accessing this regime experimentally has allowed for breakthroughs in quantum technologies such as single-photon emitters [1][2][3][4]. Beyond weak-coupling, in the strong-coupling regime the rate of decoherence is smaller than the rate of excitation exchange, allowing for the observation of vacuum Rabi oscillations: the coherent oscillatory exchange of energy between light and matter. The strong-coupling regime has helped initiate a second generation of quantum technologies [5,6]. See Fig. 1 for a simple schematic of a typical cavity-QED system with system-bath leakage. Figure 1. Schematic of a generic cavity-QED system. The optical cavity mode has quantized energy levels (in blue), with a decay rate κ. The matter system is a truncated TLS (in red), with a possible spontaneous emission decay rate γ. The two systems have a coherent coupling strength g. A coherent laser (in orange) drives the system with Rabi frequency Ω d .
Around 2005, the "ultrastrong-coupling" (USC) regime was predicted for intersubband polaritons [7]. This regime is characterized not by still lower rates of decoherence, but by a coupling strength that is a comparable fraction of the bare energies of the system. The dimensionless parameter η = g/ω 0 (i.e., the cavity-emitter coupling rate divided by the transition frequency) is used to quantify this coupling regime for cavity-QED. Typically, USC effects are expected when η 0.1, at which point the rotating wave approximation (RWA) used in the weak and strong regimes becomes invalid. Reported signs of USC emerged in 2009 with experiments involving quantum-well intersubband microcavities [8], achieving η ≈ 0.11. Terahertz-driven quantum wells have also demonstrated USC effects [9], and similar effects have been exploited to achieve carrier-wave Rabi flopping with strong optical pulses [10][11][12]. To date, many different systems have exhibited USC [13,14]. Recently, using plasmonic nanoparticle crystals, η = 1.83 has been achieved, with potential to lead to η = 2.2 [15].
With experiments pushing the normalized coupling strength continuously higher, the interest in USC effects also continues to grow, helping to improve the underlying theories of light-matter interactions [16], even at arbitrarily high coupling strengths [17]. There have also been various predictions made about what novel technologies USC will bring about, including modifications to chemical or physical properties of various systems caused by their USC to light [7,18], and the potential to create faster quantum gates and gain a high level of control over chemical reactions [13]. To push these advancements forward, it is essential to have a fundamental understanding of the physics involved with these systems and to accurately connect to experimental observables.
The cornerstone model in cavity-QED is a two-level system (TLS) interacting with a quantized cavity mode [19]. This model has been applied to atoms [20][21][22][23], quantum dots [24][25][26][27], and circuit QED [28][29][30][31]. Outside the USC regime, this model is typically represented by the canonical Jaynes-Cummings (JC) Hamiltonian [32], which makes a RWA and can be easily diagonalized. In the USC regime, however, it is necessary to retain counter-rotating terms, giving rise to the quantum Rabi model (QRM) [13,14,33]. By detecting resonance fluorescence of light emitted from the cavity as quantified by the first-order degree of coherence correlation function (CF), the spectral content of these cavity-QED models can be explored, while the second-order intensity CF is fundamental to understanding the photon statistics as probed by intensity interferometry.
The main contribution of this work is to present a self-consistent and unambiguous way to model observables in the USC regime of open system cavity-QED. Apart from addressing the subtle (and unknown) effects of dissipation, and excitation, and input-output, we show the striking influence of modelling experimentally relevant observables such as the emission spectra and quantum correlation functions. We also show how and why the form of the system-bath interactions matters, yet the form is gauge invariant, if-and only if-treated properly (in contrast to the standard master equation approaches) using gauge invariant master equations. We show equivalence between dipole gauge and Coulomb gauge master equations, if one applies gauge corrections in a consistent way, and we also demonstrate the drastic failure of currently adopted master equations in the USC regime. Our framework and formalism, to the best of our knowledge, constitutes a first way to do this, and can thus be applied to a wide range of measurements in the USC regime for open systems.

II. GAUGE INVARIANCE AND SYSTEM-RESERVOIR INTERACTIONS
It was recently shown that extra care is needed when constructing gauge-independent theories [34], for computing experimental observables for suitably strong light-matter interactions. This development started with a series of papers dealing with so-called gauge ambiguities in the USC regime [35][36][37]. As a U (1) gauge theory, different gauges in QED manifest in different representations of the Hamiltonian of a given system, but these should be unitarily equivalent and give rise to equivalent physical observables. Without proper care, gauge invariance of cavity-QED theories can break down when considering USC [38]. This is due to the truncation of the matter system's formally infinite Hilbert space to the two lowest eigenstates in forming the TLS-only keeping an infinite number of energy levels formally preserves gauge invariance [39]. Consequently, previous model predictions in the USC regime can be ambiguous since the predictions are impacted by the choice of gauge. While this issue has been known in general for several decades [40], only recently was this specific problem presented as rather insurmountable [38]. However, the issue has been resolved by using a self-consistent theory at the system Hamiltonian level [37,41], restoring gauge invariance to the theory for systems with a finite Hilbert space.
Despite this, additional subtleties occur in the USC regime regarding the interaction of the cavity-QED system with its environment. To connect to experiments, one also requires an input-output model of dissipation from the cavity to external modes, requiring an open-system model of cavity-QED. In the USC regime, complications arise with this input-output formalism associated with approximations typically made outside of the USC regime. These complications originate from the hybridization of light and matter that occurs in USC, and as such the quanta of excitations inside the cavity-QED system have different quasiparticle representations than the photons actually emitted from the system. Moreover, the separation of operators into light and matter components becomes highly gauge-specific in the USC regime, and proper care must be taken to ensure self-consistency.
To fully synthesize these considerations with the restoration of gauge invariance, we present a dissipative and gaugeinvariant master equation model, which is required to properly describe experimentally-observable quantities arising from output channels of the cavity. Key experiments to probe such observables include resonance fluorescence and two-photon detection schemes, and we make a direct connection to both of these. We also show how previous QRM master equations in the USC regime are ambiguous in general as they produce gauge dependent results for observables, and we show how to fix such problems. Moreover, our theories can be used to explore the precise form of the system-bath interactions, which in fact yield different experimental signatures in the USC regime.

III. MODEL
In the dipole gauge (namely, the multipolar gauge in the dipole approximation), we can write the system Hamiltonian, using the QRM, as ( = 1) where ω 0 (ω c ) is the TLS (cavity) transition frequency, σ + (σ − ) is the raising (lowering) operator for the TLS, and a † (a) is the cavity mode creation (annihilation) operator; g is the TLS-cavity coupling strength. We take ω c = ω 0 throughout. In contrast to the Coulomb gauge, straightforwardly truncating the dipole in the light-matter interaction to a TLS subspace does not break gauge invariance in the dipole gauge [37]. Making a RWA on Eq. (1) (i.e., neglecting counter-rotating terms a † σ + and aσ − , which do not conserve excitation number), yields the simpler JC Hamiltonian. Outside of the USC regime, the usual approach to include dissipation is with a Lindblad master equation [42], where ρ is the reduced density matrix. The dissipation term, Oρ and κ is the cavity photon decay rate. Since dissipation is usually dominated by cavity decay, we neglect direct TLS relaxation and pure dephasing [43,44]. However, the theory of how to include TLS dissipation is discussed in Appendix A. 4 The Lindbladian can be derived by following the typical approach in which one neglects the TLS-cavity interaction when considering the coupling of these systems to the environment [29]. However, when moving into the USC regime, this approach fails, and the Lindbladian must be derived while self-consistently including the coupling between the subsystems. For sufficiently strong subsystem coupling, transitions occur between dressed eigenstates of the full Hamiltonian rather than between eigenstates of the individual free Hamiltonians [43].
In the USC regime, the system has transition operators |j k| which cause transitions between the dressed eigenstates of the system {|j , |k }. To obtain these transitions for the cavity mode operator, we use dressed operators [43], and x − = (x + ) † , where the sum is over states |j and |k , with ω k > ω j , C jk = j| Π C |k , and we neglect thermal excitation effects; Π C is an operator which couples linearly to dissipation channel modes which we assume proportional to the cavity electric field operator such that Π C = i(a † − a). We then replace L bare in Eq. (2) with L dressed ρ = κ 2 D[x + ]ρ, to arrive at the dressed state (DS) master equation. One can also use a generalized master equation to capture coupling to frequency-dependent reservoirs [43,45] See Appendix A for a derivation of the generalized master equation, and Sec. V for an example application using an Ohmic bath.
Beyond this dressing transformation, it has been shown that there exists a potential gauge ambiguity in the electric field operator which causes further problems when computing observables in the USC regime [37]; namely, Π C corresponds to the Coulomb gauge electric field, but the QRM Hamiltonian is derived in the dipole gauge. The gauge transformation from the Coulomb gauge to the dipole gauge is generated by a unitary transformation, which for the restricted TLS subspace is given by the projected unitary operator [37] U = exp(−iη(a + a † )σ x ). The photon destruction operator transforms as a → UaU † = a + iησ x [34]. Thus, to "gauge-correct" the master equation in the dipole gauge, we conduct the dressing operation as above, but with where we take C jk = j| UΠ C U † |k = j| Π D |k = j| i(a † − a) + 2ησ x |k ; see Appendix A for a derivation of the master equation in the dipole and Coulomb gauges and their equivalence.
To study the quantum dynamics and spectral resonances, we excite the system with an incoherent pump term, P inc D[x − GC ]/2, or with a coherent laser drive, H drive (t) = (Ω d /2)(x − GC e −iωLt + x + GC e −iωLt ), added to H QR , where Ω d is the Rabi frequency and ω L = ω c is the laser frequency; thus, H S = H QR + H drive . Note that the QRM with a coherent drive is time-dependent and oscillates around a pseudo-steady-state. In addition, because of the driving laser, the periodic nature of the system Hamiltonian means that in principle the QRM spectra, already quite rich, are modified further; however, we use Ω d g, and neglect the influence of the coherent drive on the system eigenstates. The first few (lowest) energy eigenvalues are plotted for the QRM (dipole gauge) and JCM in Fig. 2(a) for a range of normalized coupling strengths. Three transitions are shown, which we will refer to below.

IV. GAUGE INVARIANT OBSERVABLES
We first define the system excitation number, and a quadrature operator matrix element squared, which is proportional to the photodetection rate of cavity-emitted photons from the |j → |k transition [41]. In Fig. 2(b), we show N cav versus η, using incoherent driving (cf. Fig. 3), where the solid curves show the effect of gauge corrections. Equivalent gauge-corrected results are obtained in the Coulomb gauge. With the correction, the population saturates , while the uncorrected population continues to increase superlinearly, and jumps when states 2 and 3 cross in energy, potentially related to the photon blockade [46]. With gauge corrections, we see a strong influence from the TLS operator physics. In Fig. 2(c), we show |P jk | 2 for the relevant transitions which are, for weak excitation, proportional to the transition linewidths; again, the solid lines show the gauge corrected results. Note that the corrected dipole gauge quadrature operator (Π D = i(a † − a) + 2ησ x ) causes a major modification of the transitions, significantly impacting their behavior in the nonperturbative regime.
In Appendix E, we give analytical insight into these quadrature matrix elements using a Bloch-Siegert (BS) transformation, which analytically (to lowest order in η) predicts the following changes with gauge correction: |P I | 2 =1/4(1+3η/2) → 1/4(1−5η/2), and |P III | 2 =1/4(1−3η/2) → 1/4(1+5η/2), causing a reversed asymmetry with gauge corrections. Physically, this asymmetry arises from the BS shift of cavity and TLS resonances giving rise to photon-like and atom-like polariton branches; the composition of the Π operator (which is affected by the gauge correction) ultimately determines which state is more cavity-like, and thus has a greater decay rate (see Appendix E for details).
Next, we define the cavity-emitted spectrum, where x ± GC,∆ =x ± GC − x ± GC and Ω=ω − ω L . Beyond the spectrum, which uses a first-order quantum CF, we also compute the normalized second-order quantum CF, which quantifies the likelihood of a photon being detected at (t + τ ) if one was detected at t. We also introduce the time-averaged g (2) (τ ) = t1+T t1 g (2) (t, τ )dt/T , where t end is an arbitrary time point at which the system has reached the pseudo-steady-state and T is the period of oscillation (see Appendix D). Note that without the gauge-correction, we use the uncorrected (corresponding to a Coulomb gauge representation) x ± , x ± ∆ for computing the observables, and x ± for incoherent or coherent driving (see Appendix A). All calculations use Python with the QuTiP package [47,48].
For weak incoherent pumping, Fig. 3 compares the computed spectra with and without the gauge correction (DGC: dipole-gauge-corrected and DG: dipole-gauge, respectively), for η ranging from 0.05 (strong coupling) to 0.5 (USC). For relatively small η = 0.05, the DGC (with gauge correction) spectra already begin to deviate from the DG spectra (usual QRM master equation solution).
With increasing η, notably, the DGC and DG spectra are substantially different above η = 0.1: the DGC spectra still show a reversed asymmetry, with a significant narrowing of the lower polariton resonance (I) and a broadening of the upper polariton resonance (III); the ratio of higher-lower polariton peak areas under weak excitation changes from 1 − 3η + O(η 2 ) to 1 + 5η + O(η 2 ) with gauge correction -a dramatic change even for η < 0.1 (see Appendix E). These peaks can be identified as resulting from the |1 → |0 (olive arrow on Fig. 2(a)) and |2 → |0 (brown arrow) transitions, respectively. Since |P j,k | 2 contributes to photon emission directly through the κ decay channel [41], the narrowing (broadening) of peak I (III) with increasing η can be explained with Fig. 2(c). Without the correction, the opposite trend is observed, which is again consistent with Fig. 2(c) (dashed lines). At η = 0.5, there is also a noticeable resonance (II) around ω = 0.8g, showing a deep mixing of the TLS and cavity dynamics in the USC regime. We can identify this energy difference with the |3 → |1 transition, pink arrow on Fig. 2(a), which also has reduced broadening with η, cf. Fig. 2

(c).
We have shown how the gauge correction manifests in modified linewidths and drastically different spectral weights in comparison to the usual QRM-even so far as to result in a complete reversal of the asymmetry predicted from a nongauge-corrected model [49] (Fig. 9, η = 0.5). We now demonstrate how this gauge correction manifests in the Coulomb gauge. To do this, we display results for the cavity-emitted spectrum and CFs with coherent and incoherent pumping, using the discussed dipole gauge and the Coulomb gauge master equation.  In the Coulomb gauge, the standard system Hamiltonian for the QRM is [37] where g C = g D ω 0 /ω c and D is the strength of the diamagnetic term. Using the Thomas-Reiche-Kuhn sum rule [41], then D ≥ g 2 C /ω 0 , and for our simulations we take D = g 2 C /ω 0 . Thus, with ω 0 = ω c and g D ≡ g, we have D = η 2 ω 0 . Unfortunately, this form does not satisfy the gauge principle, and produces the wrong eigenenergies and eigenstates in the USC regime [35,37,41]. Instead, the corrected Coulomb gauge uses a different system Hamiltonian [37], which contains field operators to all orders, and the C superscript indicates we are using the corrected form for the system Hamiltonian. In the Coulomb gauge, the gauge-invariant dissipator term is (see Appendix A) where x + C = j,k>j C C jk |j k| with C C jk = j|Π C |k , and we now compute the dressed states in the Coulomb gauge, using both uncorrected and corrected forms of the system Hamiltonian. Figure 4 (top) shows the coherent and incoherent spectra at η = 0.5, showing that the gauge correction results in a profound effect in either case. For coherent driving, using Ω d = 0.1g, there is a significant sharpening of the resonances. The Coulomb gauge result without the gauge correction corresponds to a minimal coupling Hamiltonian naively truncated to a TLS, which results in incorrect energy levels for the dressed-state master equation [37]. This effect of having the incorrect energy levels and eigenstates is clearly shown in the uncorrected Coulomb gauge results in Fig. 4, which is especially wrong with coherent pumping, since the system is effectively being pumped off resonance (because of the diamagnetic term). For coherent pumping, additional Rabi field strengths are shown in Appendix C, where we also show simulations with and without a RWA for the pump field.
Next, in Fig. 4 (bottom), we examine the second-order coherence, which is important for characterising the generation of non-classical light. In all cases shown, we observe photon bunching at short time-delays. With the gauge correction, there is a significant reduction in the level of bunching, and the usual USC master equations significantly overestimate the bunching characteristics. Moreover, the dynamics are qualitatively different, and thus the non-GC master equations results clearly fail in the USC regime. In all cases, we confirm full agreement between the corrected dipole gauge and corrected Coulomb gauge results, since these are the correct gauge invariant solutions, and thus produce identical results. In the simulations above, for simplicity, we used a flat density of states (DOS) for the spectral bath function; namely, the DOS was assumed to be constant relative to the energy scale of the resonances. This helps to better identify intrinsic spectral asymmetries related to gauge correcting.
For completeness, here we explicitly show an example numerical solution without invoking the approximation that κ(ω) is frequency independent. Specifically, we compute the emitted spectra when κ(ω) = κ as well as κ(ω) = κω/ω c (Ohmic bath). We use the same example as in Fig. 3 with incoherent driving at η = 0.5. These numerical solutions are obtained from the generalized master equation (A9), described in Appendix A.
As can be seen in Fig. 5, clearly the form of the spectral bath function does not affect any of our general conclusions, as the gauge correction is, in both cases, dramatic, and of course produces exactly the same result for both the dipole gauge and the Coulomb gauge. To be clear, if we plot these together, then they are indistinguishable, which also confirms that our numerical results are well converged in terms of basis size and time steps.

VI. CONCLUSIONS
We have presented a gauge-invariant master equation approach and calculations for the cavity emission spectra in the USC regime, and shown how the usual QRM in the dipole gauge fails, yielding effects that are just as pronounced (or even more pronounced) as counter-rotating wave effects in this regime. We have demonstrated how the gauge correction significantly affects the intensity CF and cavity excitation number. We have also shown how the gauge correction modifies results in the Coulomb gauge compared to typically used models. Apart from yielding new insights into the nature of system-bath interactions, and presenting gauge-invariant master equations that can be used to explore a wide range of light-matter interaction in the USC regime, our results show that currently adopted master equations in the USC regime produce ambiguous results since they do not satisfy gauge invariance.
While we have shown explicit results for the cavity spectrum and intensity CF, the gauge correction causes profound effects on any observable that is computed from the master equations in the same coupling regimes. The nature of the system-bath coupling is also very important, which must also be related to the quadrature coupling to the external fields and the observables to ensure a gauge invariant master equation. For example, it may be more appropriate to use Π C = a + a † (vector potential coupling) rather than Π C = i(a † − a) (electric field coupling) for the interaction (in the Coulomb gauge), or some linear combination of the two; this change affects the dissipators, incoherent pumping, and coherent excitation in a way that still yields gauge-independent results, but the observables are different. By unitary equivalence, the form of the quadrature coupling used in the system Hamiltonian is thus also not arbitrary, which is in stark contrast to the JC model, where both these coupling forms yield identical results. These two coupling forms are widely used in the USC literature and are assumed to lead to the same result; however, they differ significantly, which reinforces the need, highlighted recently [50,51], to go beyond the usual phenomenological formulation of systemenvironment coupling Hamiltonians in the USC regime of cavity QED in favor of a general fundamental microscopic derivation. Solutions to such problems can likely be rigorously addressed using quantized quasinormal modes [52][53][54][55][56], which even apply to cavities and media in the presence of gain [57,58]. Let us first consider a general bath (or reservoir) that interacts with the system of interest (e.g., the cavity mode) weakly, as shown schematically in Fig. 1 of the main text. The bath is described in the usual way by a collection of harmonic oscillators ( = 1), where b k and b † k are bosonic annihilation and creation operators. A simple model for a single cavity interacting with the bath can be written as follows: where λ k represent the coupling strengths (assumed real), which are model specific, and Π is a gauge-dependent system operator linear in the canonical quantization variables, the form of which we specify based on physical considerations in the following sections. In the interaction picture, we havẽ In the dressed-state basis, which is necessary to use in the ultrastrong coupling (USC) regime (as the standard dissipator fails, as discussed in the main text), we can express the lowering operators of the system excitations from where and ∆ jk = ω j − ω k , such that and any C jj terms uniformly vanish due to the parity symmetry of the quantum Rabi model. The bath operators can also be written asB Thus we can write [43],H where we have dropped all terms which oscillate at a frequency equal to a sum of positive system and reservoir frequency components which do not ultimately contribute to the master equation we will derive. Applying a Born-Markov approximation, assuming continuous bath frequencies, a zero temperature approximation (namely, neglecting thermal excitation and taking the bath to be in the vacuum state), and neglecting any Lamb-like renormalization of the quantum Rabi Hamiltonian parameters, one can derive a generalized master equation [43], which takes into account the dressed-states' coupling to all the relevant baths for each system operator: where the cavity dissipator term is The dressed-state operators, X ± , decomposed in a basis of energy eigenstates with respect to H QR , are defined through where ω = ω k − ω j > 0 and X − = (X + ) † . Note we can also derive a similar generalized master equation for other system decay channels (i.e., TLS losses), but below we concentrate on the cavity operators and relevant system-reservoir interactions, though we also briefly discuss the TLS-bath interactions. One can employ any representative bath functions for the cavity reservoir, J c (ω) = g c (ω)|λ(ω)| 2 , where g c (ω) is the bath density of states (DOS), and the decay rates are subsequently defined from Γ c (ω) = 2πJ c (ω) = 2πg c (ω)|λ(ω)| 2 . (A12) Thus, for example, in the case of an Ohmic bath (J c (ω) ∝ ω), then Γ c (ω) = γ c ω/ω c , where γ c ≡ κ. Finally, assuming a relatively flat bath function with respect to the frequency differences of interest (we will relax this approximation later, in Sec. V), so that Γ c (∆ jk ) ≈ κ, with κ = κ(ω 0 ) (nominal cavity decay rate) over the energy scales of interest, we obtain where with and the usual Lindblad superoperator term, Equation (A13) is the standard dissipator form in the USC regime for the dressed-state master equation. Without any consideration of gauge, one might naively take Π = i(a † − a) (the form we take for the non gauge-corrected form of the dipole gauge model in the main text), or perhaps Π = a + a † ; however, in contrast to usual cavity-QED systems outside of the USC regime, these choices give rise to different observables, and furthermore, lead to gauge-dependent results (in the USC regime). We address this explicitly in the following sections, and review how the breaking of gauge invariance that can be introduced by truncation to a TLS subspace is "gauge corrected" in the dipole gauge by modification of the Π operators from their naive form, and in the Coulomb gauge by modifying the Hamiltonian [37, 59, 60].

Gauge-invariant Master Equation in the Dipole Gauge
To specify our dissipation model, we must assign a specific form to the gauge-dependent system operator Π, and thus we must consider how relevant physical quantities are represented in each gauge. In the dipole gauge, it is the displacement field that is expanded in terms of bosonic creation/destruction operators, and not the transverse electric field as in the Coulomb gauge. The relevant field operator is F = D/ 0 b (r) [60][61][62], where D is the displacement field and b is the dielectric constant that the TLS is embedded (e.g., for free space this is 1). Importantly, D also includes a contribution from the TLS dipole field through the polarization. The single cavity mode field-TLS interaction is then where g c = −i ωc 2 0 µ · f c (r 0 ) and r 0 is the dipole (TLS) location. The cavity mode amplitude is real (corresponding to a normal mode), and defining g c = −ig, where g is real, then which is projected onto a two level subspace. Relating F to E (the electric field operator), and using a TLS coupling for the source of the polarization, results in the electric field being expanded in terms of transformed cavity operators, a = a + iησ x , such that [59] where A(r) = 1/2 0 ω c f c (r), the amplitude of the vector potential field. Note that the explicit coupling between a and σ x here is a direct consequence of a strict single mode approximation. If we assume a weak coupling between the cavity electric field and reservoir modes, the system-reservoir coupling takes the gauge-corrected form: where we let Π D denote the system operator to be inserted into Eq. (A2) in the dipole gauge, and we assume the b k are unchanged. Note as mentioned above we could also consider a linear coupling between the vector potentials of the cavity and reservoir fields, such that Π ∝ a + a † (which is manifestly gauge invariant in form). Outside of the USC regime, these couplings produce identical results within the rotating wave approximation (RWA), and are often assumed to be interchangable; however, in our simulations, choosing this form of system-reservoir coupling leads to significantly different observables (similar conclusions were drawn in Ref. [41]). This is because in the JC model, x + ∝ a + O(η) in any gauge, and any change in the phase of a is compensated for in the Lindblad term which pairs a and a † . In the USC regime, the counter-rotating terms in the QRM ensure that the dissipator is not invariant under such a change. Noting that a coupling of the form a + a † can be transformed into i(a † − a) by the unitary transformation U = exp [−i π 2 a † a], an important consequence of this is that a coupling in the (dipole gauge) QRM of the form ig(a † − a)σ x is not equivalent to one of the form g(a + a † )σ x when dissipation is to be considered, despite what is commonly assumed. In the USC regime, the gauge and form of the dissipators must be properly considered in conjunction with the Hamiltonian in order to ensure gauge-invariant observables. The only symmetry in the dissipative QRM is then that of parity symmetry, which ensures that the overall sign of any couplings terms can be changed. For this work, we restrict ourselves to electric field-like couplings such that Π C = i(a † − a).
Following the same steps as above, we obtain the dipole gauge result for the dressed-state dissipator, where with and the QRM system Hamiltonian is where we use ω 0 /2σ z instead of ω 0 σ + σ − (in the main text) to compare with the Coulomb forms below. The 'boxed' equations ( (A21),(A23),(A24)) represent the gauge-corrected dissipator and QRM system Hamiltonian in the dipole gauge. Note that as in the main text, to compute optical observables emitted from the cavity in this gauge we should also apply the gauge correction (i.e., use the x ± D operators), to be consistent with input-output theory [63]. More formally, before we switch to the Coulomb gauge, we should also identify g ≡ g D as being the TLS-cavity coupling rate in the dipole gauge. In the main text, we let quantities without explicit subscript/superscript refer to the dipole gauge, and in particular, Π = i(a † − a) without any gauge correction, and Π = i(a † − a) + 2ησ x with the proper gauge correction (subscripts GC in main text), both in the dipole gauge.

Gauge-Invariant Master Equation in the Coulomb Gauge
As discussed in the main text, in the Coulomb gauge Hamiltonian, we have the following system Hamiltonian for the QRM [37] where g C = g D ω 0 /ω c and D is the diamagnetic term. For ω 0 = ω c and g D ≡ g, we can use D = η 2 ω 0 as a lower bound [41].
In the USC regime, Eq. (A25) does not produce the same eigenenergies as Eq. (A24), since it fails to respect the gauge principle. Instead, the properly gauge-transformed form, which does produce the same eigenenergies, is given by the following system Hamiltonian for the QRM [37]: where U = exp(−iη(a + a † )σ x ), as in the main text. Notably, H C QR contains field operators to all orders. In the Coulomb gauge, the form of the electric field operator is proportional to i(a † − a), assuming the same bath interactions, and thus we have Π C = U † Π D U = i(a † − a), and the system-bath coupling is written as where λ C k is cavity-bath interaction in the Coulomb gauge. Following similar steps to before, we obtain the gauge-invariant dissipator term: where with and now one uses the dressed states in the Coulomb gauge, namely using H C QR . Note, to include an arbitrary spectral function, then we use and J c and λ(ω) are identical in the dipole gauge and Coulomb gauge. Below we show this explicitly for the case of an Ohmic bath. The three boxed equations ((A26),(A29),(A30)) represent the correct form for the Coulomb gauge master equation to give equivalent results to the dipole-gauge forms, which we prove in Sec. B and show explicitly in Fig. 3 of the main text.
Note also that the expressions for C jk can be rewritten via a sum rule [41]. For example, in the Coulomb gauge [41]: and in the dipole gauge:

Two Level System (TLS) Dissipator
Next, for completeness, we discuss the TLS dissipator (whose contribution is negligible in our simulations) and again show equivalence between the dipole gauge and Coulomb gauge. The σ x is invariant when transformed through the gauge correction, and thus there is no change; namely we simply have: in either gauge. However, for a specific model for the spontaneous emission decay, a more realistic model would include frequency dependent reservoirs representative of the free space emission channel (such as Ohmic). Later we show an example of how to incorporate such effects for the dominant cavity model decay channel, and it is easy to also do this for the TLS decay, if required.

Incoherent Pump Term
In a standard master equation, the incoherent pumping for the cavity mode is usually written as a reversed Lindblad decay process [64,65], which in a dressed-state decomposition is This type of excitation can be derived by input-output theory with (for example) non-vacuum inputs (e.g., a thermal state with T = 0) [63]. Thus to be consistent with our dissipation channels and the microscopic form of the systemreservoir coupling, we choose

Coherent Pump Term
Next we present a derivation of the coherent drive term in the Hamiltonian, H drive (t). We consider the general interaction picture Hamiltonian of Eq. (A2), where the gauge of this interaction is ultimately determined by the form ofS(t), which is left general here. In deriving the master equation models, we formally considered the reservoir to be in a multimode vacuum state |0 as t → 0. To model coherent driving at the level of a system-reservoir approach, where the input drive (laser field) is not significantly impacted by the dynamics of the cavity-QED system, we can instead assume the reservoir to be in a multimode (or approximately single mode) coherent state, with the input condition: where U c = k D k (β k ), D k (β k ) = exp β k b † k − β * k b k is the displacement operator, with D k (β k ) |0 = |β k (where all k = k remain in the vacuum state), and β k are substantial only for wavevectors k around the laser resonance. Since U c is unitary, we can apply a unitary transformation to the system plus reservoir density operator and Hamiltonian. Within the Born-Markov approximation, we haveρ S+B =ρ S ρ B , thus we apply the unitary transformationρ S+B → U † cρS+B U c , andH SB → U † cHSB U c . The effect of this is merely to take b k → b k + β k within the interaction picture. Thus we havẽ B(t) →B(t) + k λ k β k e −iω k t , and Since the new term in Eq. (A42) only depends on the system operators, we can call itH drive and consider it part of the system Hamiltonian. Moving back to the Schrödinger picture: In this new frame, the bath is in the multimode vacuum state |0 . Thus, the master equation can be derived in exactly the same manner as before, with the only difference in the equations being the addition of H drive (t) in the system Hamiltonian. Typically, we can make a RWA for this term, as we have separated positive and negative frequency components, but in principle we leave this general as the RWA could break down for ultrastrong coherent driving (however, in this regime, a Floquet master equation would be more accurate).
To transform Eq. (A43) into an effective single-mode drive, we move to a continuous frequency representation: Since β c (ω) is only nonzero around a very narrow window around ω = ω L (the laser center frequency), we have The form of β c (ω) is not important provided it is sharply peaked around ω = ω L ; for concreteness we can assume a Lorentzian form: where W 0 is the FWHM of the laser beam, and we find: where we have defined We assume that the laser linewidth W 0 is small enough such that the drive remains coherent over any experiment of interest. We can also choose Ω d to be real without loss of generality, as the phase factor can be absorbed into the initial phase of the drive which is not relevant. Thus we find the form used in the main text: using identical system operators as in the dissipators and incoherent pump terms. Also note, the coherent drive should not be too strong to invalidate the dressed-state representation for the system Hamiltonian, namely Ω d g, and in this regime one could make a RWA for the pump term such that H In the main text, we use this RWA pumping term and also consider a resonant drive where ω L = ω c . For completeness, in Sec. C below, we also show example calculations with and without a RWA on the pump term, and confirm that they yield essentially identical spectra, as expected (i.e., for the stated approximations).

Appendix B: Equivalence between the dipole gauge and Coulomb gauge master equations and gauge invariant expectation values
Naturally, any observables from a unitarily transformed quantum master equation should be gauge-independent; we include this section primarily to show that the x ± operators transform in the way one might expect. By gauge independent expectation values, we mean expectation values that do not and should not depend on the choice of gauge.
Ultimately, for any gauge-dependent Hermitian operator O D and O C corresponding to a physical observable, the expectation value should be gauge-independent, sȯ Beginning with the evaluation in the Coulomb gauge, Clearly, the evolution is gauge-invariant if, We have (explicitly noting the gauge of each state):   . Cavity spectra and g (2) with Ω d = 0.1g coherent pumping, using a full cosine excitation (left), and a RWA for the pumping term (right, as also shown in Fig. 4 of the main text). Solid lines show the gauge corrected master equation results. spectra in a truncated basis space. With 24 dressed states in the truncated space, we observe negligible excitation in the highest states and numerically converged results (i.e., additional dressed states make no change to our simulations and results). The eigenenergy simulation in Fig. 2(a) of the main text was conducted with 200 photon states to ensure accurate numerical convergence. Longer times are required for simulations at higher η and as such, the simulation time was between 150/g and 550/g throughout, with 20 time-steps in each period of the pseudo-steady-state oscillation. This was also carefully checked to be sufficient. Numerical calculations were performed using QuTiP under Python [48].
With coherent excitation, numerical calculations of the spectra outside the rotating wave approximation require some care. Specifically, the t integral in the spectrum definition (Eq. (7) of the main text) was completed over the last ω L time period so as to ignore turn-on dynamics, after ensuring that the system had reached its pseudo-steady-state (namely, after it evolves to a continuous oscillation dynamic with no change). Consequently, there is a potential issue with computing a Fourier transform of an oscillating function over a finite range. This is commonly done for computing USC spectra but is rarely discussed. Formally, the Fourier transform of a sin ω 0 t function over a finite range a is proportional to the difference of two sinc functions at ±ω 0 , shown explicitly below, When extending this time sampling range to infinity, the Fourier transform tends towards the sum of two Dirac delta functions, This can have a significant effect on both the total and coherent spectrum, but in all our simulations the incoherent spectrum with coherent driving is unaffected, as it does decay to zero for large time delays, and performing the quantum regression theorem over only a single period is thus adequate for our case. For example, we have checked that we obtain the same result when integrating over ten periods.

Quantum Regression Theorem
To calculate the two-time correlation function in the spectrum definition (Eq. (7) of the main text), we make use of the quantum regression theorem, where U (t+τ, t) is the total (system + environment) unitary evolution operator such that A(t+τ ) = U † (t+τ, t) [A(t)] U (t+ τ, t) for an operator A(t). Within the Born-Markov approximation, the implementation of the quantum regression theorem for Eq. (D3) is as follows: find the reduced density matrix at t, multiply on the right by x − ∆ (0), evolve this new operator from t to (t + τ ) with the master equation to form the effective density matrix, and finally take the expectation value of x + ∆ (0) with respect to this effective density matrix. Note that x − where the second term must be evaluated at t, so x − ∆ (0) does implicitly depend on t. For the positive frequency operator, we have x + ∆ (0) = x + (0) − x + (t + τ ) , but if we substitute this into Eq. (D3), we see that the second term (proportional to x + (t + τ ) ) is exactly zero, so there is no need to find the expectation value of x + at (t + τ ). In practice, we conduct the quantum regression theorem for every t in the last period (in the simulation) of the pseudo-steady-state.
For the more complex second-order correlation function in Eq. (8) of the main text, we have a more complicated version of the quantum regression theorem seen in Eq. (D3) as follows: where we have written U (t) ≡ U (t, 0), we use the notation A(0) = A for an operator A(t), and we make use of the identities U (t + τ ) = U (t + τ, t)U (t) and U (t)U † (t) = U † (t)U (t) = 1 where 1 is the identity matrix. This can be understood simply as the expectation value of the operator x − x + with respect to the effective density matrixρ(t + τ ) = , t), which is the density matrix at t multiplied on the left by x + and on the right by x − and evolved from t to t + τ .  Fig. 2(c) of the main text with (olive solid curve) and without gauge corrections (olive dashed dashed). We also show the BS models, up to first order, again with (blue solid curve, 1/4(1 − 5η/2)) and without gauge corrections (orange dashed curve, 1/4(1 + 3η/2)). The general trends are clearly qualitatively very good, especially at lower η, which is precisely when we expect the BS model to be valid.
For resonant bare dipole and cavity frequencies however, the BS Hamiltonian gives no corrections to the JC energies to order η 2 (and thus often finds more utility in describing the dispersive regime of cavity QED), but does correct the eigenstates. To first order in η, the corrected JC-like states are which in conjunction with the ground state |G ≡ |g, 0 gives the three states considered in the WEA. Note that in this section, we use a notation where the system states are identified by their composition in the transformed frame.
The effect of the counter-rotating terms eliminated in the BS transformation can be quantified by considering the transition matrix elements of the quadrature operator Π which we use to couple to the external reservoir fields. As in the main text, we use Π C to refer to the uncorrected operator in the dipole gauge (which is equivalent in form to the Coulomb gauge operator) corresponding to the electric field quadrature mode operator. Performing the BS transformation, we find Π C = i(a † − a) → P = 1 √ 2 i(a † − a) − η/2σ x , and so x + = −ia − η/2σ − , where x + is the "positive frequency" (taking higher energy states to lower energy ones) component of the transformed operator We introduce the notation P to reiterate that the operator which we assume to couple the system to the reservoir modes is proportional to the momentum quadrature operator of the cavity mode. With gauge corrections, the correct dipole gauge quadrature operator is instead Π D = i(a † − a) + 2ησ x , and so to first order in η, we find P → P = 1 √ 2 i(a † − a) + 3η/2(σ + + σ − ) , and x + → x + GC = −ia + 3η/2σ − . The transition matrix elements (modulus squared) with respect to P are, with no gauge corrections, However, with gauge corrections, we have and we can infer immediately, that the linewidth asymmetry will change with gauge correction. As shown in Fig. 8, the leading order effect in η of the effect of gauge corrections is in excellent agreement with the numerical solution (Fig. 2(c) of main text), for the perturbative regime (η 0.2). For higher values of η, then the numerically exact | G|P|j | 2 with and without gauge corrections explain the main features of the spectra, especially the different linewidths as a function of η, and how these drastically differ with gauge correction. Below we explain why, to the same order of approximations, that the change in linewidth is directly proportional to the weights of the spectral peaks in the spectra.
Solving the relevant Bloch equation with weak excitations, then the spectral linewidths (full widths at half maxima) of the first two excited states are, in the SC limit g/κ 1, simply given by the projections above multiplied by 2κ. This is the primary effect for the observed asymmetry for increasing η (as we can also see from the full numerical calculations). Thus, even in the perturbative BS regime, the asymmetry stemming from the counter rotating wave effects is qualitatively different when one properly accounts for gauge corrections. We justify this assumption in more detail below.
Within the WEA, it is possible to find an analytic expression for the emission spectrum S cav that is valid perturbatively up to order η by solving the above Bloch equations derived from the BS Hamiltonian. In the strong coupling regime, this spectrum takes on a particularly simple form, which is useful to gain qualitative insight into the spectral asymmetries which are shown to arise in our main results. In particular, to leading order in E, the steady state solutions to the Bloch equations give the very simple solution ρ + = ρ − = (1 − ρ G )/2 = E, with all other matrix elements zero.
The steady-state cavity spectrum with incoherent driving is The steady-state correlation function x − x + (τ ) can be calculated with the QRT, and the result for the spectrum after Fourier transforming is where Much simplification can be made if we assume α ≡ g/κ is large, and neglect terms of order 1/α 2 Then, G ≈ g − i(Γ . Without gauge correction, this ratio is ∼ 1 − 3η + O(η 2 ), and with corrections it is ∼ 1+5η +O(η 2 ), which quantifies the change in asymmetry to leading order. We can understand this asymmetry on physical grounds as arising from which polariton branch is more cavity-like: In the BS frame, the BS shift causes a detuning between cavity and TLS resonances, which leads to cavity-like and atom-like polariton branches. In the WEA, both polariton branches become equally populated, and thus, in the SC regime, their spectral weights are determined by the transition matrix elements P G± , or in other words, how much the operator which couples the cavity field to decay channel modes also couples the polariton-ground transition. The more cavity-like transition experiences a larger decay rate, but which branch this corresponds to is dependent on both the gauge corrections and the BS frame transformation. The interplay of these effects thus gives the overall asymmetry.
Finally, to confirm the accuracy of the analytical formula using the same material parameters as in the main text, we show a zoom in of the two main polariton peaks (near ω c ± g) using the full numerical solution versus the simple analytical formula in Fig. 3, using η = 0.05 and η = 0.1. Clearly the comparison is qualitatively excellent and the main differences with gauge corrections stem from the changing linewidth.