Exact solution of polaritonic systems with arbitrary light and matter frequency-dependent losses

In this paper we perform the exact diagonalization of a light-matter strongly coupled system taking into account arbitrary losses via both energy dissipation in the optically active material and photon escape out of the resonator. This allows to naturally treat the cases of couplings with structured reservoirs, which can strongly impact the polaritonic response via frequency-dependent losses or discrete-to-continuum strong coupling. We discuss the emergent gauge freedom of the resulting theory and provide analytical expressions for all the gauge-invariant observables both in the Power-Zienau-Woolley and the Coulomb representations. In order to exemplify the results the theory is finally specialised to two specific cases. In the first one both light and matter resonances are characterised by Lorentzian linewidths, and in the second one a fixed absorption band is also present. The analytical expressions provided in this paper can be used to predict, fit, and interpret results from polaritonic experiments with arbitrary values of the light-matter coupling and with losses of arbitrary intensity and spectral shape, in both the light and matter channels.


I. INTRODUCTION
The interaction between discrete energy levels and degrees of freedom with continuum spectra is crucial to the description of any real-world quantum system, in which the coupling with the environment eventually leads to energy and information leakage. While many powerful perturbative open quantum systems approaches have been developed 1 , non-perturbative diagonalization is possible using a procedure due to Fano 2 . In his landmark paper Fano considered the problem of one discrete level coupled to one continuum. In the same paper he then moves to consider the cases of multiple discrete levels coupled to one continuum and of one discrete level coupled to multiple continua, showing that both cases can be reduced to the initial one. A short summary of Fano's approach and its generalizations is given in Appendix A.
One important application of the Fano's theory in the many discrete levels-one continuum case is light interacting with a dissipative dielectric, originally developed by Huttner and Barnett (HB) 3 . In such a formalism light with a well-defined momentum and polarization propagating in a bulk dielectric is represented in second quantization as an harmonic oscillator. The light is coupled with a discrete optical resonance of the material, itself modeled as an harmonic degree of freedom and coupled to an harmonic reservoir leading to dissipation. By diagonalizing the light-matter Hamiltonian one finds two hybrid polaritonic branches, which in the following we will call lower (-) and upper (+) polaritons, coupled to a reservoir through their matter component. As expected, the more matter-like is the polariton the larger losses it will incur, with pure photons very detuned from the material resonance propagating unimpeded in the dielectric. a) Electronic mail: E-mail: e.cortese@soton.ac.uk b) Electronic mail: E-mail: s. de-liberato@soton.ac.uk Complications arise in systems with boundaries, as traditional cavity quantum electrodynamics (CQED) setups or surface modes. An HB-like diagonalization can still be performed in real space if the photons are supposed to be perfectly trapped in a finite volume 4 , but in the general case a novel dissipation pathway opens, this time linked with the photonic component of the polaritons: photons can escape out of the system coupling with the free-space photonic continuum.
Such a setup is thus described by two discrete resonances (the photons and the optically active resonance) coupled to two different continua (the material reservoir and the continuum of free space photons). This case was not explicitly treated in the Fano's original paper and, as we will see, it is not possible to trivially apply the method adopted in the other cases. Still, various approximate approaches have been developed to deal with open CQED. Input-output approaches integrate out the system in order to describe relations between the incoming and outgoing fields [5][6][7] . Master equations integrate out the environment 8,9 , or at least most of it 10 , to describe the internal system dynamics. Some approaches exactly solve the coupling with the propagative electromagnetic field (radiative broadening), while describing phenomenologically matter losses 11,12 . It is also possible to use quasinormal mode quantization in order to quantize directly the lossy electromagnetic field [12][13][14] .
This large interest is motivated by the increasing experimental relevance of a rigorous treatment of lossy CQED systems, including the impact of frequencydependent structured environments. Ever larger values of the light-matter interaction energy 15 have in-fact allowed us to access non-perturbative coupling regimes as the ultrastrong [16][17][18][19] or the very strong ones 20,21 . In these regimes the polaritonic modal shifts are comparable to other energy scales, and frequency-independent approximation can dramatically fail. In particular polaritonic discrete resonances can interact with continua, with multiple theoretical 22 polariton operator k, branch j = ± and ω Eq. (30) x k,j ,ȳ k,j ,w k,j ,z k,j polariton Hopfield coefficients k, branch j = ± and ω, ω Eq. (30) K k,j integral function of the photonic coefficient k, branch j = ± and ω Eq. (33) J k,j integral function of the matter coefficient k, branch j = ± and ω Eq. ing studied the possibility of strong coupling taking the continuum into account.
An analytical solutions extending the Fano's approach would be useful in this context, in order to be able to study the quantum properties of systems in the presence of generic couplings, environments, and loss channels. It would allow for quantitative modeling of the lineshape of plasmonic systems once the loss channels in the metals are known 31 . Such an approach was derived in Ref. 32 to calculate the quantum properties of the ground state at arbitrary values of the system-reservoir coupling. In this work it was shown how an unphysical degree of freedom appears in the theory due to the presence of two coupling continua, and how a solution can be obtained through an arbitrary gauge fixing. The method has also been more recently used in Ref. 30 to reproduce experimental data in which polaritonic nonlocality created a broad absorption band above the bare photonic frequency.
Our aim in this paper is to improve over such a contribution and develop a full, general, and usable analytic theory for polaritonic systems with arbitrary couplings to the environment. Such an improvement will take four different forms. The first and non-negligible one is that the theory will be clearly laid down in the paper, while in Refs. 30,32 the derivation is at most sketched. The second, more important reason is that the pure gauge nature of the extra degree of freedom was not proven but only assumed. The third is that the theory was developed in the Coulomb representation, which has since then be shown to be not always correct at arbitrarily large coupling strengths [33][34][35] for systems in which only few material resonances are considered. Finally, the theory pre-viously used, although applicable to model reservoirs of arbitrary spectral shapes, requires a renormalization procedure. It is thus not directly applicable to cases beyond the presence of a simple Lorentzian broadening.
We will implement these improvements by developing explicitly the theory in the Power-Zienau-Woolley (PZW) representation. Calculating results for the gaugeinvariant physical observables without fixing the extra degrees of freedom we will show their fully gauge nature. Our results will thus provide a proof of the approach from Ref. 32 , while demonstrating gauge fixing is not necessary. Finally we will provide a recipe to add arbitrary frequency-dependent reservoirs. Given this work's technical nature we describe the calculations in details. Equivalent results for the Coulomb representation are reported for completeness in the Appendix C. All the mathematical symbols used are listed in Table I.
The paper is organised as follows. In Sec. II we introduce and diagonalise the dissipationless polaritonic Hamiltonian. This will be useful to introduce the problem and to extract the discrete polaritonic dispersion which we will then use to interpret the results of the dissipative theory. Note that in this paper we will always start from Hamiltonians. Their derivation from Lagrangians can be found in Ref. 4 for the PZW case and in Ref. 3 for the Coulomb case. In Sec. III we introduce the light and matter reservoirs and we diagonalise independently the light and matter sectors of the full Hamiltonian, into a broadened matter resonance interacting with a broadened photonic one. In Sec. IV we solve the full Hamiltonian, describing the problem with gauge ambiguities, and derive expressions for the gauge-invariant observables. In Sec. V we specialise the model to the case of Lorentzian resonances. In Sec. VI we provide the recipe to add arbitrary reservoirs to the Lorentzian broadening, and present the results in the case of a fixed absorption band.

