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 us 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 in both the Power–Zienau–Woolley and the Coulomb representations. In order to exemplify the results, the theory is finally specialized to two specific cases. In the first one, both light and matter resonances are characterized by Lorentzian linewidths, and in the second one, a fixed absorption band is also present. The analytical expressions derived 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. A Matlab code implementing our results is provided.

## 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 system approaches have been developed,^{1} non-perturbative diagonalization is possible using a procedure developed by 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 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 a harmonic oscillator. The light is coupled with a discrete optical resonance of the material, itself modeled as a harmonic degree of freedom and coupled to a 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 the polariton, the larger losses it will incur, with pure photons very detuned from the material resonance propagating unimpeded in the dielectric.

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 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–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–14}

This large interest is motivated by the increasing experimental relevance of a rigorous treatment of lossy CQED systems, including the impact of frequency-dependent 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–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–26} and experimental^{27–30} efforts having studied the possibility of strong coupling taking the continuum into account.

An analytical solution extending 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 line shape 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 such a contribution and develop a full, general, and useable 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 and 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 been shown to be not always correct at arbitrarily large coupling strengths^{33–35} for systems in which only a few material resonances are considered. Finally, the theory previously 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.

In this paper, we will implement these improvements by developing explicitly the theory for the diagonalization of the light–matter system in the Power–Zienau–Woolley (PZW) representation, describing the interaction between the polarization field **P** and the electric displacement **D**. We will demonstrate that an exact diagonalization of the problem for arbitrary couplings and environments is possible, and we will provide analytical expressions for the physical observables. Most importantly, these expressions will be gauge-invariant, not depending upon the extra degrees of freedom. Our results will thus provide 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 detail. The equivalent results for the Coulomb representation are reported for completeness in Appendix C. All the mathematical symbols used are listed in Table I.

Symbol . | Description . | Depending on . | First appearance in (equation number) . |
---|---|---|---|

k | Composite wavevector index | (1) | |

ω_{k} | Bare photon frequency | k | (1) |

a_{k} | Photon annihilation operator | k | (1) |

$a\u0303k$ | Rotated photon operator | k | (C3) |

ω_{x} | Bare matter frequency | (1) | |

b_{k} | Matter annihilation operator | k | (1) |

$b\u0303k$ | Rotated matter annihilation operator | k | (3) |

g | Light–matter coupling strength | (1) | |

$\omega \u0303x$ | Matter frequency renormalized by the diamagnetic term | (1) | |

ω_{±,k} | Polariton frequencies | k | (8) |

α_{k} | Photonic reservoir annihilation operators | k and ω | (9) |

β_{k} | Matter reservoir annihilation operators | k and ω | (9) |

V_{k} | Photonic reservoir interaction function | k and ω | (11) |

Q | Matter reservoir interaction function | ω | (12) |

$\omega \u0304k$ | Photon frequency renormalized by the reservoir | k | (14) |

$\omega \u0304x$ | Matter frequency renormalized by the reservoir | k | (14) |

A_{k} | Broadened photon operator | k and ω | (15) |

x_{k}, z_{k} | Photonic mode Hopfield coefficients | k and ω | (15) |

y_{k}, w_{k} | Photonic reservoir Hopfield coefficients | k, ω and ω′ | (15) |

γ_{k} | Unknown function | k and ω | (21) |

χ_{k} | Self-energy term for the photonic reservoir | k and ω | (24) |

B_{k} | Broadened matter operator | k and ω | (25) |

x, z | Matter mode Hopfield coefficients | ω | (25) |

y, w | Matter reservoir Hopfield coefficients | ω, ω′ | (25) |

t | Self-energy term for the matter reservoir | ω | (27) |

ζ_{k} | Expansion coefficient of the photonic field | k and ω | (29) |

η | Expansion coefficient of the matter field | ω | (29) |

P_{k,j} | Polariton operator | k, branch j = ± and ω | (32) |

$x\u0304k,j,y\u0304k,j,w\u0304k,j,z\u0304k,j$ | Polariton Hopfield coefficients | k, branch j = ± and ω, ω′ | (32) |

K_{k,j} | Integral function of the photonic coefficient | k, branch j = ± and ω | (35) |

J_{k,j} | Integral function of the matter coefficient | k, branch j = ± and ω | (35) |

Z | Integral function of |η(ω)| | k and ω | (36) |

W_{k} | Integral function of |ζ_{k}(ω)| | k and ω | (36) |

