Virtual photons in the ground state of a dissipative system

Much of the novel physics predicted to be observable in the ultrastrong light-matter coupling regime rests on the hybridisation between states with different numbers of excitations, leading to a population of virtual photons in the system's ground state. In this article, exploiting an exact diagonalization approach, we derive both analytical and numerical results for the population of virtual photons in presence of arbitrary losses. Specialising our results to the case of Lorentzian resonances we then show that the virtual photon population is only quantitatively affected by losses, even when those become the dominant energy scale. Our results demonstrate most of the ultrastrong-coupling phenomenology can be observed in loss-dominated systems which are not even in the standard strong coupling regime. We thus open the possibility to investigate ultrastrong-coupling physics to platforms that were previously considered unsuitable due to their large losses.

In their seminal work [1] Huttner and Barnett extended the microscopic quantum theory of light-matter interaction developed by Hopfield [2] to the case of homogeneous but lossy dielectrics. They accomplished this by quantising the electromagnetic field coupled to a dispersionless matter excitation, itself coupled to a continuum reservoir of harmonic oscillators acting as a bath in which energy could be dissipated. The full Hamiltonian of such a system can be written asĤ describes the free electromagnetic field, models the matter excitation and the bath, is the dipolar interaction between light and matter, and comes from the diamagnetic A 2 part of the the minimal-coupling Hamiltonian. In Supplementary Eqs. (2)-(5)â(k), b(k), andb ω (k), are bosonic annihilation operators respectively for a photon, a matter excitation, and an excitation of the bath with frequency ω, all indexed by the wavevector k, ω 0 is the frequency of the optically active transition, V (ω) models its coupling to the bath, ω c quantifies the intensity of the light-matter coupling, and the renormalised frequencyω 0 is linked to the bare one by the formulã Note that we can chose arbitrary large values for the light-matter coupling ω c because, thanks to the presence of the termĤ A 2 , the ground state is stable and no superradiant phase transition can take place, as demonstrated by various no-go theorems [4,5]. In Ref. [1] the authors make the choice to removeĤ A 2 by effectively performing a Bogoliubov rotation in the space of theâ(k), thus puttingĤ em +Ĥ A 2 in diagonal form in terms of the Bogoliubov-rotated operatorsâ(k) with ck = c 2 k 2 + ω 2 c . Unluckily this is not acceptable for us, because it implies that the very definition of the the bare photonic operators depends upon the strength of the coupling, and as such it becomes problematic to define the number of virtual photons in the ground state. Ignoring this problem leads to a number of inconsistencies, most notably the total energy per mode diverges for k → 0. In the following we will thus sketch a derivation that keepŝ H A 2 as part of the interaction. The Hamiltonian in Supplementary Eq. (1) can be diagonalised a la Fano [3] in two steps. FirstĤ mat is put into diagonal form where theB(k, ω) operators, describing the continuously broadened optically active matter resonance, obey bosonic commutation relations They can be expressed as linear combinations of the different uncoupled matter operators aŝ where the coefficients can be written as and Introducing the quantity which obeys the normalization the interaction part of the Hamiltonian, now describing the interaction of the photonic modes with the broadened transition, takes the form Analogously to the above, the novel Hamiltonian can now be put into diagonal form through the operatorŝ with the coefficientsα where the complex dielectric function is In comparison the equivalent of the coefficient in Supplementary Eq. (21) obtained in Ref. [1] has instead the form with the others following a similar scheme.

Supplementary Note 2: Virtual photon population in the lossless case
In order to verify our calculations, it can be of use to compare the results obtained through the dissipative theory in the case of vanishing losses with those obtained using the standard nondissipative theory due to Hopfield [2]. In order to do this we start from the equivalent of Supplementary Eq. (1) neglecting the batĥ

The Hamiltonian in Supplementary Eq. (27) can be put in the diagonal form
through the introduction of the polaritonic operatorŝ where and The transformation from bare to polaritonic basis in Supplementary Eq. (29) can be inverted aŝ and the photonic population in the mode k calculated with the Hopfield theory is thus We can now write N k in the limit of vanishing losses, and thus real eigenfrequencies, as Plugging Supplementary Eqs. (30)-(31) into Supplementary Eqs. (34)-(35), and developing the heavy but straightforward algebra we can then prove that the two expressions coincide, verifying the correctness of our procedure.