II. DIAGONALIZATION IN THE LOSSLESS CASE
We start here by introducing and diagonalizing the lossless Hamiltonian for a photon field of dispersion ω k indexed by the composite wavevector index k, which incorporates also polarization and all other relevant conserved quantum numbers. Such a field is described by the bosonic annihilation operator a k coupled to a matter excitation, of fixed frequency ω x , described by the bosonic annihilation operator b k . We neglect here nonlocality due to the dispersion of the material resonance 30,36,37 .
The PZW light-matter Hamiltonian for the system described above is with g the resonant light-matter coupling strength and the operators obeying bosonic commutator relations with δ k,k a Kronecker delta. The second term in Eq. (1) is the square P 2 term which can be removed by performing a Bogoliubov rotation of the matter component of the Hamiltonian, renormalizing the bare matter frequency, as The Hamiltonian in Eq. (1) then takes the simpler form The secular equation of such an Hamiltonian reads leading to the frequencies of the two polariton branches (±) for each value of the wavevector k

III. DIAGONALIZATION OF PHOTONIC AND MATTER RESERVOIRS
We now pass to introduce the dissipation in the picture by defining two reservoirs in which photons and matter excitations can be lost. Those reservoirs, which model respectively the continuum of free-space photons and the continuum of phononic and electronic degrees of freedom in the material are modeled as ensambles of harmonic oscillators indexed by the continuum frequency ω. Their annihilation operators are α k (ω) and β k (ω) respectively, with bosonic commutator relations The total Hamiltonian now takes the form with the Hamiltonians for the photonic and matter sectors and their interaction in the form where V k (ω) and Q(ω) are the interaction functions modeling the interaction of, respectively, the photonic mode and the matter excitation with their respective reservoirs, and the bare light and matter resonances are dressed by the coupling as We can diagonalise the photonic Hamiltonian H PB in Eq. (9) introducing the bosonic operators describing broadened photons A k (ω), whose coefficients can be found via HB diagonalization, illustrated in more details in Appendix B. From the eigenequation the resulting system reads Such a system cannot be trivially solved eliminating oneby-one its unknowns because, under the hypotesis that the eigenfrequency ω falls into the photonic reservoir continuum, there will always be a value of ω = ω which makes the left-hand-side of and Eq. (17) vanish. This is in stark contrast with the discrete case in which coupled modes are never degenerate with bare resonances 38,39 .
The system can nevertheless be solved in the distribution sense as where P is the principal value and γ k (ω) is an unknown function, which can be fixed by imposing the bosonic commutation relation After some manipulations we can solve the system in Eqs. (15)- (18) arriving to the following expressions for the coefficients and we defined V k (ω) the odd analytic extension of |V k (ω)| 2 in the negative frequency range. Exactly the same procedure can be applied to the Hamiltonian in Eq. (10) describing the matter sector H MB , by introducing the bosonic operator for the broadened optically active resonance The solution is in the analogous form As for its photon counterpart, we defined Q(ω) the odd analytic extension of |Q(ω)| 2 in the negative frequency range.
We can now recover the inverse transformations for the bare operators in term of the broadened ones, and write the bare field operators as superpositions of the broadened fields where the expansion coefficients of the bare fields upon the broadened operators are given by the expressions

