Generalized Dicke model and gauge-invariant master equations for two atoms in ultrastrongly-coupled cavity quantum electrodynamics

We study a generalization of the well-known Dicke model, using two dissimilar atoms in the regime of ultrastrongly coupled cavity quantum electrodynamics. Our theory uses gauge invariant master equations, which yields consistent results in either of the standard multipolar and Coulomb gauges, including system-bath interactions for open cavity systems. We first show how a second atom can be treated as a sensor atom to measure the output spectrum from a single atom in the ultrastrong-coupling regime, and compare results with the quantum regression theorem, explaining when they can be different. We then focus on the case where the second atom is also ultrastrongly coupled to the cavity, but with different parameters from those of the first atom, which introduces complex coupling effects and additional resonances and spectral features. In particular, we show multiple resonances in the cavity spectra that are visible off-resonance, which cannot be seen when the second atom is on-resonance with the rest of the system. We also observe clear anti-crossing features particularly pronounced for when the second atom tunes through resonance.

In particular, USC exploits the nature of counterrotating wave physics and pondermotive forces [4,5], and pushes one toward a non-perturbative regime where the light and matter excitations must be treated on an equal footing, i.e., as joint/dressed states [12], where even the ground state can contain virtual photons.These features make the USC regime responsible for many intriguing phenomena including the formation of quasiparticle (e.g., exciton-polariton or plasmon-photon) collective modes with finite lifetimes, as well as hybrid and entangled states with higher degrees of controllablity [8,9,[13][14][15][16][17][18][19].
The intricate interactions between quantized cavity modes and quantum emitters can be modeled in the framework of cavity quantum electrodynamics (cavity-QED), where atoms and atom-like structures (e.g., quantum dots, molecules, superconducting circuits) interact with a (dominant) single quantized cavity mode [20][21][22].Traditionally, strong coupling occurs when the cavityemitter rate, g, exceeds any dissipation rate κ (of the cavity) or γ (decay of the emitter) [23][24][25][26][27][28][29], while the USC regime is characterized not only by the lower rates of decoherence, but also when the atom-cavity coupling strength, g, becomes a significant fraction of the bare energies, ω 0 , of the system, commonly quantified as g > 0.1ω 0 [4,5].Additionally, the hybridization of quantum states with different numbers of excitations in the USC regime results in a population of virtual photons in the collective system's ground state, also with significant dissipation (and thus not even requiring the strong coupling regime) [30].
The profound applicability of cavity-QED and its ease of modeling are derived from truncating the full emitter problem to a two-level system (TLS), which is typically coupled to a single quantized cavity mode.However, the truncation of the Hilbert space, in either the material and/or photonic part, causes problems for gauge invariance when working in the USC regime [12,[31][32][33].Recently, many of these issues have been partly fixed for the standard quantum Rabi model (QRM) Hamiltonian [31,34,35], and extended recently to ensure that dissipation and input/output is also included in a gauge invariant way [34].While we focus our work in the context of natural atoms (including molecules and quantum dots) in cavities, the general ideas also extend to other cavity-QED systems including circuit QED, e.g., see Ref. [36].More general quantization for arbitrary media (dealing with matter and field truncation), and the USC regime, is discussed in Ref. [33].

A. Gauge invariance
In the dipole gauge (specifically, the dipole approximation in the multipolar gauge), the QRM describes the arXiv:2301.02127v2[quant-ph] 31 Jul 2023 TLS-cavity system via the Hamiltonian [12,31] (in units of ℏ = 1): up to a constant (1ω c η 2 ), where ω c is the cavity transition frequency, a (a † ) is the cavity photon annihilation (creation) operator, ω a is the TLS transition frequency, σ z = σ + σ − − σ − σ + and σ x = σ + + σ − , with σ + = |e⟩⟨g| (σ − = |g⟩⟨e|) the atomic raising (lowering) operator; also, g D is the atom-cavity coupling in the dipole gauge (g D ∝ √ ω c ), and η = g D /ω c is the normalized coupling parameter.We can neglect terms proportional to the identity as these do not affect the system dynamics; they simply introduce an offset in the ground state energy, which we can normalize to any value.Equation (1) reduces to the Jaynes-Cummings model (JCM) in the rotating wave approximation (RWA), yielding [37,38] When the system is subjected to matter truncation, ĤD QR produces the correct eigenenergies [39], but the electric field operator [34,36] Ê ∝ −i(a ′ † − a ′ ), where a ′ = a + iησ x , which can be derived from several different viewpoints [31,34,36].For example, in the restricted TLS subspace, one can transform the Coulomb gauge operators to the dipole gauge operators, through the projected unitary transform [31] U = exp[−iη(a + a † )σ x ], so that a ′ → UaU † = a + iησ x [36].These transformed operators then must be used when computing cavity field observables and for deriving master equations.We note that in the Coulomb gauge, the electric field operator can be directly written in terms of the cavity photon operators.However, this is not the case for the electric field operator described in the multipolar gauge, where it takes on a matter component from the TLS [33,34].
In the Coulomb gauge, the standard system Hamiltonian for the QRM is [12,31,40] ĤC,naive 2 /ω a is the ponderomotive coupling strength [41].Unfortunately, this "naive" system Hamiltonian is wrong (which is why we use this name in the superscript) as it does not produce the correct eigenenergies in the USC regime [12], and breaks gauge invariance.The breakdown of gauge invariance here can be seen as a formation of a potential nonlocality due to the truncation of the matter Hilbert space [31,42].Instead, by applying an appropriate unitary gauge transformation to the dipole gauge-independent QRM model, the correct gauge-fixed Coulomb QRM Hamiltonian is [31] which produces identical eigenenergies to ĤD QR .

B. Gauge-invariant generalized master equation
For realistic cavities, one must also account for dissipation/losses and photon input-output channels.Generally, open-system cavity-QED problems are modelled by considering the atom and the cavity are interacting with general baths, as an open quantum system.In such situations, a master equation description is widely used, leading to a detailed understanding of the cavity spectra and other desired observables [43][44][45][46].
Commonly, the bare-state master equation formalism, where the joint basis states are constructed from the bare light states and the bare matter states before light-matter interaction, is often used in open system cavity-QED, yielding the standard Lindblad master equation.However, the bare-state master equation formalism uses the wrong states in the USC regime (including the ground state, which can now be an entangled state of photons and matter) and it has been shown that one needs a dressed-state approach to avoid unphysical transitions [47].Moreover, one also needs a "generalized" master equation (GME) approach to account for frequency-dependent baths and non-secular effects [45].A secular approximation, which relies on, e.g., |ω − ω ′ | ≫ κ for ω ̸ = ω ′ in the GME (described later), is not able to describe dissipation or decoherence in open quantum systems with mixed harmonic-anharmonic or quasi-harmonic spectra [45,48].This is especially important in the USC regime, where g dominates the eigenenergies.Specific examples include cavity QED in the dispersive regime [47], cavity optomechanics [49], resonant Raman interactions in cavity-QED [50], electron-phonon interactions [51,52], as well as frequency-dependent pure dephasing bath [53] and/or radiative decay [54].Thus, our model below with an Ohmic bath (κ(ω) ∝ ω) is also an example application of the GME approach, since a secular approximation leads to the same assumption as a spectrally flat bath (i.e., κ(ω) = κ), and we show explicitly how these differ.
Beyond these details, in the USC regime, such approaches are typically gauge relative, and again one must use a corrected a ′ for cavity mode operators with the dipole gauge or use the corrected Coulomb gauge Hamiltonian.Although such studies have so far assumed very simple models for the system-bath interactions, these approaches do produce gauge-independent results [34,35].
An advantage of using a GME approach is that realistic observables can be computed, such as the cavity-emitted spectra, typically using the quantum regression theorem [34,55].However, gauge-independent GMEs have so far only been applied to the case of one atom/TLS, and we can also expect a significant impact when applied to multiple atoms.In this regard, the Dicke model is a fundamental model of quantum optics describing the light-matter interaction where a cavity mode is coupled to a set of identical TLSs [56,57].The model is known to be an established description for a class of intriguing phenomena in cavity-QED such as superradiant phase transitions and quantum chaos [57][58][59][60][61][62][63][64][65][66][67][68][69].

C. Dicke model in the USC regime
The Dicke model has also been investigated in the USC regime [70][71][72][73][74][75][76][77].In the study of effective light-matter interactions in a circuit QED system, coupled symmetrically to multiple superconducting qubits, Ref. [71] studied a microscopic model Hamiltonian that not only describes the usual collective qubit-photon coupling but also the effect of direct qubit-qubit interactions.Various other works in the USC regime have been presented on related coupling effects, mainly at the thermodynamic limit, using simple system Hamiltonians [78,79].In these extended Dicke models, similar to the previous studies on the Dicke model or even the Hopfield model in the USC regime [72][73][74][75], the atoms are degenerate (i.e., they share the same coupling coefficient and resonant frequency).Recent studies also include gauge-invariant system Hamiltonian models [75], or discuss more exotic schemes of the Dicke model, such as the anisotropic or nonequilibrium models in which the counter-and corotating terms have different coupling strengths, but the two atoms are still identical [76,77].

D. Generalized Dicke model
It is desirable to explore a more general two-atom case where the TLS parameters can be different, and a natural extension to investigate is a system of two dissimilar atoms, which we term a generalized Dicke model (GDM), in the limit of two atoms.Extensions to multiple dissimilar atoms could be the subject of future work.From a practical viewpoint, one must also include realistic dissipation and input-output channels to the system.In this paper, we present such a study, using gauge-independent master equations valid for exploring USC dynamics.Ultimately, this GDM is a more realistic scenario for studying how atoms interact in the USC regime, as it is practically impossible to experimentally produce two identical TLSs for experimental systems [80][81][82][83].Coupling with two different atoms also leads to new coupling regimes that are not accessible with identical atoms.
With the recent technological and experimental advances, promising opportunities are emerging to study GDM [1,8,15,80,[84][85][86][87][88][89][90].While most of these are in the strong coupling regime, new designs and systems are emerging in the USC regime as well.References [4,5] give a review of the various experiments systems, including molecular as well as solid state systems.Moreover, it has been extensively shown that the efficiency of lightmatter interaction in bound systems can be enhanced in nanostructures, especially using metallic nanostruc-tures, where experiments have even demonstrated the deep USC regime (g > 1) [6].
We study the two-atom GDM by introducing a dissimilar second atom to a general one-atom-cavity USC problem, using a gauge-invariant GME description.We exploit this model in two different ways: (i) we first introduce a second TLS as a weakly coupled sensor atom for the cavity-emitted spectrum [sensor atom approach, Fig. 1(a)], and show that it produces qualitatively similar spectra to that computed with the quantum regression theorem, though only with certain types of bath coupling; we also confirm that these sensing atom results are identical in both the dipole gauge and Coulomb gauge, as they must be; (ii) we then focus on the main topic where the second atom is now also treated as an ultrastrongly coupled atom, distinct from the first atom [GDM, Fig. 1(b)], and demonstrate several new spectral features that emerge as we change the coupling parameters of the second TLS.
The rest of our paper is organized as follows: In Sec.II, we present the main theory, which includes a description of the GME, our excitation scheme, as well the various system Hamiltonians, bath interactions, and observables, including the cavity-emitted spectra.
In Sec.III, we present the main calculations and results for the sensing atom approach, and show how the sensing atom coupling can be used to model the detection of light.We also show how these results compare to calculations with the quantum regression theorem and explore the more general case of different bath couplings (for the atoms as well as the cavity).Next, in Sec.IV, we consider the case of two atoms in the USC regime, where we change the parameters of the second atom, and study the effect that this has on both the system eigenenergies as well as the cavity observables.We first show explicitly how our GME produces gauge-independent results when using the correct gauge-fixed approaches as described in the main text.Subsequently, we then present a series of investigations using the dipole gauge.Finally, we conclude in Sec.V.

II. THEORY
In this section, we present the GME, as well as the different bath models and system Hamiltonians that we will use.We also show how these can be used to compute the cavity spectra, using either the quantum regression theorem, or a sensing atom approach.

A. Generalized master equation
We first introduce the main GME that we use to compute the key observables of interest: where ρ is the composite system (composed of the cavity and the atom, or atoms) density matrix, and Ĥ ≡ ĤD/C is the system Hamiltonian in either gauge (dipole, 'D', or Coulomb, 'C').
The Lindbladian for each dissipation channel is of the same form, so we write it generally as [45] L Since we now have several possible dissipation channels for the cavity, and the atoms, Λ indexes the cavity and the atom, or atoms.The dressed operators are defined from with ] † , and we assume that Π has electric field coupling, such that ΠC = i(a † −a) in the Coulomb gauge, and ΠD = i(a ′ † −a ′ ) in the dipole gauge [34].We note that the dressed eigenstates {|j⟩} are required to construct the correct dressed operators utilized in the GME; these are the eigenstates of the full light-matter system Hamiltonian including the interaction term [34,45,54].The dressed states are naturally gauge-dependent, but the observables are not.
Modeled by a (continuous) superposition of damped bosonic harmonic oscillators, baths are generally described by their correlation functions and, in turn, their spectral densities of states which contain information on the frequencies of the baths' modes and their coupling to the system [45].For our purpose, the frequencydependence of the baths is modeled as either a flat bath, or an Ohmic bath, However, in the case of a sensor atom, we use since in reality the sensor will also have a center frequency at the main detection frequency of interest, while we assume is at the cavity resonance frequency.If the open system also includes a sensing element, special considerations for the sensing atom's bath must be taken into account.Essentially, we must add the dissipation channel for this sensor atom in an analogous way to the primary atom.However, in principle, we require that the inclusion of the sensor should act as a noninvasive measurement.Therefore, we must ensure that γ s ≪ κ in either the flat or Ohmic shape of Γ sen (ω), where γ s is the sensor atom decay rate, or else the sensor atom introduces additional broadening to the existing peaks in the spectra.Careful attention is also needed as the dissipation rate of the sensor puts a limit on the coupling strength between itself and the cavity.We will cautiously take into account these considerations in our results.
For a cavity-QED system in the USC regime, the γ ≪ κ process is usually negligible; however, γ plays an important role in the sensing atom approach (for its light detection), so we keep the bath functions general for such a study.However, in the case of two atoms in the USC regime, we will use Ohmic baths throughout, where only Γ cav (ω) is generally important.
For the excitation process, we also include the incoherent driving through the pump Lindbladian, with where , and P inc is the incoherent driving strength.

B. Observables
Now that our main master equation model is established, we next present the key observables with which to explore the dynamics of the system.These can also be used to ensure we have properly enforced gauge invariance.We will focus on the cavity-emitted spectrum.
The cavity spectrum is typically computed from the Fourier transform of the two-time cavity correlation function, which exploits the quantum regression theorem.In such an approach, the cavity spectrum is defined from [91] S cav (ω) where ω is the emission frequency.With incoherent steady-state driving, this simplifies to a single time integration, carried out after the system dynamics has reached steady state.
An alternative method for computing the spectra is to include a sensing atom, and compute its excitation flux.Reference [92] showed how normal-order correlation functions, used to compute the spectrum and other observables, can be computed from a frequency-tunable sensing atom in the limit of small coupling with the field.In the USC regime, such methods have been discussed at the system Hamiltonian level [36], and here we test how well such an approach recovers the same sort of spectra as the quantum regression theorem.If such a model is correct and is gauge invariant, it naturally extends to allow us to explore two atoms in the USC regime.In the latter case, we will use the quantum regression theorem again, primarily for convenience.However, we remark that the sensing atom approach has several potential advantages: (i) it does not require a Born-Markov approximation to be valid, and (ii) it can easily be used to model pulsed excitations, without the need for a double-time integral.

C. Spectrum detected from the sensing atom approach
The detected spectrum from the sensing atom approach (SAA) is defined through where ω s is the sensing atom frequency and, according to Eq. ( 7), X + sen = ⟨j|σ x,s |k⟩ |j⟩ ⟨k|, with σ x,s being the Pauli x-matrix for the sensor.In contrast to the quantum regression theorem, such an approach does not require any two-time correlation functions.Moreover, since we are exciting the system with a steady-state drive, then which is (again) computed when the system reaches a steady state.This method has a simple physical interpretation; the sensing atom excitation number is proportional to the photon flux of the cavity-QED system which contains another atom in the USC regime.The sensing atom then "detects" the cavity output flux.
For the SAA to be valid, the sensing parameters should generally be noninvasive to avoid affecting the detected spectrum.To be specific, the sensor atom should have a vanishing coupling strength, g s ≪ g.In order to determine acceptable parameters for the "sensor atom" (non-perturbative coupling), we use parameters that provide constant results over a range of frequencies, with acceptable run times, which are guided by the criteria g s ≪ γ s R/2, where R is any rate in the system [92].

D. Photodetection rate of cavity photons
Another useful quantity to calculate is the photodetection rate of cavity photons, emitted from the |j⟩ → |k⟩ transition [54], which is proportional to the P ( Π/ √ 2) quadrature matrix elements squared, namely: This is the main system-level quantity that affects the spectral transition rates; however, the transition rates, T ij , also must be multiplied by a factor D 2 (ω jk ), where D(ω) is the density of states of the relevant bath [41], so that, in the case of cavity emission, T jk = 2πD 2 cav (ω jk )|P jk | 2 ; this can be derived from Fermi's golden rule.Moreover, in the presence of a sensing atom, there is additional filtering through the sensing atom's density of states, as will show below.
For all our numerical calculations, we will use Python and we also exploit the QuTiP module for quantum objects and operations [93,94].

E. System Hamiltonians and gauge-fixing for the sensor atom approach
We first introduce a second atom (TLS) as a sensor for the cavity-emitted spectrum.While the sensor atom approach does not introduce any further observables to probe or any vastly new physics, it provides a check for gauge invariance and verifies that the second atom is included correctly in the model, before elevating it to a second atom also in the USC regime.It also demonstrates the influence of an additional bath coupling, which is relevant for other types of detection, including cavity detection.As mentioned earlier, this method for simulating spectra also holds some potential advantages over the quantum regression theorem, which requires the calculation of the two-time correlation function.In particular, the sensing atom approach is potentially more powerful when used to compute various multi-time correlation functions [92] and to model short-pulse excitation.The sensor atom approach is also a valid physical model for the detection of photons emitted from the cavity, and a similar approach could be used to describe a sensing cavity as well.
In order to not affect the spectrum, this sensor atom should have a vanishing coupling strength, g s ≪ g.In practice, however, it may have a minor influence on the computed spectrum, and a qualitatively different one if it also has a non-trivial bath function.
The system Hamiltonian in such a case may be naively constructed by the addition of two terms: (ω s /2) σ z,s and ig D s (a † − a)σ x,s , namely the bare Hamiltonian of the sensor, and the interaction Hamiltonian between the sensor and the cavity, respectively [92], to Eq. ( 1); this is analogous to that of the primary atom, in the dipole gauge.Perhaps counter-intuitively, after our discussions on the dipole gauge Hamiltonian, here the gauge correction, including only the main atom, also needs to be applied at the Hamiltonian level for the sensor interaction [36].This is because the sensor atom couples to the electric field of the cavity with its coupling to the primary atom already included, which explicitly contains the corrected â′ operators (similar to the cavity bath operators terms in the GME).Therefore, the naive choice is incorrect, and gauge fixing must be applied so that the interaction Hamiltonian between the sensor atom and the cavity becomes since η s ≪ η and also σ 2 x,s = 1 (which only gives an energy offset).
Thus, in the dipole gauge, the gauge-corrected full system Hamiltonian reads [36,54,95] Applying a RWA, then we have where η is the normalized coupling for the primary atom.
Clearly, in the sensor atom approach, the sensor atom does not need to modify the principal atom-cavity coupling and related observable operators, though its bath interactions can play a qualitatively important role.In the USC regime, of course, the RWA does not work, but it is useful to highlight the effects of counter-RWA terms, at least at the level of how these affect the system eigenfrequencies.
In the Coulomb gauge, the Hamiltonian can be also obtained from the dipole-gauge one, by conducting the unitary transform ĤC where , and again we have neglected terms proportional to the identity operator.Note here that, as in the Coulomb Hamiltonian without the sensor, there is no separable bare sensor Hamiltonian when including the gauge correction.However, without the gauge correction, we have, (ω s /2)σ z,s +g D s (a+a † )σ y,s for the bare sensor and cavitysensor interaction Hamiltonian.
Naturally, we can apply this transformation in reverse, i.e., to the cavity operators in the Coulomb gauge, to find the two-atom version of the gauge correction in the dipole gauge.As expected, this results in the corrected operators in the dipole gauge, and yields As discussed before, computing field observables in the dipole gauge must then be done with these corrected electric-field operators.The GME for the sensor atom approach has the same format as Eq. ( 5).However, one must use the correct definition of Hamiltonian Ĥ ≡ ĤD/C SAA , to ensure gauge-invariant results, and also one has additional bath coupling terms, which need not have the same spectral function (i.e., the density of states seen by the sensing atom could be different to the density of states seen by the cavity).

F. System Hamiltonians and gauge-fixing for the generalized Dicke model
The Dicke model in the USC regime takes the QRM and adds a second identical atom to the system, also in the USC regime.This also has to be done with the gauge-corrections, to ensure the results are gauge invariant.The difference between the regular (or previously extended) Dicke model [75,[96][97][98] and the GDM is that we allow the two ultrastrongly coupled atoms to vary in their frequency and coupling strength.We have already described how to include a second atom (TLS) into our system in both the dipole gauge and the Coulomb gauge, using the sensor atom approach, where gauge invariance is also ensured.The only difference here is that we use Eq. ( 12) to compute the spectrum since the second atom is now participating in energy exchange with the cavity and its population cannot be relied on to obtain the spectrum, as it is no longer acting as a weak sensor.Of course, one could bring in a third atom as a sensor, but in general, the quantum regression theorem is efficient and accurate for the problems we study below, especially with time-independent incoherent driving.
Conveniently, the required Hamiltonians for two atoms in the USC regime are similar to the above, except that for the operators and quantities associated with the original atom we assign a subscript 'a' and for those of the sensor atom we assign a subscript 'b', and we no longer assume η b ≪ η a .Thus we have and for reference, with the application of the RWA, one has then the correct Coulomb-gauge Hamiltonian is (24) Simulations can then proceed as before, e.g., if using the dipole gauge, then one must use the corrected cavity operator to compute cavity observables in the dipole gauge, and these are also used for deriving the dissipation terms in the GME.The GME [Eq.( 5)] then uses the appropriate Hamiltonians, Ĥ ≡ ĤD/C GDM , which are the total gaugecorrected Hamiltonians of the system in their respective gauge for the GDM.In the Coulomb gauge, no modification is needed on the field operators.

III. RESULTS AND DISCUSSIONS OF THE SENSOR ATOM APPROACH
Thus far, we have presented a gauge-corrected model for the system Hamiltonian for the sensing atom approach; namely, one that yields the same eigenfrequencies for either the dipole gauge or the Coulomb gauge.This leads to the proper understanding of the eigenenergies and eigenstates of the closed system, which is basically an extension of the original QRM.We must also consider dissipation for this sensor atom, including it in an analogous way to the primary atom, under the condition γ s ≪ κ.The chosen dissipation rate of the sensor also puts a limit on the coupling strength between itself and the cavity.In general, we require that the coupling must be small enough to ensure that losses from the cavity into the sensor and the sensor back-action into the main system are negligible.This crucial argument leads to an acceptable range of parameters discussed in Ref. [92], and also discussed in more details below.
For computing spectra in either considered gauge (dipole and Coulomb), we then allow the system to evolve to steady state, also including an incoherent pump term.Once the steady state has been reached, we take the expectation value of the sensor excitation.We do this for a range of frequencies of interest to form the sensor atom spectra, performing a calculation for each scanned ω s .
For this sensor atom section, we will focus our attention on the sensing atom interacting with the primary atom in the USC regime, using a fixed coupling parameter of η = 0.5 to the primary atom.This USC regime has been studied recently using gauge-independent master equations, and shown to yield identical results in the dipole and Coulomb gauges, and is thus an excellent testbed to also compare with a sensing atom simulation [34].
We must first ensure an approximately vanishing coupling rate compared to the coupling between the cavity and the main atom.We take g s = 0.001g to satisfy this condition.Then, to obtain a lower limit on γ s , we find the smallest transition rate in our system at η = 0.5 to be R ≈ 0.3g, and we chose a slightly smaller value than this, to satisfy the previously-mentioned criterion g s ≪ γ s R/2 [92].If we then use γ s = 0.0025g, we obtain γ s R/2 ≈ 0.02 ≫ 0.001g = g s .Therefore, we use κ ≫ γ s ≥ 0.0025g, as the acceptable range of values to choose from.

A. Dressed eigenenergies/eigenstates and example transitions
In Fig. 2(a), we plot the eigenenergies of the singleatom cavity-QED without (blue solid curves) and with (red dashed curves) the RWA.This helps to highlight the role of the counter-rotating wave terms for increasing η.The computed energies are gauge-independent (namely, the Coulomb and dipole gauge results yield identical results), as they should be.In the sensor atom approach, one expects no difference between the main eigenenergies with a single TLS-cavity system, which we have confirmed to be the case; however, additional states naturally appear because of the sensing atom states, which depend on ω s .
The three significant optical transitions are identified with the downward arrows and the letters 'A' (ω 10 /ω c ≈ 0.5), 'B' (ω 31 /ω c ≈ 8.2), and 'C' (ω 20 /ω c ≈ 1.45), for η = 0.5 (primary atom).These transitions are responsible for the significant peaks that will appear in the incoherent spectra shown in Fig. 3, discussed below.1) and ( 2), respectively, and (b) the sensor atom approach Hamiltonian, where both the primary and the sensor atoms are in resonance with the cavity, ωa = ωs = ωc, and the sensor atom has the coupling of g D s = 0.001g D ; and (c) the sensor atom approach Hamiltonian, where the primary atom is in resonance with the cavity, ωa = ωc, and the sensor atom has ωs = ωc/2 and g D s = 0.001g D , obtained from Eqs. ( 18) and ( 19), respectively.In panel (a), the three significant transitions are identified with the downward arrows and the letters 'A' (ω10/ωc ≈ 0.5), 'B' (ω31/ωc ≈ 8.2), and 'C' (ω20/ωc ≈ 1.45), for η = 0.5.These resonances are highlighted for reference when explaining the key features of the cavity spectrum.
When the primary atom and sensor atom are both on resonance with the cavity, as expected, the eigenenergy lines start together, at low η, and then split from the same initial points at multiples of the cavity transition energy, as shown in Fig. 2(b).When the primary atom is on resonance but ω b = 0.5ω c in Fig. 2(c), the addition of extra eigenenergy lines starts at half-multiples of the cavity transition energy.In all three panels of Fig. 2, the deviation of the JCM eigenenergies from the full QRM Hamiltonian eigenenergies is apparent as the normalized coupling parameter increases, a failure that is of course fully expected [4,5].

B. Cavity-emitted spectra via incoherent driving
We next study the cavity-emitted spectra via incoherent driving, and consider weak drives so as not to perturb the system eigenstates too much.We will compare our sensor results to those obtained using the quantum regression theorem [34].Additionally, we investigate the effects of changing the various bath spectral functions.
In the top two rows of Fig. 3, panels (a-d), we consider a flat bath for the cavity (i.e., Γ cav (ω) = κ) and show the effect of changing the atomic baths from flat to Ohmic.Using the quantum regression theorem, changing the atomic bath has almost no visible effect on the spectrum when the cavity bath is flat.However, the spectrum detected by the sensor atom is drastically modified if the atomic bath of the sensor atom has a non-trivial frequency dependence.This can be viewed as an additional filtering process, e.g., in the case of an Ohmic sensor bath (i.e., Γ sen (ω) = γ s ω/ω c ), there is increasing dissipation at higher frequencies, thus reducing the strength of the peak on the right and increasing the relative strength of the peak on the left.
Next, as shown in the bottom half of Fig. 3, panels (eh), we consider again an Ohmic cavity bath and look at the effect of changing the atomic baths.The first result to note is that the Ohmic cavity bath produces the largest change of any of the models explored here.This is not too surprising, as the cavity dissipation is the largest by far to begin with (κ = 0.25g vs. γ = γ s = 0.005g), and we are modeling cavity emission.Thus, the dissipation is overwhelmingly dominated by κ, and any frequency dependence included with it will have a larger effect.Also, since the two models (quantum regression theorem and sensor atom) here have the same dependence on the single cavity, we see the change in the cavity bath having a similar effect on both spectra (specifically, reversing the asymmetry and modifying the relative peak heights to a similar extent).When we now also change the atomic baths to be Ohmic, we see similar effects to those above.The quantum regression theorem results are now slightly affected, and we again see a large effect on the sensor approach model.
Obviously, the gauge-uncorrected results are not expected to produce the correct physical results in the USC Cavity-QED incoherent spectra for one atom in the ultrastrong coupling regime.Cavity emitted spectra were computed using the sensor atom approach (left column) and the full quantum regression theorem (right column).On the right, we list the type of bath used for the cavity and the two atoms (which use the same bath type).Gauge-corrected (not corrected) results are shown with solid blue (dashed orange) curves.We use incoherent driving with Pinc = 0.01g.Other system parameters are κ = 0.25g, γ = γs = 0.005g, gs = 0.001g, ωc = ωa, η = 0.5.The three significant peaks A, B, and C represent the major transitions shown in Fig. 2.
regime, and thus one obtains different plots when computed by different approaches.Beyond this, as shown and discussed above, the gauge-fixed results, with the two different detection models, may produce different results depending on the model for the bath function of the detection.Consequently, the various bath functions are generally important in the USC regime (e.g., the cavity bath may alter the transition rates by its spectral density).Depending on the context, one may argue that in certain detection models, the results from the SAA would be more experimentally relevant.In terms of highlighting the main physics for our two-atom USC below, either model is adequate, and we will use the quantum regression theorem approach as it is simpler and more computationally efficient.

IV. RESULTS AND DISCUSSIONS OF THE GENERALIZED DICKE MODEL
In the description of the GDM, we extended our sensor atom approach to be a primarily part of the coupled system (i.e., no longer weakly coupled but also in the USC regime).Using this approach, we now allow the second atom's properties to vary relative to the first atom, but now when the second atom is also in the USC regime.Our two-atom Hamiltonian in the dipole gauge [Eq.(22)] is equivalent to the extended Dicke model in Ref. [71] in the case that the two atoms are degenerate (i.e., g a = g b and ω a = ω b ).However, our main focus will be on analyzing spectra obtained with dissimilar atoms, where the coupling parameters and resonant frequencies need not be the same.
Similar to the sensor atom approach, one has to first identify the dressed operators which are now found using the eigenstates of the full gauge-corrected GDM Hamiltonians, including the second atom.We can also consider dissipation for this second atom, including it in the same way as for the primary atom, but these rates are basically negligible, as cavity decay is the main source of loss.In either the Coulomb gauge or the dipole gauge, we then allow the system to evolve to a steady state, again including an incoherent pump term.From now on, we use the quantum regression theorem to compute the spectra.Our first calculations will show explicitly the effect of gauge fixing and confirm gauge invariance for the spectra, and then we just choose the dipole gauge, since both gauge results yield identical results.

A. System characterization: dressed eigenenergies/eigenstates and transitions
In Fig. 4, we show the first seven eigenenergies of the GDM without (blue solid curves) and with (red dashed curves) the RWA.In panel (a), we display the eigenenergies with respect to the variation of the equal normalized coupling parameter of the two atoms, η ≡ η a = η b .In the GDM, one expects a considerable difference between the resulting eigenenergies compared to the single TLS-cavity system (or with the sensor atom approach) as compared to Fig. 2. In particular, one observes significant hybridization of the two TLSs in the system leading to the production of the splitting of the eigenstate curves.As expected from our previous observation in the sensor atom approach, when ω b = 0.5ω c in Fig. 4(a), the extra eigenenergy lines start at half-multiples of the cavity transition energy.As opposed to a regular Dicke model with similar atoms, the different starting and splitting point of eigenenergies results a different hybridization and crossing/anti-crossing.Therefore, the possibility of the existence of more exotic transitions and spectra, in comparison to the one-atom spectra and regular Dicke model, appears.
In panel (b) of Fig. 4, the GDM eigenenergies are plot-  ted when the coupling parameters have a phase difference.Letting g b = g a exp[iπϕ], we equate their amplitude but vary their phase via the sweep of ϕ from 0 to 1.The eigenenergies in both models show symmetry and crossing at ϕ = 0.5, as the sign of the real part of the coupling does not change the physics.In Fig. 5, we show further details about the eigenstate properties and transitions.In particular, in Fig. 5(a is the total excitation number (in the dipole gauge).We label the states in Fig. 5(a) even (odd) if their parity is positive (negative).We see that for the considered range of frequencies, the first three excited states have odd parity while the ground state and the 4 − 6 th excited states have even parity.Correspondingly, in panel (b), we plot the energy eigenvalues of these lowest seven states, where we distinguish their parity and label the main transitions that will show up in the cavity spectra.Also, in Fig. 5(c), |P ij | 2 for these transitions are shown which are related to the transition rates (which also depend on the density of states of the cavity bath).All of these properties clearly explain the main spectral peaks that emerge in the computed spectrum.
As mentioned previously, and defined in Eq ( 16), we relate the rate of a transition from state |j⟩ to state |k⟩ [34,41], as proportional to , in the dipole gauge.

B. Cavity emitted spectra via incoherent driving
Our first set of GDM spectra to study is with ω a ̸ = ω b , assuming that the two atoms have the same coupling strength.We will first highlight that our current models do indeed ensure gauge invariance.
In Fig. 6(a), we compare the spectra obtained in the dipole gauge and Coulomb gauge and also show the naive non-gauge-corrected counterparts.Throughout this section, we use Ohmic baths for the cavity and atoms.We display the spectra as a function of the second atom's frequency, while the first is held on resonance with the cavity mode, and both atoms are in the USC regime.It can easily be seen that, while the non-gauge-corrected spectra do pick up some of the correct features, they clearly do not satisfy gauge invariance.The corrected spectra are not only clearly gauge invariant but are also much richer, with additional features including a visible anti-crossing around ω b /ω c ≈ 1 and the disappearance of a main peak in this regime as well.
Beyond the numerical success of our model and codes, we can once again identify key differences in the behavior of our system with and without gauge correction.In all subsequent calculations, we will just choose the dipole gauge, since the results produce observables that are clearly gauge-invariant.
In Fig. 6(b), we plot selected spectra at a few selected ω b values of interest, along with the spectra in the absence of the second TLS.The anticrossing behavior in Fig. 6(a) is at its closest at ω b /ω c ≈ 1; examining the 2D spectra at this frequency in Fig. 6(b), we can see that the splitting is about 0.063ω c = 0.126, or about g/8.The location of this minimal splitting can be understood by looking at Fig. 5(b) and noting that the eigenvalue of the third excited state is closest to twice that of the first excited state at this point, thereby making the frequencies of the B and A transitions closest.
The change in |P ij | 2 (which is proportional to the transition rates) is shown in Fig. 5(c), which correlates with a change in the associated peak's height, e.g., peak C, which is absent from the spectra (peak height is zero) on resonance; this can be explained by the transition rate going to zero in this regime, as shown in Fig. 5(c).This is a strong indicator that there are features of the GDM that can only be accessed when the atoms are dissimilar (i.e., when the second atom is off-resonance).However, one cannot solely rely on the P quadrature matrix element squared to determine the relative heights of the peaks.If this were the case, we would expect peak G to be the largest by far, and peak A to be very small, whereas it dominates above ω b /ω c = 0.5.This is primarily due to the increased damping of higher states [99], and the fact that, as mentioned before, one must take into account the effect of the cavity Ohmic bath in the definition of the transition rate as T jk ∝ ω 2 ij |P jk | 2 .Indeed, this clarifies that transitions involving higher states (D, E, F, G) are significantly more broadened than those involving just the lower states (A, B, C), which appear very sharp on the spectra.Furthermore, transition C has a lower |P jk | 2 value than B, but since B involves higher states, C dominates.Even peak A, with a far lower |P jk | 2 value than B, is larger due to the effect of damping.To sum up, this trend largely depends on the spectral bath function so that one can expect more broadening at higher energy levels, as well as the proper definition of the transition rate for a general (though Ohmic here) relevant bath.
From the eigenvalues in Fig. 5(b), we can identify the peaks with their transitions.As labelled in Fig. 5(c), the visible peaks when η = 0.5 are caused by transitions varying from |1⟩ → |0⟩ to |6⟩ → |0⟩.Hence, it is clear that the incoherent drive (even though weak) excites states up to at least |6⟩.Moreover, it appears that most of the peaks are due to relaxation to the ground state.This is partly because higher-order photons are already part of the lower hybrid states in the USC regime.The key transitions are summarized in Table I.Now that these peaks have been identified with specific transitions, below we next vary the coupling strength of the second atom to determine how the spectra are affected.

Peak Transition
Table I.Identification of the key transitions causing some of the peaks.Note that not all peaks are visible at all values of η.At η = 1, there are other peaks present that we have not labeled.
Next, we vary η from near the threshold of the USC regime (η = 0.1), to the verge of the deep strong coupling regime (η = 1), plotted in Fig. 7.At η = 0.1, we see a reasonable level of symmetry around ω = ω c , yet we also see the appearance of a new resonance which anticrosses with the lower polariton peak, near ω/ω c ≈ 0.9.As we increase η, this symmetry significantly reduces and the anticrossing peaks shift to lower frequencies, while the general Rabi splittings increase as expected, in addition to various Stark shifts.We also see reduced broadening (sharper peaks) with increasing η for the lower frequency peaks, as expected from the GME baths.At η = 1, we discern some of the background peaks becoming the main peaks, and the apparent anti-crossing at lower η appears to become a true crossing, i.e., near ω/ω c ≈ 1.One of the peaks we can identify through the entire range is the one that appears forbidden (or highly reduced) when the second atom is near resonance.This peak can be identified as the C (|2⟩ → |0⟩) transition.However, states |2⟩ and |3⟩ cross in energy between η = 0.5 and η = 1, and are degenerate up to about ω b /ω c = 0.5 at η = 1.For simplicity, we retain the label |2⟩ even after it crosses with |3⟩, so that this feature is indeed due to the same transition throughout.

Influence of amplitude variation of g b
In the above investigations, we chose a few values of ω b to study in detail.We now extend this study further, by examining the role of g b when it varies from zero to g a .First, in Fig. 8, we show how the spectra change when increasing the second atom's coupling strength at a few interesting values of the second atom's frequency.The main feature we can identify is a splitting of some peaks with increased g b and, interestingly, the merging of some other peaks.Some peaks also shift in frequency without any other behavior appearing (Stark shifts).
In the first example, at ω b /ω c = 0.5 [panel (a) of Fig. 8], we see one peak splitting into three at low frequency.Since we have already identified the origin of these peaks and given them labels at g b = g a , we can easily explain where this splitting comes from by examining the change in energy eigenvalues as we increase g b .In Fig. 9, we can see that states |1⟩ and |2⟩, initially near degenerate at g b = 0, split in energy.Recalling from Table I that peaks A, C, and B are due to transitions |1⟩ → |0⟩, |2⟩ → |0⟩, and |3⟩ → |1⟩, respectively, we can see why these peaks decrease, increase, or remain roughly unchanged in energy respectively over the range of g b considered here.Turning now to the peaks involving higher states, these are more complex due to the anticrossing of states |4⟩ and |6⟩ (labeled according to the order at g b = g a , to be consistent with the previous sections) around g b /g a = 0.8.Considering g b = g a , the energy differences in Fig. 9 do explain the peaks with the same transitions as in Table I.Below the anticrossing, however, the peaks can only be explained by different transitions, namely switching |4⟩ with |6⟩.
Next, consider the case of the resonant second atom [panel (b) of Fig. 8].Once again the bright left-most peak can be trivially associated with the A transition.Similarly, the next peak, which almost merges with the first, is identified as transition B, as expected.Transition C is not visible in this regime, but transition D is seen as a broad peak at high g b .Transition E is also not visible, but transition F is visible throughout and transition G is visible at high g b .
Finally, at ω b /ω c = 1.5 [panel (c) of Fig. 8], we again see a significant dressing of the resonances as we change g b , and all of the peaks can be identified as aligning with the transitions in Table I throughout.The differences here are that transition C is strongly visible throughout and merges with F at low g b and that peaks D, E, and G are not visible, except D at higher g b .

Influence of phase variation of g b
The transition dipole of the two TLSs might not be necessarily in the same direction, e.g., if the atomic dipoles are anisotropic and/or the field-polarization is different at the different atom locations.This will change the nature of the couplings in our GDM from pure real to generally complex quantities.Such effects have implications in real-world nanoengineered photonic systems, to manipulate the quantum states and control quantum optical interference effects [100].
To investigate the effects of a phase-dependent GDM, we next allow g b to be complex, and vary its phase.We take g b = g a exp[iπϕ] and sweep ϕ from 0 to 1, similar to Fig. 4(c).In Fig. 10, we show the spectra for g b ranging from g a to ig a to −g a .Apart from being gauge independent for all results, we mention that all three of the 3D   spectra (contours) in Fig. 10 can be simulated in a matter of minutes on a standard desktop computer, where we use typically 200 bare photon states and 12 dressed states.Moreover, a single 2D spectra can be calculated at arbitrary coupling strengths, including complex coupling from the second atom, in typically a few 10s of seconds.Thus the dressed-state truncation is not only necessary for the GME, but it considerably simplifies the numerical Hilbert space from a bare state basis.
As observed above in Fig. 6(a), the gauge-fixed model produces peaks that are absent without the gauge correction, or vice versa.Indeed, at ω b /ω c = 1.5, we can identify peak C as persisting throughout the change in ϕ.Conversely, at ω b /ω c = 0.5, peak C is completely absent at ϕ = 0.5 (but, persists throughout the range without the gauge-correction [54], not shown here).Extending the second atom's coupling strength to a complex quantity increases the separation between the first three peaks The broadening is quite low (comparatively) throughout, but the phase change does act to increase the strength of the A transition at ω b /ω c = 0.5, and increases the broadening of peak C at ω b /ω c = 1.5.We finally note that the spectra are symmetric about ϕ = 0.5, meaning that the spectra are invariant to changes in the sign of the real part of the coupling strength.All these features are in accordance with the eigenenergy lines in Fig. 4(c).

V. CONCLUSIONS
We have presented a gauge-invariant GME approach to model two atoms in ultrastrong coupling regimes of open system cavity-QED, where the atoms are modelled as Fermionic TLSs.We first analyzed the applicability of a sensor atom approach for computing the detected spectra from the cavity-QED system.This is an alternative approach to using the quantum regression theorem, allowing for the computation of spectra even when driving with extremely short pulses or with multiple timedependent fields, when the spectra solution from the quantum regression theorem may break down.This twoatom model also provides confirmation of the gauge independence of the general theory for light-matter interaction in the USC regime [34], when more than one atom is included in the system, which is a non-trivial task, even without including dissipation and optical excitations.
Using incoherent driving, we demonstrated the ability of the sensor approach to produce spectra that match well with the quantum regression theorem results, when using spectrally flat baths.We also showed the influence on the spectra when changing the bath function for both the cavity and atomic baths.We compared the Ohmic and flat baths for each case and demonstrated that the spectra only agree well when the atomic baths are flat.This is, however, not a realistic model for real-world detection over very large frequencies.
For the main part of the article, we then presented results obtained using a generalized Dicke model, in the limit of two atoms.Previous studies on the Dicke model have used identical atoms, only varying properties of both at the same.However, it is practically impossible to produce this situation in a physical lab environment.Motivated by this fact, our studies presented results obtained with dissimilar atoms, extending previous works.We first showed that our model produces gauge invariant results when including the gauge correction terms.We also showed that the gauge corrected spectra (correct spectra) are much richer than naive models and with more striking features.
We then examined the effect of allowing the resonant frequency of the second atom to vary, and showed that there are significant peaks visible off-resonance that cannot be seen when the second atom is on-resonance with the rest of the system.We demonstrated that this effect holds for a large range of normalized coupling strengths even down to the verge of USC.This shows that this first extension, namely the ability to model two atoms with dissimilar resonant frequencies, has important implications not just in the usual USC regime.We also identified the main transitions for these visible spectral peaks.
Next, we chose a few values of the second atom's frequency, including resonant with the first atom, to explore the second extension of our model, where we changed the coupling strength of the second atom relative to the first.We observed that some of the separate peaks can only be identified as separate peaks due to the coupling of the second atom.Indeed, a single peak without the second atom's coupling splits into three when the coupling is introduced in one of the regimes considered.Finally, we allowed the second atom's coupling to have a phase difference relative to the first, and showed how the relative phase can substantially tune the spectral energy levels.
Qualitatively, the degree of tunability of the system to shift, produce or nullify the resonances in the emission spectra is much richer in a GDM.A clear observable to probe is the nullification of the main cavity resonant radiation when one of the dissimilar atoms is off-resonant with the cavity, as well as new resonances unique to the USC regime.With the naturally entangled states (even the ground state) in the USC regime, when a second dissimilar atom is added to a single-atom cavity system, the transitions between the states are highly modified so that some of the natural cavity-QED radiation modes are absent, and/or higher-order photons are already part of the lower hybrid states in that regime.In addition, the effects of the separate atomic baths may be intensified in a GDM where, for example, they can massively broaden the higher-order photons.This can relate to emerging experimental systems for probing the USC regime in the near future [4,5].

Figure 1 .
Figure 1.Cavity-QED schemes with two atoms.Schematics of the cavity-QED model with a second atom, including: (a) the sensor atom approach and (b) the generalized Dicke model in the USC regime.In the sensor atom approach (a), the addition of a second TLS shown as a sensor atom is weekly coupled to the cavity (hence shown outside of the cavity).In the generalized Dicke model (b), the second atom is also considered to be ultrastrongly coupled to the cavity (similar to the first atom, but it can have different coupling parameters).

Figure 2 .
Figure2.Cavity-QED eigenenergies with and without a sensor atom.Computed first ten eigenenergies of the QRM (blue solid lines) and the JCM (red dashed lines) for the cavity-QED system with (a) a single atom that is in resonance with the cavity, i.e., ωa = ωc, obtained from Eqs. (1) and (2), respectively, and (b) the sensor atom approach Hamiltonian, where both the primary and the sensor atoms are in resonance with the cavity, ωa = ωs = ωc, and the sensor atom has the coupling of g D s = 0.001g D ; and (c) the sensor atom approach Hamiltonian, where the primary atom is in resonance with the cavity, ωa = ωc, and the sensor atom has ωs = ωc/2 and g D s = 0.001g D , obtained from Eqs. (18) and(19), respectively.In panel (a), the three significant transitions are identified with the downward arrows and the letters 'A' (ω10/ωc ≈ 0.5), 'B' (ω31/ωc ≈ 8.2), and 'C' (ω20/ωc ≈ 1.45), for η = 0.5.These resonances are highlighted for reference when explaining the key features of the cavity spectrum.

Figure 3 .
Figure 3.Cavity-QED incoherent spectra for one atom in the ultrastrong coupling regime.Cavity emitted spectra were computed using the sensor atom approach (left column) and the full quantum regression theorem (right column).On the right, we list the type of bath used for the cavity and the two atoms (which use the same bath type).Gauge-corrected (not corrected) results are shown with solid blue (dashed orange) curves.We use incoherent driving with Pinc = 0.01g.Other system parameters are κ = 0.25g, γ = γs = 0.005g, gs = 0.001g, ωc = ωa, η = 0.5.The three significant peaks A, B, and C represent the major transitions shown in Fig.2.

Figure 4 .
Figure 4. Selected eigenenergies using the generalized Dicke model.We show the first seven eigenenergies of the GDM.(a) Full quantum model without a RWA (blue solid lines) and with a RWA (red dashed lines), displaying eigenenergies with ωa = 2ω b = ωc versus η ≡ ηa = η b are plotted.(b) Full model eigenenergies for ωa = 2ω b = ωc and |g b | = ga = 0.5ωc (η = 0.5 = |ηa| = |η b |, with g b = ga exp[iπϕ]) versus the variation of the relative phase between the coupling parameters are represented; here we do not show RWA results, as they are all clearly wrong and also gauge dependent.The random color coding in panel (b) is to help distinguish the different eigenenergies.
) the parities of the first seven eigenstates are shown for a range of interest in the second atom's normalized frequency, ω b /ω c .We define the parity of a state |j⟩ as ⟨j| P |j⟩ where P = exp[iπ N ] and N

Figure 5 .
Figure 5. Generalized Dicke model state parities, optical transitions and (normalized) transition rates.(a) Parities of the first seven states of the GDM versus the second atom's normalized frequency.(b) Eigenenergies of the first seven states of the GDM versus the second atom's normalized frequency with positive (blue) and negative (orange) parity.(c) P quadrature matrix element squared of the selected transitions in panel (b) versus the second atom's normalized frequency obtained via Eq.(16).The position of the arrows in panel (b) is irrelevant.The plots are for ηa = η b = 0.5.

Figure 6 .
Figure 6.Cavity spectra for the GDM, using two different atoms that are both ultrastrongly coupled to the cavity.(a) Cavity spectra computed using the quantum regression theorem with two USC atoms in GDM, the first on resonance (ωa = ωc) and the second (ω b ) we sweep through resonance.(b) Selected cavity spectra as in the lower left panel of (a) for ω b /ωc = {0.5, 1, 1.5}.The labeled peaks in ω b /ωc = 0.5 plot in (b) reflect on the corresponding transitions in Fig. 5c.We also plot the spectra obtained in the absence of the second TLS.All baths are Ohmic and here ηa = η b = 0.5.We use the following parameters throughout the section: κ = 0.25g, γa = γ b = 0.005g, and Pinc = 0.01g.

Figure 7 .
Figure 7. Influence of coupling strength on the GDM spectra.Cavity spectra as in Fig. 6 for various η values, with both atoms in the same coupling regime (i.e., η = ηa = η b ).As before, we keep the first atom on resonance with the cavity and sweep the second atom resonance.

Figure 8 .
Figure 8. Influence of relative coupling amplitude variation on the GDM cavity spectra.Spectra at selected ω b , where we now sweep g b from negligible coupling to the same level as the first atom, at η b = ηa = 0.5.

Figure 9 .
Figure 9. Relative coupling amplitude variation in GDM eigenenergies.Eigenvalues of the lowest seven states at selected ω b values, as a function of g b .The parity of the states is again given by the line color: blue (dark) represents even parity and orange (light) represents odd parity.

Figure 10 .
Figure10.Influence of relative coupling phase variation in GDM cavity spectra.Spectra at selected ω b , where we now sweep the phase of the second atom.Notably, when we lower the coupling strength below the USC regime ηa ≤ 0.1, there is no dependence on the phase.