$Xj,k,Yj,k,Zj,k$, $Wj,k$ | Hopfield coefficients of broadened polaritons | k, branch j = ± and ω | (41) |

γ_{P} | Photonic reservoir loss rate | (50) | |

ω_{P} | Photonic reservoir cut-off frequency | (50) | |

γ_{M} | Matter reservoir loss rate | (50) | |

ω_{M} | Matter reservoir cut-off frequency | (50) | |

$Vk1$ | Interaction function between the light and absorption band | k and ω | (54) |

κ | Interaction strength to the reservoir | (55) | |

$\omega \u0304k$ | Photonic frequency renormalized by the reservoir and continuum | k | (55) |

$\omega \u0303k$ | Photonic frequency renormalized by the continuum | k | (57) |

F | Normalized coupling density to the reservoir | ω | (55) |

Ω_{k} | Effective central frequency for the added reservoir | k and ω | (61) |

Γ | Effective loss rate | ω | (61) |

Symbol . | Description . | Depending on . | First appearance in (equation number) . |
---|---|---|---|

k | Composite wavevector index | (1) | |

ω_{k} | Bare photon frequency | k | (1) |

a_{k} | Photon annihilation operator | k | (1) |

$a\u0303k$ | Rotated photon operator | k | (C3) |

ω_{x} | Bare matter frequency | (1) | |

b_{k} | Matter annihilation operator | k | (1) |

$b\u0303k$ | Rotated matter annihilation operator | k | (3) |

g | Light–matter coupling strength | (1) | |

$\omega \u0303x$ | Matter frequency renormalized by the diamagnetic term | (1) | |

ω_{±,k} | Polariton frequencies | k | (8) |

α_{k} | Photonic reservoir annihilation operators | k and ω | (9) |

β_{k} | Matter reservoir annihilation operators | k and ω | (9) |

V_{k} | Photonic reservoir interaction function | k and ω | (11) |

Q | Matter reservoir interaction function | ω | (12) |

$\omega \u0304k$ | Photon frequency renormalized by the reservoir | k | (14) |

$\omega \u0304x$ | Matter frequency renormalized by the reservoir | k | (14) |

A_{k} | Broadened photon operator | k and ω | (15) |

x_{k}, z_{k} | Photonic mode Hopfield coefficients | k and ω | (15) |

y_{k}, w_{k} | Photonic reservoir Hopfield coefficients | k, ω and ω′ | (15) |

γ_{k} | Unknown function | k and ω | (21) |

χ_{k} | Self-energy term for the photonic reservoir | k and ω | (24) |

B_{k} | Broadened matter operator | k and ω | (25) |

x, z | Matter mode Hopfield coefficients | ω | (25) |

y, w | Matter reservoir Hopfield coefficients | ω, ω′ | (25) |

t | Self-energy term for the matter reservoir | ω | (27) |

ζ_{k} | Expansion coefficient of the photonic field | k and ω | (29) |

η | Expansion coefficient of the matter field | ω | (29) |

P_{k,j} | Polariton operator | k, branch j = ± and ω | (32) |

$x\u0304k,j,y\u0304k,j,w\u0304k,j,z\u0304k,j$ | Polariton Hopfield coefficients | k, branch j = ± and ω, ω′ | (32) |

K_{k,j} | Integral function of the photonic coefficient | k, branch j = ± and ω | (35) |

J_{k,j} | Integral function of the matter coefficient | k, branch j = ± and ω | (35) |

Z | Integral function of |η(ω)| | k and ω | (36) |

W_{k} | Integral function of |ζ_{k}(ω)| | k and ω | (36) |

$Xj,k,Yj,k,Zj,k$, $Wj,k$ | Hopfield coefficients of broadened polaritons | k, branch j = ± and ω | (41) |

γ_{P} | Photonic reservoir loss rate | (50) | |

ω_{P} | Photonic reservoir cut-off frequency | (50) | |

γ_{M} | Matter reservoir loss rate | (50) | |

ω_{M} | Matter reservoir cut-off frequency | (50) | |

$Vk1$ | Interaction function between the light and absorption band | k and ω | (54) |

κ | Interaction strength to the reservoir | (55) | |

$\omega \u0304k$ | Photonic frequency renormalized by the reservoir and continuum | k | (55) |

$\omega \u0303k$ | Photonic frequency renormalized by the continuum | k | (57) |

F | Normalized coupling density to the reservoir | ω | (55) |

Ω_{k} | Effective central frequency for the added reservoir | k and ω | (61) |

Γ | Effective loss rate | ω | (61) |