IV. DIAGONALIZATION OF THE FULL LIGHT-MATTER HAMILTONIAN
After substituting the field operators in Eq. (27) into the coupling Hamiltonian from Eq. (11), the full lightmatter Hamiltonian can be written in the form which describes the broadened photonic mode coupled to the broadened material resonance. Similarly to what done previously, the Hamiltonian can be diagonalised by introducing the operators for two polariton branches defined as arbitrary linear superpositions of the bare modes at all frequencies. By exploiting the eigenequation we arrive to a system of equations analogous to Eqs. (15)-(18) In order to put the system in Eq. (32) in a form apt to be manipulated and solved we introduce two unknown integral functions of the diagonalization coefficients, which as we will see play the role of photonic and material amplitudes of the polaritonic field and two known integral functions of the coupling coefficients Note that, notwithstanding the apparent symmetry, the functions related to the photonic component |K k,j (ω)| 2 and W k (ω), have different units from those of the matter part |J k,j (ω)| 2 and Z(ω). The former are pure numbers, while the latter are times squared. This is due to the specific dependence of the light and matter fields upon their frequency in the PZW representation, clearly visible in Eq. (27).
We then solve Eq. (32) for the unknown coefficients The crucial difference between this system of equations and the one obtained in the single-continuum case in Eqs. (15)- (18) is that here both bare modes have continuum spectra and thus both the first and the second equations in Eq. (32) diverge. We have therefore to introduce two unknown functions s x,k,j (ω) and s y,k,j (ω) for each value of frequency, wavevector, and for each polaritonic branch. From Eq. (36) we can see that those functions allow us to arbitrarily fix the equal-frequency mixing between coupled and uncoupled modes. We are thus led to add four different functions at fixed wavevector and frequency, but the commutator relations represent only three new equations, for j, j = − and j, j = + (normalization) and j = −, j = + (orthogonality). This leaves a free function corresponding to a kand ω-dependent rotation in the space of the two degenerate polaritonic modes. More generally in the presence of L continua, we would add L 2 unknown functions, and obtain L normalization conditions and L(L−1) 2 orthogonality conditions, leaving L(L−1) 2 quantities to be determined, which is the dimension of the O(L) group. An element of such a group corresponds to a rigid rotation in the space of the L P k,j (ω) modes at fixed k and ω, which would leave Eq. (37) unchanged.
We can also express the polariton operators in terms of the bare modes inserting Eq. (13) and Eq. (23) into Eq. (30) with the most relevant coefficients having the form Writing explicitly the coefficients as in Eqs. (39) we can at this point find the relations linking the light and matter component of the polaritons to the relevant Hopfield coefficients The bare photonic and matter field operators can then be expressed in terms of the broadened polaritons recognising in the expressions in Eq. (33) the photonic and matter components of the broadened polaritonic modes. By expressing the coefficientỹ k,j (ω, ω ) as in Eq. (36), we also arrive to a relation between the three quantities defined in Eq. (33) Inserting the second of Eq. (36) into the first of Eq. (33), we derive the two relations (j = ±) while from Eq. (37) we derive the three equations g 2 |η(ω)| 2 π 2 + s y,k,j (ω)s * y,k,j (ω) J k,j (ω)J * k,j (ω)+ (45) |ζ k (ω)| 2 π 2 + s x,k,j (ω)s * x,k,j (ω) K k,j (ω)K * k,j (ω) = δ j,j , for j, j = −, j, j = + and j = −, j = + respectively. We are thus left, as anticipated, with an underdetermined set of 5 equations (2 from Eq. (44) and 3 from Eq. (45)) in six unknowns: the K k,j (ω) and the 2 functions s x,k,j (ω), s y,k,j (ω) for each of the two values of j. In Ref. 32 this problem was solved arbitrarily fixing s x,k,− (ω) = 0 and then solving the remaining equations. Here instead we adopt a different approach, solving directly Eqs. (44)-(45) for the gauge-invariant quantities. Due to the arbitrariness of the gauge choice, it is in-fact meaningless to distinguish between lower and upper polariton operators as only the gauge invariant quantities are the total photonic and material polaritonic components Although the solution is algebraically cumbersome, it is possible to find the analytic expressions for the total light and matter polaritonic components, which are the central result of this paper .
Note that the gauge extra variable disappears as expected from Eq. (47), thus proving its pure gauge nature and as a consequence the correctness of the results from Ref. 32 , where those same quantities had been calculated by an arbitrary gauge fixing.