Supplementary Note 3: Lorentz dielectric function
Here we will prove that the Lorentz dielectric function where ω M is a cutoff frequency that we will eventually send to infinity and q a positive constant frequency. We will limit ourselves to consider the imaginary part of the dielectric function as the real part will be fixed by Kramers-Kronig relations. Plugging Supplementary Eq. (37) into Supplementary Eq. (6), and introducing the cutoff in the integration we obtaiñ From Supplementary Eq. (15), using the Sokhotski-Plemelj theorem where the P indicates the principal part of the integral. Sending the cutoff to infinity we thus have which is the imaginary part of the Lorentz dielectric function upon the identification

Supplementary Note 4: Virtual photon population with arbitrary losses
In the general situation in which photons can be lost both through absorption and leakage, we need to consider a more general HamiltonianĤ witĥ The Hamiltonian in Supplementary Eq. (47) is similar to the one in Supplementary Eq. (1) but the external bath, modelling in this case a continuum of photonic modes represented by bosonic operatorsâ ω (k), is now coupled through the coupling Q k (ω) to the photonic excitation whose renormalised frequency is given by The Hamiltonian in Supplementary Eq. (47) can be put in the diagonal form in terms of the operators describing the continuum spectrum of a leaky resonator where the coefficients areα and Introducing the quantity which obeys the normalization the interaction part of the Hamiltonian, describing two coupled continua, takes the form Contrary to the case treated previously, describing a single continuum coupled to a discrete resonance, the diagonalization of two coupled continua was not treated in the original paper by Fano [3] and it does not seem to be readily available in the literature. We will thus describe the procedure in more details. We aim to diagonalize the system in terms of two branches of continuum polaritonic modeŝ obeying bosonic commutation relations The eigenequation leads to the system where the other two coefficients obey Defining the functions where P denotes the principal part, we can formally solve for the coefficients as with s x,j,k (ω) and s w,j,k (ω) functions to be determined. Multiplying Supplementary Eq. (71) by ω ′ ω+ω ′ , integrating, and exploiting Supplementary Eq. (67), we finally get to the equation In order to fix the extra function introduced in Supplementary Eqs. (70)-(71) we can now exploit Supplementary Eq. (61), which after some algebra leads to ω 2 c |χ k (ω)| 2 π 2 + s * w,j ′ ,k (ω)s w,j,k (ω) We have at this point determined five equations, two from Supplementary Eq. (72) and three from Supplementary Eq. (73), in the six unknown functions K j,k (ω), s w,j,k (ω), s x,j,k (ω), with j = ±. The two eigenmodesP ± (k, ω) are degenerate for each value of k and ω. This implies a basis in such a degenerate subspace needs to be chosen through gauge fixing, leading to a sixth equation which then allows to algebraically solve the system and complete the diagonalization. We will fix this gauge freedom by choosing the basis in which w +,k (ω, ω) = 0, which from Supplementary Eq. (71) implies s w,+,k = 0. The number of ground state photons can at this point be readily determined by the formula In order to get quantitative results from Supplementary Eq. (74) we now need to fix the coupling function Q k (ω). The ultrastrong coupling regime has been achieved using many different kinds of photonic resonators, each one described by a different Q k (ω). In order to achieve an acceptable level of generality, as done for the matter losses, we specialise the theory to the case of a Lorentzian resonance, which is usually an acceptable approximation for most real resonators in a quite broad parameter range. We thus consider a coupling of the form of the one in Supplementary Eq. (37) Calculating all the relevant integrals, and letting ω M → ∞ at the end, we recover the normalised density for the photonic mode with broadening γ P ω|χ k (ω)| 2 = 2γ P π ω 2 (c 2 k 2 − ω 2 ) 2 + γ 2 P ω 2 . (77) In the Supplementary Fig. 1 we plot the virtual photonic population of the resonant mode as a function of γ L and γ P for different values of ω c . We can clearly see that in the considered parameter range N k essentially depends on the total broadening γ L + γ P , and that a sizeable virtual photon population, equal to roughly the 50% of the lossless one, remains even when γ L = γ P = γ max and light and matter resonances are both overdamped.