The paper is organized as follows: In Sec. II, we introduce and diagonalize the dissipationless polaritonic Hamiltonian. This will be useful in introducing the problem and extracting 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 independently diagonalize 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 specialize 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 also incorporates 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* being the resonant light–matter coupling strength and the operators obeying bosonic commutator relations,

where *δ*_{k,k′} is the Kronecker delta, and here and in the following, we will consider *ℏ* = 1. 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,

with renormalized matter frequency

where the rotated operators still obey bosonic commutation relations,

The Hamiltonian in Eq. (1) then takes the simpler form

The secular equation of such a 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 ensembles 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

The frequency shifts in Eq. (14), usually referred to as Lamb shifts due to their similarity with the well-known atomic physics effect, can be derived from the Lagrangian of a discrete resonance linearly coupled with a continuum, as done in Ref. 3.

We can diagonalize the photonic Hamiltonian *H*_{PB} in Eq. (11) introducing the bosonic operators describing broadened photons *A*_{k}(*ω*), written as linear superposition of the uncoupled operators,

whose coefficients can be found via HB diagonalization, illustrated in more detail in Appendix B. We will generally refer to the coefficients of the linear decompositions of bosonic operators as Hopfield coefficients.^{38} From the eigenequation

the resulting system reads

Such a system cannot be trivially solved eliminating one-by-one its unknowns because under the hypothesis that the eigenfrequency *ω* falls into the photonic reservoir continuum, there will always be a value of *ω*′ = *ω*, which makes the left-hand-side of Eq. (19) vanish. This is in stark contrast with the discrete case in which coupled modes are never degenerate with bare resonances.^{39,40} 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

The need to introduce the unknown function *γ*_{k}(*ω*) stems from the fact that *y*_{k}(*ω*, *ω*), function of *ω* and *k*, is not determined by Eq. (19) as its left-hand-side identically vanishes.

After some manipulations, we can solve the system in Eqs. (17)–(20) arriving at the following expressions for the coefficients

with

and we defined $Vk(\omega )$ 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. (12) describing the matter sector *H*_{MB}, by introducing the bosonic operator for the broadened optically active resonance,

The solution is in the analogous form

with

As for its photon counterpart, we defined $Q(\omega )$ the odd analytic extension of |*Q*(*ω*)|^{2} in the negative frequency range.

We can now recover the inverse transformations for the bare operators in terms 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. (29) into the coupling Hamiltonian from Eq. (13), the full light–matter 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 diagonalized by introducing the operators for two polariton branches *P*_{j}(*ω*) with *j* = ±,

defined as arbitrary linear superpositions of the bare modes at all frequencies. By exploiting the eigenequation

In order to put the system in Eq. (34) 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. (29).

We then solve Eq. (34) for the unknown coefficients,

The crucial difference between this system of equations and the one obtained in the single-continuum case in Eqs. (17)–(20) is that here both bare modes have continuum spectra and thus both the first and the second equations in Eq. (34) diverge. Therefore, we have to introduce two unknown functions *s*_{x,k,j}(*ω*) and *s*_{y,k,j}(*ω*), analogous to *γ*_{k}(*ω*) introduced in Eq. (21), for each value of frequency, wavevector, and each polaritonic branch. From Eq. (38), 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 *k*- and *ω*-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\u22121)2$ orthogonality conditions, leaving $L(L\u22121)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. (39) unchanged.

We can also express the polariton operators as superpositions of the bare ones inserting Eqs. (15) and (25) into Eq. (32),

with the most relevant Hopfield coefficients having the form

representing the single resonant and anti-resonant matter and photon components of the polariton modes.

Writing explicitly the coefficients as in Eq. (41), we can at this point find the relations linking the light and matter weight of the polaritons to the Hopfield coefficients,

The bare photonic and matter field operators can then be expressed in terms of the broadened polaritons,

recognizing in the expressions in Eq. (35) the photonic and matter components of the broadened polaritonic modes. By expressing the coefficient $y\u0303k,j(\omega ,\omega \u2032)$ as in Eq. (38), we also arrive at a relation between the three quantities defined in Eq. (35),

while from Eq. (39), we derive the three equations

for *j*, *j*′ = −, *j*, *j*′ = +, and *j* = −, *j*′ = +, respectively. We are thus left, as anticipated, with an underdetermined set of five equations [two from Eq. (46) and three from Eq. (47)] in six unknowns: the *K*_{k,j}(*ω*) and the two 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. (46) and (47) 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,