V. LORENTZIAN RESONANCES
In order to analytically evaluate the functions introduced above, we need to specify a model for the coupling between the bare modes and the photonic and matter reservoirs. We start by writing a model for homogeneously broadened resonances, thus reproducing the optical response of a Lorentz dielectric model. There have been claims that such a dielectric model cannot be modeled in the HB theory 40 . This is nevertheless false as we will now demonstrate constructively. The problem in Ref. 40 is that the authors do not recognise the need to introduce a divergent renormalized frequency. We will consider couplings of the form with q P = πω 2 k /2γ P and q M = πω 2 x /2γ M , where we have introduced cut off frequencies ω P and ω M , the photonic and matter loss rates, γ P and γ M and Θ is the Heaviside function. In Eq. (48) also appear the renormalised frequencies as from Eq. (12) The form of the couplings in Eq. (48) has been chosen to recover, in the limit of diverging cutoff frequencies, Lorentzian resonances centered at the bare excitation frequencies. This is shown in more details in Appendix D, but by inserting Eq. (48) in Eq. (28) we obtain lim from which we recognise the two Lorentzian lineforms which would be obtained from a classical Lorentz dielectric model with center frequency ω k orω x and width γ P or γ M . Notice that the resonances are centered at the bare frequencies, not the ones renormalised from the interaction with the reservoirs in Eq. (49), which instead diverge with the cutoff frequency. Hence, according to Eq. (34) we can calculate the other two required functions Using these explicit forms for the couplings in Fig. 1 we plot the photonic and material component of each polaritonic branch obtained with the gauge-fixing s x,k,− (ω) = 0 used in Ref. 32 . As we can see the distinction between the operators of the two polaritonic branches, represented in the first two columns, is arbitrary and we cannot identify them with specific resonances. Only their sum, the total intensity of the photonic or matter components, shown in the third column, has physical meaning. In all the figures the black dotted line represents the resonant frequency of the material resonanceω x , the red dashed line the photonic frequency ω k , and the black dash-dotted lines the polaritonic resonances in the lossless case from Eq. (6). In Fig. 2 and Fig. 3 we plot the gaugeindependent results from the Eq. (47) for different values of, respectively, the losses and the coupling strength, showing that the theory behaves as predicted from inputoutput theories 5-7 , with two polaritonic Lorentzian resonances with linewdiths proportional to the weighted average of the light and matter respective linewidths.

VI. DIAGONALIZATION WITH COLOURED RESERVOIRS
The theory we developed allows to model polaritonic systems with arbitrary coloured reservoirs, including the scientifically and technologically relevant case of a continuum absorption band, a case object of multiple theoretical [22][23][24][25][26] and experimental 27-30 studies. Focusing for the sake of definiteness on a compactly-supported reservoir interacting with the photonic component of the excitation (an absorption band), we can include it in the theory by choosing frequency-dependent coupling functions V k (ω) with support in the chosen frequency range. Exactly the same procedure would couple to the matter component by using Q(ω) instead.
The presence of an absorption band will generally not substitute other loss channels influencing the excitations and ω +,k are marked by a dot-dashed blue lines. Although we maintain the notation K ±,k and J ±,k , it is clear that it is no longer possible to isolate the contributions of the two polariton branches lifetime and the necessity of keeping both can cause some formal problem given that, as we saw before, the mod-  elling of a Lorentzian resonance requires us to renormalize an otherwise infinite resonant frequency. In this Section we will provide the recipe to add an arbitrary coloured reservoir on the top of the Lorentzian one. The microscopic reservoir modes leading to the Lorentzian lineshape of the uncoupled photonic resonance are a priori completely uncorrelated from those leading to an absorption band. They will in many case have even different physical origins, e.g., phonon interaction for the former and interaction with electronic bands for the latter. As better explained in Appendix A this means their effects sum incoherently and we can thus write the full interaction of the photonic component with its environment using the coupling function where V 0 k (ω) and V 1 k (ω) model, respectively, the interaction with the photonic reservoir and the absorption band, defined by a normalised density F (ω) with F (ω < 0) = 0, ∞ 0 dωF (ω) = 1, and a dimensionless coupling intensity κ. The two coupling functions take the form where |V 0 k (ω)| 2 , which will lead to the Lorentzian broadening as in the previous case, has the same form as in Eq. (48), and the other term is chosen in order to provide the required absorption band after the renormalization.
The dressed frequency of the photonic mode is now renormalized by both terms. It is practical to write the impact of the Lorentzian broadening as a renormalization over the frequency of the photonic resonanceω k dressed only by the absorption band and define the latter as We then follow the same diagonalization procedure as in Section III, with function χ k (ω) now taking the form where V 0 k (ω) and V 1 k (ω) are the odd analytical extensions in the negative frequency range of, respectively, |V 0 k (ω)| 2 and |V 1 k (ω)| 2 . The squared function |ζ k (ω)| 2 can be written as in Eq. (28) leading to, after some algebra, and after letting ω P → ∞, to with effective central frequency and effective losses where the linewidth of the Lorentzian losses is now defined as a function of the frequency dressed by the absorption band γ P = πω 2 k 2q P . Note that in the low-frequency regime ω → 0 from Eq. (59) we have Ω k ≈ ω k , showing that the presence of the absorption band does not change the background permittivity.
From Eq. (58) we see that the photonic losses increase in the frequency region in which F (ω) = 0, an effect already observed in Ref. 30 , and it also leads to a resonance effect in the central frequency. This can be understood in the light of recent works on strong coupling with the continuum 24,26 , and we expect it to model the possibility of the photonic resonance becoming strongly coupled with the absorption band. From the renormalised expressions in Eq. (59), the integral function W k (ω) in Eq. (34) can be simply evaluated numerically.
In Fig. VI we plot example results obtained using a sharp absorption band of center frequency ω c and width ∆ The boundaries of the band are marked by horizontal green dashed lines. We recognise the effects expected from our analytical results: the band acts as a localised absorber for the photonic component of the polariton and eventually the polaritonic mode gets strongly coupled to the band, an effect better visible when the band width becomes comparable with the intrinsic photonic linewidth.

VII. CONCLUSIONS
In this article we exactly solved the polaritonic problem with a quantum formalism in the case of arbitrary dissipative couplings for both the bare photonic and matter excitations. In order to do this we discussed the extension of the Fano and HB theories to the case of multiple discrete levels coupled to multiple continua, showing how a gauge indeterminacy emerges. While a previous approach to this problem had been to peform an arbitrary gauge choice, here we analytically calculated the gauge-invariant observables. We thus demonstrated both the self-consistency of our theory and provided an analytical, albeit cumbersome formula allowing to exactly calculate the resonance lineshape for strongly coupled resonances interacting with reservoirs of arbitrary spectral shapes. We then showed how coloured features can pratically be added to an homogeneous resonance linewidth. Note that while our approach is based on a purely bosonic spectrum of the material resonance, generally correct for plasmonic and phononic systems, saturation and finite size effects can a priori be included in the theory as higher order terms 41 .
We hope these results will be of use in the subfields of polaritonic in which losses have an important impact. This is true for example in plasmonic systems, characterised by important Ohmic losses, nonlocal systems where photonic excitations couple to a continuum of propagative modes, and in systems in which the extremely large coupling between light and matter pushes the polaritonic resonances into other spectral features. In this Appendix we briefly discuss the basic idea behind the Fano diagonalization and its extenstions to the case of many cantinua or many discrete levels. Following Ref. 2