## 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 HB theory.^{41} This is, nevertheless, false as we will now demonstrate constructively. The problem in Ref. 41 is that the authors do not recognize the need to introduce a divergent renormalized frequency. We will consider couplings of the form

with $qP=\pi \omega k2/2\gamma P$ and $qM=\pi \omega \u0303x2/2\gamma M$, where we have introduced cut-off frequencies *ω*_{P} and *ω*_{M}, and the photonic and matter loss rates *γ*_{P} and *γ*_{M}, and Θ is the Heaviside function. In Eq. (50) also appear the renormalized frequencies as from Eq. (14)

The form of the couplings in Eq. (50) 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 detail in Appendix D, but by inserting Eq. (50) into Eq. (30), we obtain

from which we recognize the two Lorentzian lineforms that would be obtained from a classical Lorentz dielectric model with center frequency *ω*_{k} or $\omega \u0303x$ and width *γ*_{P} or *γ*_{M}. Hence, according to Eq. (36), we can calculate the other two required functions,

Note that, as made clear from Eq. (51), the frequencies $\omega \u0304k$ and $\omega \u0304x$ appearing in the Hamiltonians in Eqs. (11)–(13) diverge with the cutoffs. Following the basic idea behind renormalization, we can recognize this not to be an issue because they are not observables of our system, while the proper physical observables, *Z*(*ω*) and *W*_{k}(*ω*) in Eq. (53), feature only the finite frequencies *ω*_{k} and $\omega \u0303x$.

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 Figs. 1–5, the black dotted line represents the resonant frequency of the material resonance $\omega \u0303x$, the red dashed line represents the photonic frequency *ω*_{k}, and the blue dashed-dotted lines represent the polaritonic resonances in the lossless case from Eq. (8). In Figs. 2 and 3, we plot the gauge-independent results from Eq. (49) for different values of, respectively, the losses and the coupling strength, showing that the theory behaves as predicted from input–output theories,^{5–7} with two polaritonic Lorentzian resonances with linewidths proportional to the weighted average of the light and matter respective linewidths.

## VI. DIAGONALIZATION WITH COLORED RESERVOIRS