we consider the Hamiltonian
where in the original paper a and b(ω) in Eq. (A1) are normalised Hilbert space vectors, but in the present context they can as well be interpreted as second-quantized annihilation operators. Under the assumption that all the coupled eigenvalues fall inside the initial continuum range Fano showed how the system can be exactly diagonalised in term of an hybridised continuum The discrete mode thus gets dressed by a cloud of continuum excitations, translating into a spectral broadening of the resonance. Notice that this set of solution is not necessarily complete, as known from the study of the Friederichs-Lee model 42 . In the boundto-continuum strong coupling regime discrete modes can emerge from the continuum, as theoretically and experimentally demonstrated in the case of two-dimensional electron gases 24,29 . After having completed the diagonalization procedure, Fano passes to consider the case in which there are N discrete levels and one continuum. Such a problem can be reduced to the one treated above by initially performing a partial diagonalization of one discrete level coupled to the continuum, leading to a novel Hamiltonian in the same form as the initial one but this time with N − 1 discrete levels. Proceding by iteration the system can be solved in term of a single hybridised continuum of the form x n (ω)a n + dω y(ω, ω )b(ω ). (A3) Finally, the case of a single discrete state coupled to N continua is treated, described by the Hamiltonian Such a system can be solved by performing the transformationb where U mn (ω) is a unitary matrix whose first row is given by transforming Eq. (A4) into the Hamiltonian of one discrete level coupled to a single continuum, plus other N −1 uncoupled continua

Appendix B: Huttner-Barnett Diagonalization
In this Appendix we perform the HB diagonalization of the Hamiltonian H PB in Eq. (9), describing the interaction between the discrete photon mode and the photonic reservoir, modeled as an ensamble of harmonic oscillators indexed by the continuum frequency ω.
We introduce the bosonic operators describing broadened photons A k (ω), A k (ω) = x k (ω)a k + z k (ω)a † k + (B1) dω y k (ω, ω )α k (ω ) + w k (ω, ω )α † k (ω ) , whose coefficients are chosen so that the operators satisfy the eigenequation This equation leads to the system between the coefficients This set of equation can be solved to obtain z k (ω), y k (ω, ω ) and z k (ω, ω ) in terms of x k (ω). By subtracting Eq. (B4) from Eq. (B3), we obtain which can be substituted in Eq. (B5) and Eq. (B6) to obtain The function γ k (ω) can be found, after some algebra, replacing both the equations in Eq. (B8) in Eq. (B3) where we assume that the analytic extension in the negative frequency range V k (ω) of |V k (ω)| 2 is an odd function. In order to calculate x k (ω) we impose the commutation relation By using the expression for A k (ω) in Eq. (B1) and for the coefficients in Eq. (B4) and Eq. (B8), in terms of x k (ω), Eq. (B10) leads to definition of x k (ω) up to a phase factor Exploiting the expression for γ k (ω) in Eq. (B9), Eq. (B11) can be written as with The final expressions for the coefficients are thus obtained as .
Appendix C: Diagonalization in the Coulomb representation FIG. 5. Coupled photonic K 2 k (a) and matter J 2 k (b) fields as a function of the bare cavity frequency ω k , calculated starting from a Coulomb representation Hamiltonian. The light-matter couling is g = 0.3ωx, while the loss rates are γP = γM = 0.05ωx. All the other parameters remain as in Fig. 1.
In this Appendix, we show how our approach gets modified if one wants to work in the Coulomb representation. Given that only relatively few quantities are affected by the change, we just provide the expressions for the affected ones. The Hamiltonian in the Coulomb representation can be written as in which we can see the appearance of a diamagnetic A 2 term. After reabsorbing such a diamagnetic term by performing a Bogoliubov transformation, the Hamiltonian takes the form with renormalised cavity frequencỹ