The theory we developed allows us to model polaritonic systems with arbitrary colored reservoirs, including the scientifically and technologically relevant case of a continuum absorption band, a case object of multiple theoretical^{22–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 lifetime and the necessity of keeping both can cause some formal problem given that, as we saw before, the modeling 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 colored reservoir on the top of the Lorentzian one. The microscopic reservoir modes leading to the Lorentzian line shape of the uncoupled photonic resonance are *a priori* completely uncorrelated from those leading to an absorption band. They will in many cases 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 $Vk0(\omega )$ and $Vk1(\omega )$ model, respectively, the interaction with the photonic reservoir and the absorption band, defined by a normalized density *F*(*ω*) with *F*(*ω* < 0) = 0, $\u222b0\u221ed\omega F(\omega )=1$, and a dimensionless coupling intensity *κ*. The two coupling functions take the form

where $|Vk0(\omega )|2$, which will lead to the Lorentzian broadening as in the previous case, has the same form as in Eq. (50), 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 $\omega \u0303k$ dressed only by the absorption band,

and define the latter as

We then follow the same diagonalization procedure as in Sec. III, with function *χ*_{k}(*ω*) now taking the form

where $Vk0(\omega )$ and $Vk1(\omega )$ are the odd analytical extensions in the negative frequency range of $|Vk0(\omega )|2$ and $|Vk1(\omega )|2$, respectively. The squared function |*ζ*_{k}(*ω*)|^{2} can be written as in Eq. (30),

leading to, after some algebra, and after letting *ω*_{P} → ∞,

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 $\gamma P=\pi \omega \u0303k22qP$. Note that in the low-frequency regime *ω* → 0 from Eq. (61), we have Ω_{k} ≈ *ω*_{k}, showing that the presence of the absorption band does not change the background permittivity.

From Eq. (60), 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 light of recent studies 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 renormalized expressions in Eq. (61), the integral function *W*_{k}(*ω*) in Eq. (36) can be simply evaluated numerically.

In Fig. 4, we plot the example results obtained using a sharp absorption band of center frequency *ω*_{c} and width Δ,

The boundaries of the band are marked by the horizontal green dashed lines. We recognize the effects expected from our analytical results: The band acts as a localized 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 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 perform 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 line shape for strongly coupled resonances interacting with reservoirs of arbitrary spectral shapes. We then showed how colored features can be practically added to a 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 finite size effects can *a priori* be included in the theory as higher order terms.^{42} The same is true for generic terms present in the PZW Hamiltonian nonlinear in the light and matter fields, which can be written using Eq. (43) as nonlinear polaritonic interactions.

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, characterized 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.^{26,28} In order to facilitate the use of our results by the community, we uploaded on GitHub a Matlab code implementing them.^{43}

## ACKNOWLEDGMENTS

S.D.L. is a Royal Society Research Fellow and was partly funded by the Philip Leverhulme Prize of the Leverhulme Trust. The authors acknowledge funding from the RGF\EA\181001 grant from the Royal Society and the Leverhulme Trust under Grant No. RPG-2019-174.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no Conflict of Interest to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request. The program used to generate all the data is available online in Ref. 43.

### APPENDIX A: THE FANO DIAGONALIZATION AND ITS EXTENSION

In this appendix, we briefly discuss the basic idea behind the Fano diagonalization and its extensions to the case of many continua or many discrete levels. Following Ref. 2, we consider the Hamiltonian

where in the original paper *a* and *b*(*ω*) in Eq. (A1) are normalized 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 diagonalized in terms of a hybridized 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.^{44} In the bound-to-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. Proceeding by iteration, the system can be solved in terms of a single hybridized continuum of the form

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 transformation

where *U*_{mn}(*ω*) is a unitary matrix whose first row is given by

with

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. (11), describing the interaction between the discrete photon mode and the photonic reservoir, modeled as an ensemble of harmonic oscillators indexed by the continuum frequency *ω*.

We introduce the bosonic operators describing broadened photons *A*_{k}(*ω*),

whose coefficients are chosen so that the operators satisfy the eigenequation

This equation leads to the system between the coefficients

This set of equations 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

The function *γ*_{k}(*ω*) can be found, after some algebra, replacing both the equations in Eq. (B8) with Eq. (B3),

where we assume that the analytic extension in the negative frequency range $Vk(\omega )$ 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 Eqs. (B4) and (B8), in terms of *x*_{k}(*ω*), Eq. (B10) leads to the definition of *x*_{k}(*ω*) up to a phase factor,

with

The final expressions for the coefficients are thus obtained as

We have thus demonstrated that the operators *A*_{k}(*ω*) are eigensolutions of *H*_{PB}. The proof these operators form a complete set of solutions in the case of an unbounded continuum can be found in the original HB paper^{3} in which the linear transformation in Eq. (B1) is inverted, expressing the discrete *a*_{k} operator as a linear superposition of *A*_{k}(*ω*).

### APPENDIX C: DIAGONALIZATION IN THE COULOMB REPRESENTATION

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 rotated photon operators

and renormalized cavity frequency

The dispersion relation obtained diagonalizing the Hamiltonian in Eq. (C2) is

with solutions

In the Coulomb representation, the expressions of the field operators in Eq. (29) are modified as

where

By substituting the bare operators with the dressed ones, we arrive at the Hamiltonian as in Eq. (31), which can be diagonalized by the same procedure described in the main text. The dressed photonic and matter field operator can be finally written as superpositions of polaritonic broadened modes as

where

Note that in the Coulomb representation, the functions related to the photonic component and those related to the matter part have exchanged units from those in the PZW representation. This is due to the inverted dependence of the light and matter fields upon their frequency.

### APPENDIX D: DERIVATION OF ANALYTICAL RESULTS FOR A LORENTZIAN BROADENING

In this appendix, we will derive the analytical expression of the functions *ζ*_{k}(*ω*) and *η*(*ω*) in the case of a Lorentzian broadening. This is not completely trivial, as testified by the existence of a published paper claiming it is impossible.^{41} The key issue is that in order to recover a frequency-independent broadening, the shift in the mode due to the coupling with the reservoir has to diverge. As such, the result can only be found by a renormalization procedure allowing to cancel such a divergence.

Assuming that the coupling potentials to the photonic and matter reservoirs take the form in Eq. (50), we can calculate the dressed resonance frequencies

and the real and imaginary parts of functions *χ*_{k}(*ω*) and *t*(*ω*) as

In the limit of infinite cut-off frequency *ω*_{P} → ∞ and *ω*_{M} → ∞, the resonant frequencies diverge, but the intensity of the coupling vanishes, and we arrive at the finite results,

By inserting Eq. (D3) into Eq. (30), we finally obtain the Lorentzian form for the functions *ζ*_{k}(*ω*) and *η*(*ω*) as