# Relaxation breakdown and resonant tunneling in ultrastrong-coupling cavity
QED

Daniele De Bernardis1 1Pitaevskii BEC Center, CNR-INO and Dipartimento di
Fisica, Università di Trento, I-38123 Trento, Italy.

###### Abstract

We study the open relaxation dynamics of an asymmetric dipole that is
ultrastrongly coupled to a single electromagnetic cavity mode. By using a
thermalizing master equation for the whole interacting system we derive a
phase diagram of the Liouvillian gap. It emerges that the ultrastrong coupling
inhibits the system’s relaxation toward the equilibrium state due to an
exponential suppression of the dipole tunneling rate. However, we find that
polaronic multi-photon resonances restore fast relaxation by a cavity-mediated
dipole resonant tunneling process. Aside of the numerical evidences, we
develop a fully analytical description by diagonalizing the Rabi model through
a generalized rotating-wave approximation, valid in the so-called polaron
frame. The relaxation physics of such ultrastrong-coupling systems is then
reduced to a multi-photon polaron version of the standard text-book dressed
states picture. At the end we discuss an extension to a multi-well dipole that
can set the basis of a cascaded resonant tunnelling setup in the ultrastrong
coupling regime.

##  I Introduction

Relaxation from a metastable state toward equilibrium is a central problem in
many branches of physics, such as chemical reactions, radiaoactive decay and
electronic transport, to name a few Landau and Lifshitz (2013). The energy
barrier separating a local minimum from the stable equilibrium, i.e. the
activation barrier of chemical reactions Piskulich  _et al._ (2019), can be
overcome by thermal fluctuations, for which, after an initial absorption of
energy from the bath, the system is kicked out the metastable state, rolling
down to its absolute equilibrium state and releasing the energy excess.

In quantum physics similar processes can happen driven only by quantum
fluctuations. When the temperature is too small to kick the system over the
metastable energy barrier, relaxation is then dominated by the _tunnel effect_
(or quantum tunneling), which is one of the first surprising consequences of
the quantum theory Merzbacher (2002).

Following the hand-wavy intuition that quantum fluctuations replace thermal
ones in kicking the system out of the metastability, one might speculate that
including in these systems a supplemental quantum reservoir could sensibly
alter the tunneling dynamics. The work of Leggett et al. on tunneling-systems
coupled to an environment Caldeira and Leggett (1981); Leggett  _et al._
(1987) has shown that this is actually the case, and tunneling can be sensibly
changed as a function of the environment parameters. Since for most systems
the natural environment is provided by the electromagnetic radiation, here
quantum tunneling is crossing its path with another fundamental concept of
quantum physics: the non empty vacuum of quantum electrodynamics (QED) Milonni
(1994), rising the question: can vacuum fluctuations of the electromagnetic
field affect tunneling and relaxation in material systems?

Considering only free-space radiation the immediate answer is that the
electromagnetic vacuum does not give any sensible contribution in the
tunneling dynamics, because vacuum effects scale with the small fine-structure
constant $\alpha_{\rm fs}\simeq 1/137$, giving rise to only minor
modifications of materials Sáez-Blázquez  _et al._ (2022).

However, when the electromagnetic field is confined, such as in cavity QED
setups, vacuum effects are boosted by an effective fine-structure constant
$\alpha_{\rm eff}>1$ Devoret  _et al._ (2007); De Bernardis  _et al._ (2018a),
entering in the field of the so-called _ultrastrong coupling regime_ Frisk
Kockum  _et al._ (2019); Forn-Díaz  _et al._ (2019a). Here the electromagnetic
vacuum is not empty anymore, being heavily populated by virtual photons Ciuti
_et al._ (2005). This suggests an important interplay between the coupling to
a cavity field and tunneling relaxation.

Experiments have already shown that this intuition may be substantial, in
particular in chemical reaction and electronic transport where strong
differences are observed when molecules, atoms or electrons are enclosed in
such extreme resonant electromagnetic environment Hutchison  _et al._ (2012);
Paravicini-Bagliani  _et al._ (2019a); Valmorra  _et al._ (2021); Appugliese
_et al._ (2022).

All these exciting observations have stimulated multiple theoretical debates
in various communities opening new research lines such as: polaritonic
chemestry Flick  _et al._ (2017); Fregoni  _et al._ (2022); Schäfer  _et al._
(2022), cavity QED control of electronic transport in mesoscopic devices or in
quantum Hall systems Paravicini-Bagliani  _et al._ (2019b); Arwas and Ciuti
(2023), cavity QED modification of ferromagnetism, ferroelectricity and
superconductivity De Bernardis  _et al._ (2018a); Schuler  _et al._ (2020);
Ashida  _et al._ (2020); Andolina  _et al._ (2020); Román-Roche  _et al._
(2021); Schlawin  _et al._ (2019), all with the general aim to explore and
understand up to which degree the quantum vacuum of cavity QED can be a
resource to modify and control properties of matter Schlawin  _et al._ (2022);
Bloch  _et al._ (2022); Garcia-Vidal  _et al._ (2021).

Because of the complexity in treating the non-perturbative nature of these
light-matter ultrastrong-coupling systems, at the beginning of this debate the
focus was mainly on closed systems, trying to understand the correct way to
write down these Hamiltonians with the only intention of exploring their
ground-state physics. However, recently there has been a growing interest also
in their out-of-equilibrium dynamics Hausinger and Grifoni (2008); Mercurio
_et al._ (2022); Settineri  _et al._ (2018); Schaeverbeke  _et al._ (2019);
Magazzù  _et al._ (2021); Chen  _et al._ (2022), which is still a largely
unexplored territory. Here even the simplest case regarding relaxation toward
thermal equilibrium can exhibit surprising non-trivial features.

In this article we address exactly this problem studying how the
electromagnetic vacuum of ultrastrong-coupling cavity QED can affect the
relaxation toward equilibrium of a polarizable material. In order to isolate
every single different effect we consider a simple paradigmatic setup: an
asymmetric double well dipole in a single-mode resonant cavity. Its low-energy
dynamics can be approximated to the _quantum Rabi model_ , which is the
simplest theoretical framework to study light-matter interactions. We complete
the description of the model including two basic dissipative mechanism: Ohmic
cavity dissipation and dipole radiative losses. Under these circumstances we
describe the system relaxation dynamics through a master equation valid for
arbitrary light-matter coupling values, whose steady state is the correct
thermal equilibrium state. From the spectral gap $\lambda$ of its Liouvillian
operator we derive a phase diagram describing how relaxation toward
equilibrium is modified by the coupling to the cavity.

Differently from previous studies Hausinger and Grifoni (2008); Settineri  _et
al._ (2018); Mercurio  _et al._ (2022), here we exploit a generalized
rotating-wave approximation of the Rabi model from which we analytically
derive the transition rates of the master equation in the ultrastrong coupling
regime. From this calculation we obtain a complete and simple picture on how
relaxation and thermalization work in terms of polaronic dressed states, valid
in the ultrastrong coupling regime.

As a result, we show that the ultrastrong coupling exponentially suppresses
the relaxation rate of the dipole, being a prototype for the so-called
localization transition in the spin-boson model Leggett  _et al._ (1987).
However, special resonances exist where relaxation is restored thanks to a
cavity-mediated multi-photon resonant tunneling process, in very close analogy
to Frank-Condon physics describing electron tunneling assisted by vibrational
transitions Koch  _et al._ (2006); Leturcq  _et al._ (2009); Cui  _et al._
(2015); Vdovin  _et al._ (2016). After showing that this mechanism is already
observable in current experimental platforms such as superconducting circuits
we comment on the possible consequence for cavity assisted quantum transport
and cascaded ultrastrong-coupling setups with multi-well dipoles.

The article is organized as follows. In Sec. II we introduce the physical
system and its approximated description in terms of the asymmetric quantum
Rabi model. By considering the Liouvillan gap of its open dynamics, in Sec.
III we study how the ultrastrong coupling regime changes the relaxation and
thermalization rate. By using a generalized rotating-wave approximation to
diagonalize the Rabi model we explicitly show an exponential slow-down of the
system’s relaxation due to the ultrastrong coupling regime. In Sec. IV we show
that the fast relaxation can be restored by a cavity assisted multi-photon
resonant tunnelling process. Exploiting again the generalized rotating-wave
approximation we develop the discussion in terms of multi-photon polaron
dressed states. In Sec. V we extend this setup to the extended Dicke model
leading to a cascaded resonant tunnelling device. Finally, in Sec. VI we draw
our conclusions.

##  II Model

Figure 1: (a) Cavity QED system. The cavity is modelled as an LC-circuit,
where the inductor magnetic flux $\Phi$ takes the role of the dynamical
variable of the electromagnetic field, usually given by the vector potential
$\vec{A}$. The dipole inside the capacitor couples to the voltage drop between
the plates, $U=\dot{\Phi}$. (b) The dipole is described as a particle in a
tilted double-well potential. The position $x$ represent the dipole
displacement and thus its dipole moment. When the central well of the
potential is large enough the system is approximated by only the two lowest
levels (two-level approximation). (c) Open-system schematic view. The cavity
QED system can be interpreted as an element of a dissipative circuit.

We consider the paradigmatic cavity quantum electrodynamics (cQED) setup
described in Fig. 1(a), where a single electrically polarizable object (a
dipole) is placed into the planar capacitor of a resonant LC circuit. This
system is among the simplest toy models able to reproduce most of the features
of cQED in all various coupling regimes, and is particularly important in
giving a simple and intuitive description of many solid-state or circuit cQED
setups relevant for experiment in the ultrastrong coupling (USC) regime in the
GHz or THz range Jaako  _et al._ (2016); De Bernardis  _et al._ (2018a); Forn-
Díaz  _et al._ (2019b); Frisk Kockum  _et al._ (2019).

###  II.1 A dipole in a asymmetric double well potential

The dipole dynamics is described as a single particle with mass $m$ in a
potential

$H_{\rm dipole}=\frac{p^{2}}{2m}+V(x).$ (1)

The dipole has charge $+q$ on an extreme and $-q$ on the other extreme, and
$qx$ is its dipole moment. So, the dipole displacement $x$ is its main
dynamical variable and $p$ its canonical momentum.

As depicted in Fig. 1(b), we consider only the paradigmatic case of a dipole
described by a double-well potential, very similarly to Refs. De Bernardis
_et al._ (2018b); Di Stefano  _et al._ (2019). Considering only its low energy
dynamics we are basically studying the electromagnetism of quantum tunnelling.

Differently from Ref. De Bernardis  _et al._ (2018a) here we introduce a
linear tilt that makes the height of the two wells asymmetric. The total
dipole’s potential is then given by

$V(x)=-\frac{\mu_{2}^{2}}{2}x^{2}+\frac{\mu_{4}^{4}}{4}x^{4}+q\mathcal{E}x.$
(2)

The linear tilt $\sim q\mathcal{E}x$ is physically implemented by a bias
static external electric field of amplitude $\mathcal{E}$ and has no influence
on the LC resonator dynamics.

Whenever $\mu_{2}/\mu_{4}\gg 1$ and $q\mathcal{E}/\mu_{2}\ll 1$, the two
lowest levels are below the central barrier and well separated from the other
energy levels. We can then truncate the dipole’s Hilbert space keeping only
these two lowest energy levels De Bernardis  _et al._ (2018a). We perform the
two level approximation (TLA) on the dipole Hamiltonian projecting on the
eigenstates without tilt $\mathcal{E}=0$, and we obtain

$H_{\rm dipole}^{\rm TLA}=\hbar\omega_{d}s_{z}+\hbar\epsilon s_{x}.$ (3)

Here we have introduced the pseudospin operators $s_{a}=\sigma_{a}/2$
($\sigma_{a}$ are the usual Pauli matrices). The dipole frequency $\omega_{d}$
is the energy difference between the two lowest states without the tilt
$\mathcal{E}=0$, and $\epsilon=2q\mathcal{E}x_{10}/\hbar$, where
$x_{10}=\braket{1}{x}{0}$ is the dipole matrix element between the two lowest
dipole states. The dipole operator now is given by $x\approx
x_{10}\sigma_{x}$.

When the tilt is on, the energy splitting between the two eigenstates of Eq.
(3) is given by $\omega_{\epsilon}=\sqrt{\omega_{d}^{2}+\epsilon^{2}}$, while
the eigenfunctions are partially localized on the left or right well, with a
small, but non-negligible, overlap with the opposite well, Fig. 1(b). In the
two-level language these states are given by

$\begin{split}&|{\rm
L}_{\epsilon}\rangle=\cos\frac{\theta_{\epsilon}}{2}|\downarrow\rangle+\sin\frac{\theta_{\epsilon}}{2}|\uparrow\rangle\\\
&|{\rm
R}_{\epsilon}\rangle=-\sin\frac{\theta_{\epsilon}}{2}|\downarrow\rangle+\cos\frac{\theta_{\epsilon}}{2}|\uparrow\rangle,\end{split}$
(4)

where $\tan\left(\theta_{\epsilon}\right)=\epsilon/\omega_{d}$.

###  II.2 General cavity QED Hamiltonian

The Hamiltonian of the full cavity QED system can be written summing up the
dipole energy and the total energy stored in the electromagnetic field

$H_{\rm cQED}=H_{\rm em}+H_{\rm dipole}.$ (5)

For an LC-resonant system, the electromagnetic energy is described by

$H_{\rm em}=\frac{CU^{2}}{2}+\frac{\Phi^{2}}{2L},$ (6)

where $U$ is the total voltage drop across the capacitor $C$, $\Phi$ is the
magnetic flux through the inductance $L$.

When the dipole is inside the capacitor the total voltage $U$ is no more the
right canonical variable conjugate to $\Phi$. In order to have the correct
canonical description we need to introduce the total capacitor charge variable
$Q$, such that $\left[\Phi,Q\right]=i\hbar$. Without the presence of the
dipole the total charge and the voltage drop are directly proportional through
the usual relation $U=CQ$. When the dipole is inside the capacitor the charge
responsible for the voltage drop is modified by the presence of the charge
induced by the dipole on the metallic plates. This induced charge does not
contribute to any voltage drop, and must be removed Jaako  _et al._ (2016); De
Bernardis  _et al._ (2018a). The voltage drop becomes

$U=C\left(Q-Q_{\rm ind}\right).$ (7)

For an ideal capacitor we have $Q_{\rm ind}\simeq qx/d$ De Bernardis  _et al._
(2018a), where $d$ is the distance between the capacitor plates.

We introduce now the cavity creation/annihilation operators through the
relations

$\begin{split}&\Phi=i\frac{\Phi_{0}}{\sqrt{2}}\left(a-a^{{\dagger}}\right)\\\
&Q=\frac{Q_{0}}{\sqrt{2}}\left(a+a^{{\dagger}}\right).\end{split}$ (8)

Here $\Phi_{0}=\sqrt{\hbar Z_{LC}}$, $Q_{0}=\sqrt{\hbar/Z_{LC}}$ and
$Z_{LC}=\sqrt{L/C}$. The cavity QED reads

$\begin{split}H_{\rm cQED}=\hbar\omega_{c}a^{{\dagger}}a+H_{\rm
dipole}+F_{0}x\left(a+a^{{\dagger}}\right)+\frac{F_{0}^{2}}{\hbar\omega_{c}}x^{2},\end{split}$
(9)

where $\omega_{c}=1/\sqrt{LC}$ and we introduced the zero-point electric force
$F_{0}=\sqrt{\hbar\omega_{c}/(2Cd^{2})q^{2}}$.

As detailed in Refs. De Bernardis  _et al._ (2018b); Di Stefano  _et al._
(2019), implementing the TLA described in Section II.1, we can now approximate
the cavity QED Hamiltonian with the so-called _quantum Rabi model_ ($\hbar=1$)

$\begin{split}H_{\rm cQED}\approx H_{\rm
Rabi}=\omega_{c}a^{{\dagger}}a+\omega_{d}s_{z}+\epsilon
s_{x}+g\left(a+a^{{\dagger}}\right)s_{x},\end{split}$ (10)

where the light-matter coupling is given by $g=2F_{0}x_{10}$.

We will see in particular that the transverse term $\sim\epsilon s_{x}$, that
breaks the $\mathbb{Z}_{2}$ symmetry of the usual Rabi model, is fundamental
in our development, becoming a switch between slow and fast dissipation of the
dipole. In circuit QED this term emerges quite naturally through a bias in the
external magnetic flux and is typically used in the observation of the
spectral features of the USC regime Yoshihara  _et al._ (2017).

###  II.3 Cavity dissipation

A standard way to introduce dissipation in an LC circuit is to couple it to a
transmission line Gardiner and Zoller (2015). When is traced out from the
dynamics, the transmission line plays the role of a resistive element,
effectively realizing the scheme described in Fig. 1(c). An input voltage can
inject current in the system, and the resistance damps the excited
oscillations of the LC-circuit. This scheme provides a basic input/output
theory describing the system’s read out from the energy dissipated into the
resistance.

The formal description of this setup assumes a linear coupling to a multi-mode
bath of harmonic modes, as described in Eq. (50) of Appendix A. To realize the
resistive circuit, since $\Phi=LI$ is linked to the current passing through
the inductor (and so through the whole circuit) the system operator coupled to
the bath is Gardiner and Zoller (2015)

$\hat{X}=\Phi.$ (11)

To reproduce a standard Ohmic resistance we assume the standard linear
spectral density

$J_{\rm Ohm}(\omega)\simeq\gamma\omega.$ (12)

Eliminating the bath’s dynamics, we obtain the equations of motion of the
circuit in terms of the Langevin equation in Eq. (55), and, considering
$A=\dot{\Phi}$ in Eq. (55), we have that

$\ddot{\Phi}\sim-\gamma\dot{\Phi},$ (13)

correctly matching the standard Kirchhoff equations of a resistive circuit.

With the supplemental assumption of small coupling between the circuit and the
resistance, the same type of calculation leads to the Lindblad master equation
for the density matrix of the system, which will be used in the remaining of
the paper.

###  II.4 Dipole dissipation

In our simplified picture, the main source of dissipation of an oscillating
electric dipole is given by radiative emission into free space electromagnetic
modes. Indeed, even if the dipole is strongly coupled to the sub-wavelength
mode of the LC cavity, still interacts with the all the other transverse
electromagnetic modes. These have typically a small effect on the coherent
dynamics Sáez-Blázquez  _et al._ (2022), but they provide a decay channel for
the dipole.

Because of the harmonic dynamics of the electromagnetic field, and its general
linear coupling with the dipole, we can again model the dipole dissipation as
a linear damping, as described in Appendix A. In this case the system operator
coupled to the dissipative bath is given by the dipole moment

$\hat{X}=x\approx 2x_{10}s_{x}$ (14)

The bath spectral function is given from the spectral function of the
transverse electromagnetic modes. In free-space this would be given
approximately by

$J_{\rm rad}(\omega)\simeq\kappa\omega^{3}.$ (15)

Considering the Langevin equation in Eq. (55) for the dipole moment velocity
$A=\dot{x}$, using this spectral density we find

$\ddot{x}\sim\kappa\dddot{x},$ (16)

recovering the Abraham-Lorentz formula for a radiating dipole Dalibard  _et
al._ (1982).

More generally for our developments, one can choose any spectral density for
the dipole dissipative bath of the shape $J_{\rm rad}\sim\omega^{\nu}$, with
$\nu\geq 1$. We will see that having $\nu>1$ gives particularly simple
results.

##  III Relaxation regimes of cavity QED

In this section we will explore the combined effect of light-matter coupling
$g$ and dipole asymmetry $\epsilon$ on the open relaxation dynamics of the
system. To this aim we will use a thermalizing master equation, as detailed in
Appendix B, implementing the dissipation mechanisms described before, by
including $J_{\rm Ohm}(\omega)$ and $J_{\rm rad}(\omega)$ to model the
cavity/dipole bath spectral densities.

Moreover, in order to have a master equation that correctly describes
relaxation toward the equilibrium state, the cavity/dipole jump operators are
rewritten in the total system eigenbasis, reading

$\begin{split}&c_{nm}^{c}=\braket{n}{a-a^{{\dagger}}}{m}|n\rangle\langle
m|,\\\ &c_{nm}^{\rm dip}=\braket{n}{s_{x}}{m}|n\rangle\langle m|.\end{split}$
(17)

Here $|n\rangle$ and $|m\rangle$ are the eigenstates of the Rabi Hamiltonian
in Eq. (10).

Since in our theory the coefficient $\gamma,\kappa$ are free parameters, we
absorb the dimensional quantities $\Phi_{0}/\sqrt{2}\,,\,2x_{10}$ in the
definition of the spectral densities, such that

$\frac{\omega_{c}|\Phi_{0}|^{2}}{2}J_{\rm
Ohm}(\omega)\,,~{}4\omega_{d}|x_{10}|^{2}J_{\rm rad}(\omega)\longmapsto J_{\rm
Ohm}(\omega),J_{\rm rad}(\omega).$ (18)

The full system master equation can be written as the sum of three Liouvillian
super operators

$\partial_{t}\rho=\mathcal{L}_{H}(\rho)+\mathcal{L}_{c}(\rho)+\mathcal{L}_{\rm
dip}(\rho).$ (19)

Here $\mathcal{L}_{H}(\rho)=-i\left[H_{\rm sys},\rho\right]$ generates the
coherent time evolution, while $\mathcal{L}_{c}$ and $\mathcal{L}_{\rm dip}$
are given by Eq. (63) with the respective transition rates

$\begin{split}&\Gamma^{c}_{nm}=\gamma\frac{|\omega_{mn}|}{\omega_{c}}|\braket{n}{a-a^{{\dagger}}}{m}|^{2},\\\
&\Gamma_{nm}^{\rm
dip}=\kappa\frac{|\omega_{mn}|^{3}}{\omega_{d}^{3}}|\braket{n}{s_{x}}{m}|^{2},\end{split}$
(20)

where $\omega_{mn}=\omega_{m}-\omega_{n}$ is the difference between the
eigenfrequencies of the Rabi Hamiltonian in Eq. (10).

For the sake of simplicity, through the whole manuscript we only focus on the
relevant dipole-cavity resonant case where $\omega_{c}=\omega_{d}$.

###  III.1 Zero temperature Liouvillian gap

Figure 2: (a) Phase diagram of the Liouvillian gap $\lambda$ as a function of
the light-matter coupling $g$ and the dipole asymmetry $\epsilon$. Parameters:
$\gamma=\kappa/4=0.05\omega_{c}$, $\omega_{d}=\omega_{c}$. (b) A cut of the
phase diagram at $\epsilon=0$ as a function of the light-matter coupling $g$,
in logscale. (c) A cut of the phase diagram at $g/\omega_{c}=3$ as a function
of the dipole asymmetry $\epsilon$, in logscale.

To study the relaxation toward equilibrium of the system we consider the
Liouvillian gap $\lambda={\rm Re}[\lambda_{1}]$ Kessler  _et al._ (2012);
Minganti  _et al._ (2018), obtained from the spectrum $\\{{\lambda_{n}\\}}$,
$n=0,1,2\ldots$, of the total Liouvillian operator
$\mathcal{L}=\mathcal{L}_{H}+\mathcal{L}_{c}+\mathcal{L}_{\rm dip}$ defined
from Eq. (19) Johansson  _et al._ (2013). Here we consider the zero
temperature case $T=0$ only, which is the relevant case for superconducting
cavity QED setups Yoshihara  _et al._ (2017). The same picture holds also for
finite temperature, provided that
$k_{b}T\lesssim\hbar\omega_{c},\hbar\omega_{d}$, where $k_{b}$ is the
Boltzmann constant. When the temperature grows larger, and
$k_{b}T>\hbar\omega_{c}$ USC effects are pushed to much larger light-matter
coupling values Pilar  _et al._ (2020).

In Fig. 2(a) we show the Liouvillan gap as a function of the light-matter
coupling and the dipole asymmetry, $\lambda(g,\epsilon)$.

At small light-matter coupling $g\sim 0$, the effect of increasing $\epsilon$
is to make the dipole an eigenstate of $s_{x}$, decreasing the value of the
matrix element in the dipole transition rate in Eq. (20). However the
vanishing matrix element is compensated by the increasing energy difference
between the dipole levels, giving larger controbution from the radiative
spectral density of the bath, $J_{\rm rad}\sim\omega^{3}$. This can be seen by
explicitly computing the dipole transition rate using the bare uncoupled
dipole states in Eq. (4), for which we have

$\begin{split}\Gamma_{\rm LR}^{\rm
dip}&=\frac{\kappa}{4}\left(1+\frac{\epsilon^{2}}{\omega_{d}^{2}}\right)^{3/2}\cos\left(\tan^{-1}\left(\frac{\epsilon}{\omega_{d}}\right)\right)\\\
&=\frac{\kappa}{4}\sqrt{1+\frac{\epsilon^{2}}{\omega_{d}^{2}}}.\end{split}$
(21)

In the specific case $\kappa=4\gamma$ the slowest time scale is set by the
cavity losses, and so $\lambda=-\gamma/2$. It is worth noticing that a Ohmic
spectral density increases too slowly to compensate the effect of the
vanishing matrix element of the dipole transition rate. As a result we would
have that $(\Gamma_{\rm LR}^{\rm dip})^{\rm
Ohm}=\kappa/(4\sqrt{1+\epsilon^{2}/\omega_{d}^{2}})$, and so the overall
relaxation rate would decrease as a function of $\epsilon$ at very weak
coupling $g/\omega_{c}\simeq 0$ to then increase again at slightly larger
coupling.

In the USC regime, $g/\omega_{c}\gg 1$, for small dipole asymmetry
$\epsilon\simeq 0$, the Liouvillian gap goes to zero monothonically with an
exponential behaviour $\lambda\sim-\exp[-g/\omega_{c}]$, as is clearly visible
from Fig. 2(b). Increasing the dipole asymmetry, $\epsilon$, we observe the
emergence of lobes where the Liouvillian gap approaches zero $\lambda\sim 0$,
separated by a narrow region where relaxation is partially restored and
$\lambda\sim-\gamma/2$. This is shown in Fig. 2(c), where we fixed
$g/\omega_{c}=3$ and we plot $\lambda$ as a function of $\epsilon$. Quite
surprisingly, these narrow gaps between the lobes appear only when
$\epsilon\simeq\omega_{c}\times k$, where $k=1,2,3\ldots$ is an integer
number.

###  III.2 Relaxation breakdown in the USC regime

The thermalization slow-down emerging from the spectral analysis of the
Liouvillian $\mathcal{L}$ can be understood as a combination between the
spectral properties and the dressing of the jump operators in Eq. (17) of the
Rabi model in the USC regime Rossatto  _et al._ (2017); Beaudoin  _et al._
(2011). Here we analyze in detail the symmetric case, when $\epsilon=0$, which
will provide the basic tools to understand the whole phase diagram of Fig.
2(a).

We start by transforming the original Rabi Hamiltonian through the unitary
transformation $U_{\rm
pol}=\exp\left[g/\omega_{c}(a-a^{{\dagger}})s_{x}\right]$, and obtaining the
_Rabi polaron_ Hamiltonian ($\hbar=1$)

$\tilde{H}_{\rm Rabi}=\omega_{c}a^{{\dagger}}a+\epsilon
s_{x}+\frac{\omega_{d}}{2}\left[\mathcal{D}(g/\omega_{c})\tilde{s}_{+}+\mathcal{D}^{{\dagger}}(g/\omega_{c})\tilde{s}_{-}\right].$
(22)

Here $\tilde{s}_{\pm}=s_{z}\pm is_{y}$ are the raising/lowering operators
along the $s_{x}$-axis, while
$\mathcal{D}(g/\omega_{c})=\exp\left[g/\omega_{c}(a-a^{{\dagger}})\right]$ is
the usual displacement operator.

Since both cavity and dipole dissipative operators are unaffected by the
polaron transformation $U_{\rm pol}\Phi U_{\rm pol}^{{\dagger}}=\Phi,~{}U_{\rm
pol}s_{x}U_{\rm pol}^{{\dagger}}=s_{x}$, the general master equation in Eq.
(19) is still valid, with the only difference that the eigenstates
$|n\rangle,|m\rangle$ appearing in the jump operators in Eq. (17) are now
replaced with the eigenstates of the Rabi polaron Hamiltonian in Eq. (22).

As reported in Irish (2007) and detailed in Appendix C, the polaron Rabi
Hamiltonian supports a generalized rotating-wave approximation (gRWA) and thus
follows the structure of the Jaynes-Cummings model, with the approximated
conservation of the polaron excitation number $\hat{N}_{\rm
exc}^{z}=a^{{\dagger}}a+s_{z}$. Its eigenstates are then given by the usual
dressed states

$\begin{split}&|+,n\rangle=\cos\frac{\theta_{n}}{2}|\downarrow,n\rangle+\sin\frac{\theta_{n}}{2}|\uparrow,n-1\rangle,\\\
&|-,n\rangle=-\sin\frac{\theta_{n}}{2}|\downarrow,n\rangle+\cos\frac{\theta_{n}}{2}|\uparrow,n-1\rangle,\end{split}$
(23)

where $\theta_{n}$ is given in Appendix C. The ground-state of the system is
simply the uncoupled vacuum state

$|{\rm GS}\rangle=|\downarrow,0\rangle.$ (24)

In order to appreciate the quality of this approximation, in Fig. 3(a) we
compare the spectrum obtained from the exact diagonalization (solid lines) and
from the gRWA analytical formula reported in Appendix C (yellow dots), from
which is quite clear that the gRWA gives very good results.

Relaxation can be then understood from the dressed state perspective Claude
Cohen-Tannoudji (1998) and in Appendix D we explicitly compute the transition
rate in Eq. (20).

Figure 3: (a) Spectrum of the Rabi model as a function of the light-matter
coupling $g$ at fixed $\epsilon=0$. The solid lines are the result of full
diagonalization, and the color red/blue are only meant to match the color-code
in (c). The yellow dots are given by the analytic Eq. (72) in Appendix C. (b)
$\cos\theta_{n}/2$, $\sin\theta_{n}/2$ given by Eq. (74) for each $n$ block as
a function of the light-matter coupling $g$. (c) Scheme of the relaxation
mechanism. The cavity relaxes jumping mainly between $++$ or $--$ dressed
states, while for the dipole is mainly between $+-$ states. In the USC limit
the dipole does not relax anymore. Parameters: $\epsilon=0$,
$\omega_{c}=\omega_{d}$.

From the explicit expression for the Hopfield coefficients $\sin,\cos$
(present in the Appendix C), we found that in USC

$\begin{split}&\lim_{g/\omega_{c}\rightarrow\infty}\cos\frac{\theta_{n}}{2}=1\\\
&\lim_{g/\omega_{c}\rightarrow\infty}\sin\frac{\theta_{n}}{2}=0,\end{split}$
(25)

which is clearly shown in Fig. 3(b) where we plot the Hopfield coefficient
analytically computed through the gRWA for a few lowest eigenstates. Using
this observation together with the matrix element computed in Appendix D we
can build the transition rates in Eq. (20), arriving to the conclusion that
the only non-negligible transitions are

$\begin{split}&\lim_{g/\omega_{c}\rightarrow\infty}\Gamma_{(+,n)(+,n-1)}^{c}=\gamma\frac{\omega_{+,n}-\omega_{+,n-1}}{\omega_{c}}\approx\gamma\\\
&\lim_{g/\omega_{c}\rightarrow\infty}\Gamma_{(-,n)(-,n-1)}^{c}=\gamma\frac{\omega_{-,n}-\omega_{-,n-1}}{\omega_{c}}\approx\gamma\\\
&\lim_{g/\omega_{c}\rightarrow\infty}\Gamma_{(-,n)(+,n-1)}^{\rm\,dip}=\kappa\left(\frac{\omega_{-,n}-\omega_{+,n-1}}{\omega_{c}}\right)^{3}\approx
0.\end{split}$ (26)

The fact that the dipole transition rate goes to zero
$\Gamma_{(-,n)(+,n-1)}^{\rm\,dip}\approx 0$ follows from the approximate
degeneracy of the states $|-,n\rangle,|+,n-1\rangle$ in the USC limit, for
which $\omega_{-,n}-\omega_{+,n-1}\approx 0$, while
$\omega_{+,n}-\omega_{+,n-1}\approx\omega_{-,n}-\omega_{-,n-1}\approx\omega_{c}$.

Here we realize that the USC relaxation breakdown observed in Fig. 2(a-b)
affects the dipole only, while the cavity transition rates return to their
bare uncoupled values when the USC regime is reached. The schematic
representation of the remaining relaxation channels is shown in Fig. 3(c).

This result can be directly interpreted from a polaronic perspective: the USC
cavity vacuum heavily dresses the dipole with virtual photons, which inhibit
its ability to tunnel from one side to the other of its double well potential.
Because of the radiative nature (but in the Ohmic case as well) of its
dissipation mechanism, the dipole can loose energy only moving between the two
wells (i.e. tunneling), and the faster it moves the stronger it dissipates. In
this regime of heavy virtual-photon dressing tunneling becomes extremely slow
and so the dipole’s rate to release energy in the bath.

##  IV USC multi-photon resonant tunneling

In this section we are going to explore more in detail the nature of the gaps
between the relaxation-breakdown lobes in Fig. 2(a). In these narrow regions
the system can relax as is almost unaffected by the USC suppression of
tunneling described in the previous section. However, if we artificially
remove the cavity dissipation, $\gamma=0$, we see that these narrow gaps
disappear. This suggests that the suppression of tunneling described above is
still present for $\epsilon\neq 0$, but a new resonant mechanism appears,
allowing the dipole to tunnel again by exchanging photons with the cavity.
This effect is the cavity analogous of resonant tunneling in electronic setups
interacting with vibrational degrees of freedom, and thus establishing a
connection between USC cavity QED and Franck-Condon physics in molecular-
electronic setups Koch  _et al._ (2006); Cui  _et al._ (2015); Vdovin  _et
al._ (2016).

###  IV.1 Diagonalization of the asymmetric Rabi Hamiltonian

Figure 4: Spectrum of the Rabi model as a function of the dipole asymmetry
$\epsilon$ for various $g/\omega_{c}=0.1,1,2.5,3.5$ light-matter couplings.
The solid lines are the result of exact diagonalization while the yellow dot
are given by the analitical formula in Eq. (31). Parameters
$\omega_{d}=\omega_{c}$.

As in the symmetric case $\epsilon=0$, also the asymmetric Rabi model,
$\epsilon\neq 0$ is approximately block-diagonal, as a consequence of the
general form of the displacement operators. However here the situation is more
complicated and we cannot find a unique formula that fits the whole spectrum
for every $(\epsilon,g)$, but we can only have analytic expressions valid near
to each resonance.

We start by noticing that Eq. (22) is written in a form that calls for the
gRWA, provided that the system has an _asymmetric resonance_
$\epsilon\simeq\omega_{c}\times k$, with $k=1,2,\ldots$. Differently from the
usual Jaynes-Cummings model, and the gRWA of the symmetric Rabi model, here we
need to take the dipole basis as an eigenstate of $s_{x}$. Moreover,
considering higher resonances at
$\epsilon=\omega_{c},2\omega_{c},3\omega_{c}\ldots$ is well motivated by the
fact that the displacement operator contains all power of
creation/annihilation operators, giving access to multi-photon processes with
higher frequencies. This is indeed well visible considering the normal-order
expansion Cahill and Glauber (1969)

$\mathcal{D}(x)=e^{-x^{2}/2}\sum_{n,m=0}\frac{(xa^{{\dagger}})^{n}}{n!}\frac{(-xa)^{m}}{m!}.$
(27)

From this expression is also clear that the non-linear interaction term in the
polaron Hamiltonian in Eq. (22) is exponentially suppressed by the factor
$\sim\omega_{d}e^{-g^{2}/(2\omega_{c})}$. As a consequence, when

$\epsilon,\omega_{c}>\omega_{d}e^{-g^{2}/(2\omega_{c})}$ (28)

the polaron light-matter interaction becomes perturbative, and we can adopt
the gRWA. Notice that this correspond to keep only the terms $n<m$ with
$\omega_{c}(m-n)\simeq\epsilon$ in Eq. (27), so, even if the interaction is
perturbative, is still multi-photon and thus highly non-linear.

The asymmetric polaron Rabi Hamiltonian can then be approximately diagonalized
around each $k$-resonance by projecting it on the states
$\\{{|\leftarrow,n\rangle,|\rightarrow,n-k\rangle\\}}$ and the ground-state is
simply given by $|{\rm GS}\rangle\approx|\leftarrow,0\rangle$.

The Hamiltonian can be then expressed succinctly in a matrix form, as the sum
of $2\times 2$ blocks

$\begin{split}&\tilde{H}_{\rm
Rabi}^{k}\approx\sum_{n=1}^{\infty}\frac{\omega_{c}k-\epsilon}{2}\sigma_{x}^{(n,k)}+\frac{\omega_{d}}{2}\mathcal{D}_{n\,n-k}\,\sigma_{z}^{(n,k)}\\\
&+\frac{2\omega_{c}n-\omega_{c}k}{2}\mathds{1}_{(n,k)},\end{split}$ (29)

where $\sigma_{x,y,z}^{(n,k)}$ are the Pauli matrices for each $n=1,2,\ldots$
block for the $k$-resonance, while

$\mathcal{D}_{n\,n-k}=\frac{g^{k}}{\omega_{c}^{k}}e^{-\frac{g^{2}}{2\omega_{c}^{2}}}L_{n-k}^{(k)}\left(g^{2}/\omega_{c}^{2}\right)\sqrt{\frac{(n-k)!}{n!}}$
(30)

is the $n,n-k$ matrix element of the displacement operator Cahill and Glauber
(1969). Here $L_{m}^{(l)}(x)$ is the special Laguerre polynomials. The excited
eigenfrequencies are then given by

$\begin{split}&\omega_{k,n,\pm}^{\rm
R}=\omega_{c}\left(n-\frac{k}{2}\right)\pm\frac{1}{2}\sqrt{\left(\omega_{c}k-\epsilon\right)^{2}+\omega_{d}^{2}\mathcal{D}_{n\,n-k}^{2}}.\end{split}$
(31)

Since the displacement operator has diagonal matrix element different from
zero $\mathcal{D}_{nn}\neq 0$, one should consider the dipole basis states
composed by dipole states oriented along $\sim\cos\phi s_{x}+\sin\phi s_{z}$,
with a certain angle $\phi$ given by $\mathcal{D}_{nn}$. Including these
corrections makes the analytical formula in general quite complicated, having
a simple expression only for the ground-state, which is

$\omega^{\rm
R}_{0,0}=-\frac{\sqrt{\epsilon^{2}+\omega_{d}^{2}e^{-g^{2}/\omega_{c}^{2}}}}{2}.$
(32)

However, in the USC regime $\phi\sim 0$ is a small angle and we can thus
neglect it, proceeding with the simple $s_{x}$ picture developed above.

In Fig. 4 we compare the real spectrum to the one obtained from the gRWA at
each resonant point. In the USC limit, when $g/\omega_{c}\gg 1$ the agreement
is very good.

The eigenstates are now given in terms of a multi-photon version of the
$s_{x}$-polarized Jaynes-Cummings dressed states, fully characterized by the
Hopfield coefficients $\cos\theta_{(k,n)}/2,\sin\theta_{(k,n)}/2$,
generalizing the symmetric case in Appendix C. When the resonance condition
$\epsilon=\omega_{c}\times k$ is satiesfied, the system eigenstates become

$\begin{split}&|+_{(k,n)}\rangle=\frac{1}{\sqrt{2}}\left(|\leftarrow,n\rangle+|\rightarrow,n-k\rangle\right)\\\
&|-_{(k,n)}\rangle=\frac{1}{\sqrt{2}}\left(|\leftarrow,n\rangle-|\rightarrow,n-k\rangle\right),\end{split}$
(33)

and the ground state is $|{\rm GS}\rangle=|\leftarrow,0\rangle$.

Repeating the analysis on the matrix elements in Appendix D, we realize that
the $s_{x}$ operator can only connect dressed states of the same $(k,n)$
block, for which the only non-diagonal non-zero matrix element is

$\braket{+_{(k,n)}}{s_{x}}{-_{(k,n)}}=\cos\frac{\theta_{(k,n)}}{2}\sin\frac{\theta_{(k,n)}}{2}.$
(34)

Each block is disconnected by the others and the ground-state is disconnected
from all other states. Therefore relaxation toward equilibrium is still
suppressed from the USC also when $\epsilon\neq 0$.

On contrary, the cavity is still able to efficiently dissipate. So, when
hitting a $k$-resonance, also the dipole can lose energy by exchanging
$k$-photons with the cavity, which are consequently flushed out. This resonant
tunneling effect provides a relaxation channel for the dipole, as depicted in
Fig. 5(a) and give the proper explanation for the gaps between the lobes
observed in Fig. 2(a).

We conclude this subsection by highlighting that: in the polaron frame, the
USC open dynamics is mainly given by a polaronic version of the standard text-
book dressed state master equation dynamics Claude Cohen-Tannoudji (1998).

###  IV.2 Multi-photon oscillations and cavity-mediated relaxation

Figure 5: (a) Schematic view of the resonant tunnel mechanism. In the USC
regime, the dipole can switch well by exchanging $k$-photons with the cavity.
(b) Rabi oscillations data collapse. Each curve is labelled by its resonant
index $k$ and represents
$\braket{\tilde{s}_{x}}=e^{k\gamma/2t}(\braket{s_{x}}+1/2)-1/2$. For each
$k$-curve the time is normalized on its respective $k$-Rabi frequency,
$\Omega_{(k,k)}/(2\pi)$. In this way is clear how our analytical description
fits very well the full numerics. Parameters: $\omega_{d}=\omega_{c}$,
$g/\omega_{c}=3$, $\gamma=\kappa/4=0.002\omega_{c}$.

Here we illustrate how the polaronic dressed state picture emerges clearly in
the full open dynamics. As a striking example we show that the system
undergoes to damped Rabi oscillations, as in traditional cavity QED systems
described by the Jaynes-Cummings model. However here, depending from the
resonance condition, the Rabi oscillations involve multiple photons Garziano
_et al._ (2015); Ma and Law (2015) and must be interpreted as tunneling
oscillations for the dipole.

From the block-Hamiltonian in Eq. (29) we can derive the Rabi frequency of the
$k$-resonance multi-photon Rabi oscillations reading

$\Omega_{(k,n)}=\omega_{d}\frac{g^{k}}{\omega_{c}^{k}}e^{-\frac{g^{2}}{2\omega_{c}^{2}}}L_{n-k}^{(k)}\left(g^{2}/\omega_{c}^{2}\right)\sqrt{\frac{(n-k)!}{n!}},$
(35)

where $n\geq k$ is the total number of photons involved.

Differently from usual Rabi oscillations in cavity QED, here the dipole
oscillates between the right and left states of its asymmetric double well
potential, for which we can call them _tunneling oscillations_. The relevant
quantity to follow is then $\braket{s_{x}}(t)$ (on contrary to traditional
Rabi oscillations, visible looking at $\braket{s_{z}}(t)$, in standard
notation).

We then numerically simulate $\braket{s_{x}}(t)$ starting from the initial
state $|\psi_{0}\rangle=|\rightarrow,0\rangle$ (in the polaron frame). When
$\gamma<\Omega_{(k,k)}$, we observe a very good fit on the curve

$\braket{s_{x}}(t)\approx
e^{-\frac{k\gamma}{2}t}\frac{\cos\left[\Omega_{(k,k)}t\right]+1}{2}-\frac{1}{2}.$
(36)

Notice that the overall decay rate is given by $\sim k\times\gamma/2$, with a
factor $k$. This takes into account that the photon decay increase linearly
with the number of photons involved, which, in this case is properly $k$. In
Fig. 5(b) we show the Rabi oscillations data collapse for various resonant
values $\epsilon=\omega_{c},2\omega_{c},3\omega_{c},4\omega_{c}$. For each $k$
the curve is plotted against its normalized time
$\tilde{t}=\Omega_{(k,k)}t/(2\pi)$ and is normalized to remove the exponential
decay accordingly to
$\braket{\tilde{s}_{x}}=e^{k\gamma/2t}(\braket{s_{x}}+1/2)-1/2$.

###  IV.3 Response functions and higher-order processes

Figure 6: Current transmission $|\mathcal{T}(\omega)|$ as a function of the
dipole asymmetry $\epsilon$ and the probe frequency $\omega$, for various
$g/\omega_{c}=0.1,0.5,2,2.5$ light-matter couplings. Parameters:
$\omega_{d}=\omega_{c}$, $k_{b}T=0.2\hbar\omega_{c}$,
$Q=\omega_{c}/\gamma=10^{2}$.

Here we investigate how the spectral features analyzed so far manifest
themselves through standard transmission measurements.

We start considering a weak probe current entering in the LC-circuit and we
look for the transmitted current. With the help of linear response theory (see
Appendix E), the current response is mainly given by the cavity structure
factor

$\mathcal{S}_{c}(\omega)=\frac{\hbar
Z_{LC}}{2}\sum_{n,m}\frac{e^{-\hbar\omega_{n}/(k_{b}T)}}{\mathcal{Z}}\left|\braket{n}{a-a^{{\dagger}}}{m}\right|^{2}\delta(\omega-\omega_{mn}).$
(37)

Here $\mathcal{Z}=\sum_{n}e^{-\hbar\omega_{n}/(k_{b}T)}$ is the thermal
equilibrium partition function of the system. From this quantity is possible
to define the system impedance

$Z_{\rm sys}(\omega)=-\frac{i\omega\mathcal{S}_{c}(\omega)}{\hbar},$ (38)

and consequently the current transmission function

$\frac{I_{\rm out}}{I_{\rm
in}}=\mathcal{T}(\omega)=\frac{Q^{-1}}{Q^{-1}+Z_{\rm LC}/Z_{\rm
sys}(\omega)}.$ (39)

Here $Q=\omega_{c}/\gamma$ is the LC cavity quality factor.

In Fig. 6 we show the current transmission $|\mathcal{T}(\omega)|$ as a
function of the dipole asymmetry $\epsilon$ and the probe frequency $\omega$.
To mimic experimental conditions, we consider a fixed temperature
$k_{b}T\simeq 0.2\hbar\omega_{c}$.

For coupling strength up to $g/\omega_{c}\lesssim 0.5$ the transmission
spectrum exhibits the usual Jaynes-Cummings polaritonic (or dressed state)
behaviour. At $\epsilon\neq 0$ the response mainly follows the bare LC circuit
response, while the maximum hybridization is at $\epsilon=0$, with maximum
Rabi splitting between the upper and lower dressed state branches. This is
well visible in the first panel of Fig. 6.

For increasing coupling strength, as in the second panel of Fig. 6, where
$g/\omega_{c}=1$, the upper branch of the transmission spectrum starts to
vanish exactly at $\epsilon=0$, signalling that we are entering the USC
regime.

At even larger couplings the transmission is drastically changed. This is well
visible from the two lower panels in Fig. 6, where $g/\omega_{c}=2.5,3$. In
the region around $\epsilon=0$ the transmission becomes much smaller, and the
two branches related to the Jaynes-Cummings dressed states are gone. Instead
the $k=1$ avoided crossing due to the resonant tunneling is well visible
around $|\epsilon|/\omega_{c}=1$. The $k>1$ higher resonances, on contrary,
are not well visible, since they are covered by the bare photon resonance
between the lower levels $n<k$. It is worth noticing that the ultrastrong-
coupling spectral features shown here, and in particular the $k=1$ resonance,
are already visible in recent experiments with superconducting circuits
Niemczyk  _et al._ (2010); Yoshihara  _et al._ (2017).

Another interesting quantity to probe the spectrum of the system is provided
by the dipole structure factor

$\mathcal{S}_{\rm dip}(\omega)=2\hbar Z_{\rm
dip}\sum_{n,m}\frac{e^{-\hbar\omega_{n}/(k_{b}T)}}{\mathcal{Z}}\left|\braket{n}{s_{x}}{m}\right|^{2}\delta(\omega-\omega_{mn}).$
(40)

Here we introduced the dipole impedance as $Z_{\rm dip}=\hbar/q^{2}f_{10}$,
and $f_{10}=2m\omega_{0}|x_{10}|^{2}/\hbar$ is the oscillator strength of the
two-level dipole transition.

From the linear response theory perspective, $\mathcal{S}_{\rm dip}(\omega)$
quantifies the dipole radiation response to a direct drive of the dipole.
Because of the consideration done in Sec. IV.1, is clear that this quantity is
strongly suppressed in the USC regime, at low temperature. If we only stick to
the dressed state picture we should observe vanishing transitions for
$T\rightarrow 0$, due to the fact that the dipole matrix element between the
ground-state and the first block is zero

$\braket{\pm_{(k,1)}}{s_{x}}{\rm GS}=0.$ (41)

This is thus a perfect framework where to look for transitions beyond the
dressed state picture.

Figure 7: Dipole radiation impedance $|Z_{\rm rad}(\omega)|$ as a function of
the dipole asymmetry $\epsilon$ and the probe frequency $\omega$ for various
$g/\omega_{c}=2,2.5$ light-matter couplings. Parameters:
$\omega_{d}=\omega_{c}$, $k_{b}T=0.5\hbar\omega_{c}$. In this plot we assumed
a linewidth $\gamma_{\mathcal{S}_{\rm dip}}=0.05\omega_{c}$.

Similarly to the system circuit impedance we can define a dipole radiation
impedance as

$Z_{\rm rad}(\omega)=-\frac{i\omega\mathcal{S}_{\rm dip}(\omega)}{\hbar}.$
(42)

In Fig. 7 we show the radiation impedance $Z_{\rm rad}(\omega)$ in logscale,
as a function of the dipole asymmetry $\epsilon$ and the probe frequency
$\omega$ (with an artificial linewidth $\gamma_{\mathcal{S}_{\rm dip}}$ to
smear out the delta function in Eq. (40)). On contrary to the previous case of
the cavity response, at frequencies $\omega\sim\omega_{c}$ the dipole response
is strongly suppressed in favour of higher frequencies transitions that emerge
with a diamond-like pattern.

The dipole matrix elements giving the amplitude for these transitions are much
weaker than the cavity-current matrix elements for the $k$-resonant
transitions between dressed states, since are given by beyond gRWA
corrections. They are the USC cavity equivalent of the vibronic transitions
responsible of the Coulomb diamond structure in the Franck-Condon blockade
voltage/current characteristic Koch  _et al._ (2006).

##  V Cascaded relaxation in multi-well dipole

Finally we comment on the possibility to extend our results to the case of a
multi-well dipole. The relaxation dynamics of this system is particularly
interesting because it can be interpreted as a prototype of a transport
problem through an extended system. Despite that a full coverage of cavity-
modified relaxation or transport in an extended system is well beyond the
scope of this paper, we can still use the concepts developed above to give an
initial intuition which sets the basis for future investigations.

###  V.1 The extended Dicke model

We model the multiple-well dipole generalizing the two-level approximation to
$(N+1)$-level, where each level represent a potential well. In this way, the
dipole is simply described by spin-$N/2$ operators, $S_{x,y,z}$ that
generalizes the spin description of Sec. II.1. In particular the eigenvalues
of $S_{x}$, $|m_{x}\rangle$, are interpreted as localized states in the
$m_{x}$-th well, while the $S_{z}$ operator creates some tunnelling between
them.

While the dissipations, and the master equation, are derived in the same way
as before, just replacing $s_{x,y,z}\mapsto S_{x,y,z}$ everywhere, the light-
matter Hamiltonian is no more given by the Rabi model of Eq. (10). Indeed,
when taking the two-level approximation of Eq. (9) we had discarded the
$x^{2}$-term, which is only a constant within the two-level subspace,
$x^{2}\approx 4x_{10}^{2}s_{x}^{2}=4x_{10}^{2}\mathds{1}$. For a multi-level
dipole, described by a spin-$N/2$ system if $N>1$ we have that
$S_{x}^{2}\neq\mathds{1}/4$, and thus the correct cavity QED Hamiltonian
within the $(N+1)$-level subspace is given by the so-called _extended Dicke
model_ (EDM) Jaako  _et al._ (2016); De Bernardis  _et al._ (2018a)

$H_{\rm
EDM}=\omega_{c}a^{{\dagger}}a+\omega_{d}S_{z}+g\left(a+a^{{\dagger}}\right)S_{x}+\frac{g^{2}}{\omega_{c}}S_{x}^{2}+\epsilon
S_{x}.$ (43)

Performing the polaron transformation in the same way as for the Rabi model,
we arrive to the polaron (EDM) De Bernardis  _et al._ (2018a):

$\tilde{H}_{\rm EDM}=\omega_{c}a^{{\dagger}}a+\epsilon
S_{x}+\frac{\omega_{d}}{2}\left[\mathcal{D}(g/\omega_{c})\tilde{S}_{+}+\mathcal{D}^{{\dagger}}(g/\omega_{c})\tilde{S}_{-}\right],$
(44)

where, again, $\tilde{S}_{-}=S_{z}-iS_{y}$. In the USC regime
$g\gg\omega_{c}$, and for non-negligible asymmetry $\epsilon\neq 0$, in the
limit of large spin, $N\gg 1$, we can use the Holstein-Primakoff approximation
Holstein and Primakoff (1940) on the $S_{x}$-direction, for which
$S_{x}\approx-N/2+b^{{\dagger}}b$, and $\tilde{S}_{-}\approx\sqrt{N}b$. The
polaron EDM can then be approximated by

$\begin{split}\tilde{H}_{\rm EDM}\approx\tilde{H}_{\rm EDM}^{\rm
HP}&=\omega_{c}a^{{\dagger}}a+\epsilon b^{{\dagger}}b+\\\
&+\frac{\omega_{d}\sqrt{N}}{2}\left[\mathcal{D}(g/\omega_{c})b^{{\dagger}}+\mathcal{D}^{{\dagger}}(g/\omega_{c})b\right].\end{split}$
(45)

###  V.2 Relaxation dynamics of the EDM

The considerations done for the Rabi model in Sec. IV.1 are still valid, in
particular regarding the possibility to discard the counter-rotating terms in
the displacement operators and the suppression of tunneling. It is then clear
that one can use the same dressed state approach to diagonalize the polaron
EDM as well. In particular, considering the relaxation from the initial state
$|0,m\rangle=(b^{{\dagger}})^{m}/\sqrt{m!}|0,0\rangle$ the resonant tunnelling
effect gives rise to a cavity-mediated cascaded dynamics to the ground-state,
where, depending from the resonance condition $\epsilon=k\times\omega_{c}$,
$n_{\rm ph}\approx k\times m$ photons are released.

To have a more quantitative understanding we consider the limit of strong
cavity dissipations with respect to the k-resonance splitting,
$\gamma\gg\Omega_{(k,k)}$. In this regime we can adiabaticaly eliminate the
cavity in favor of an effective master equation for the dipole only Gardiner
and Zoller (2015). Following the previous analysis on cavity and dipole
transition rates, we completely neglect the dipole dissipations, while we take
as a jump operator of the cavity its bare annihilation operator $c=a$. Again,
this is well motivated by the analysis performed above. Using the approximated
form of the EDM in Eq. (45) and assuming that at each time the total density
matrix of the system is $\rho(t)\approx\rho_{d}(t)\otimes\rho^{\rm th}_{c}$
(here $\rho^{\rm th}_{c}$ is the thermal density matrix for the bare cavity at
temperature $T$) we have that

$\begin{split}&\partial_{t}\rho_{d}=-i\left[\epsilon
b^{{\dagger}}b,\rho_{d}\right]+\frac{\Gamma_{T}(\epsilon)}{2}\left(2b\rho_{d}b^{{\dagger}}-\left[b^{{\dagger}}b,\rho_{d}\right]_{+}\right)\\\
&+\frac{\Gamma_{T}(-\epsilon)}{2}\left(2b^{{\dagger}}\rho_{d}b-\left[bb^{{\dagger}},\rho_{d}\right]_{+}\right).\end{split}$
(46)

Similarly to non-linear optomechanics setups Rabl (2011); Minoguchi  _et al._
(2019), the cooling and heating rates are given by

$\begin{split}&\Gamma_{T}(\omega)=\frac{\omega_{d}^{2}N}{2}\times\\\
&\times{\rm
Re}\left[\int_{0}^{\infty}dt\left(\braket{\mathcal{D}\left(t,x\right)\mathcal{D}^{{\dagger}}\left(x\right)}-\braket{\mathcal{D}(x)}^{2}\right)e^{i\omega
t}\right],\end{split}$ (47)

where $H_{c}=\omega_{c}a^{{\dagger}}a$, $x=g/\omega_{c}$ and
$\mathcal{D}(t,x)=e^{iH_{c}t}\mathcal{D}\left(x\right)e^{-iH_{c}t}$. Since
that the average $\braket{\cdot}$ is intended over the cavity thermal state
$\rho_{c}^{\rm th}$, we can explicitly compute this quantity Rabl (2011);
Minoguchi  _et al._ (2019); Pilar  _et al._ (2020)

$\begin{split}&\Gamma_{T}(\omega)=\frac{\omega_{d}^{2}N}{\gamma}e^{-x^{2}\left(1+2N_{T}(\omega_{c})\right)}\times\\\
&\times\sum_{q,r\neq
0}\frac{x^{2r}N_{T}^{r}(\omega_{c})}{r!}\frac{x^{2q}(1+N_{T}(\omega_{c}))^{q}}{q!}\frac{\gamma^{2}/4}{\left(\omega-\omega_{c}(q-r)\right)^{2}+\frac{\gamma^{2}}{4}}.\end{split}$
(48)

Here $N_{T}(\omega_{c})=1/(e^{\hbar\omega_{c}/(k_{b}T)}-1)$ is the cavity
thermal population. From this expression it is particularly evident the multi-
photon character of this cavity assisted relaxation mechanism, where the
dipole can relax by emitting $q$-photons into the cavity, and, at the same
time, can be re-excited by absorbing $r$-photons from the cavity (if the
temperature is non-zero $T>0$).

Figure 8: Total relaxation rate $\Gamma^{\rm tot}_{T}(\omega)$ as a function
of frequency $\omega$ for various $g/\omega_{c}=0.1,0.5,1,2$ coupling
strengths (normalized over $\Gamma_{d}=\omega_{d}^{2}N/\gamma$). The blue line
corresponds to $T=0$, while the red line is for $k_{b}T=2\,\hbar\omega_{c}$.
Parameters: $\gamma/\omega_{c}=0.1$.

As in the standard theory of laser cooling the total relaxation rate is given
by $\Gamma^{\rm tot}_{T}=\Gamma_{T}(\epsilon)-\Gamma_{T}(-\epsilon)$. At $T=0$
and close to resonance $\epsilon\simeq\omega_{c}\times k$, the total
relaxation rate is approximately given by $\Gamma_{T=0}^{\rm
tot}\approx\Omega_{(k,k)}^{2}N/\gamma$, which is the USC version of the
Purcell effect. In Fig. 8 we show some examples of the total relaxation rate
for $T=0$ and for $T>0$ at different coupling strengths. Interestingly larger
temperature may help in activating the higher $k$-resonances even in the non-
USC regime, $g/\omega_{c}\ll 1$, resembling the behaviour of optomechanical
laser cooling setups Rabl (2010). We highlight the fact that at weak coupling
this description does not hold, since it is based on the assumption that
dipole tunneling is suppressed by the USC regime. However we found interesting
to show the total relaxation rate $\Gamma^{\rm tot}_{T}$ also in this regime
for completness.

As anticipated in the beginning of this section, the relaxation dynamics of
this multi-well setup can be seen as a way to study how the incoherent
transport is modified by the cavity. Staying at the single particle level we
can interpreted an excitation produced by the dipole operator $b^{{\dagger}}$
as the particle moving one well up in energy, and the dipole ground-state as
the state where only the lowest energy well is occupied. Following this line
of thoughts, we can say that the system has good transport properties if it
rapidly thermalizes sufficiently close to its ground-state. Since that the
saturation number of the steady state of Eq. (46),
$\braket{b^{{\dagger}}b}_{\rm ss}=\mathcal{N}_{0}$, around each $k$-resonance
$\epsilon=\omega_{c}\times k$ is the dipole thermal occupation

$\mathcal{N}_{0}=\frac{\Gamma_{T}(-\epsilon)}{\Gamma_{T}(\epsilon)-\Gamma_{T}(-\epsilon)}\approx
N_{T}(\omega_{c}\times k),$ (49)

we also need that the temperature $T$ is small enough so thermal photons
cannot push the particle (the dipole excitation) on an upper energy level.

##  VI Conclusion

In conclusion, we studied the relaxation properties of a simple (but
paradigmatic) cavity QED setup in the ultrastrong coupling regime where the
cavity is provided by an LC resonant circuit and the matter is given by an
asymmetric double-well dipole inside the capacitor of the LC circuit. By
considering the LC circuit coupled to an Ohmic transmission line we introduced
current dissipations for the cavity, while considering the coupling to
radiating modes we introduced dissipations for the dipole. In this way the
system is provided with two among the simplest and most intuitive relaxation
channels. After having defined the basic framework, we derived a thermalizing
master equation, valid at arbitrary light-matter coupling strengths and
arbitrary dipole asymmetry. From the Liouvillian gap we obtained the longest
relaxation rate of the system, that we can also consider its thermalization
rate. From this quantity emerges clearly that the effect of the USC is to
slow-down the system’s thermalization by an exponential suppression of the
Liouvillian gap. However, for special values of the dipole asymmetry the
standard relaxation is restored and the system can thermalize accordingly to
its bare relaxation rates.

In order to get more insights on this behaviour of the Liouvillian gap, we
employed a generalized rotating-wave approximation (gRWA) Irish (2007), valid
in the so-called polaron frame. Within this approximation we showed that is
possible to diagonalize the asymmetric Rabi model, obtaining analytical
expressions for its spectrum and eigenstates, valid in the USC regime.

Specifically it emerges that in the polaron frame the eigenstates are given by
a polaronic multi-photon version of the usual Jaynes-Cummings dressed states.
Within this simple picture we were able to compute the relaxation rates in the
USC regime analytically, explicitly showing the exponential slow-down of
thermalization due to an effective suppression of the dipole tunnelling
dynamics, while the cavity, remaining effectively uncoupled from the dipole,
can still efficiently relax.

However when the dipole asymmetry is resonant with the cavity, the dipole
dynamics can be revitalized. Here the dipole dynamics is dominated by cavity
assisted tunnelling, where the dipole can resonantly tunnel from one well to
the other by releasing multiple photons. Since photons can then relax out of
the cavity, this process gives an effective relaxation channel also for the
dipole.

After establishing a connection with the Franck-Condon physics of electronic
transport through a molecular dotKoch  _et al._ (2006); Valmorra  _et al._
(2021), we commented on the possibility of extending this non-linear resonant
process to a multi-well dipole. A simple toy model to describe this situation
is provided by the extended Dicke model, introduced originally to study a
multiple qubit ultrastrongly coupled to a single LC cavity Jaako  _et al._
(2016); De Bernardis  _et al._ (2018a). From this setup is clear that the USC
resonant tunnel dynamics can affect also a multi-well system, giving rise to a
resonant cascaded multi-photon process. We argue that this cascaded effect
could be observed in cavity modified transport experiments with multiple
electronic quantum dots, or in supercunding circuit devices with only minor
modifications of the already existing platforms Yoshihara  _et al._ (2017);
Valmorra  _et al._ (2021). This suggests that these findings could thus
provide an interesting playground to study an implementation for cascaded-
laser electronic devices operating in the USC regime in GHz or THz range.

###### Acknowledgements.

We thanks Gianluca Rastelli, Iacopo Carusotto, Peter Rabl, Alberto Biella,
Fabrizio Minganti, Alberto Nardin and Gian Marcello Andolina for very helpful
and insightful discussions. We acknowledge financial support from the
Provincia Autonoma di Trento from the Q@TN initiative.

##  Appendix A Linear damping

We consider a generic system, described by the hamiltonian $H_{\rm sys}$,
coupled to a bath of harmonic oscillators (which may represent the
electromagnetic field outside of a cavity, or a resistance in a circuit):

$H=H_{\rm
sys}+\sum_{k}\left[\frac{P_{k}^{2}}{2m_{k}}+\frac{1}{2}m_{k}\omega_{k}^{2}\left(Y_{k}-\frac{c_{k}}{m_{k}\omega_{k}^{2}}X\right)^{2}\right].$
(50)

The equations of motion for a generic system operator $A$ are given by

$\begin{split}&\partial_{t}A=-i\left[A,H_{\rm
sys}\right]+i\sum_{k}\frac{c_{k}}{2}\left(Y_{k}\left[A,X\right]+\left[A,X\right]Y_{k}\right)+\\\
&-i\sum_{k}\frac{c_{k}^{2}}{2m_{k}\omega_{k}^{2}}\left[A,X^{2}\right]\end{split}$
(51)
$\partial_{t}Y_{k}=\frac{P_{k}}{m_{k}},~{}~{}~{}~{}\partial_{t}P_{k}=-m_{k}\omega_{k}^{2}Y_{k}+c_{k}X.$
(52)

The formal solution of the bath’s equations is given by

$Y_{k}=Y_{k}^{homg.}(t)+\frac{c_{k}}{m_{k}\omega_{k}}\int_{t_{0}}^{t}dt^{\prime}\sin(\omega_{k}(t-t^{\prime}))X(t^{\prime}),$
(53)

where

$Y_{k}^{homg.}=Y_{k}(t_{0})\cos(\omega_{k}(t-t_{0}))+\frac{P_{k}(t_{0})}{m_{k}\omega_{k}}\sin(\omega_{k}(t-t_{0})).$
(54)

Plugging back this solution is (51), and integrating by part, we get

$\begin{split}\partial_{t}A&=-i\left[A,H_{\rm
sys}\right]+\frac{i}{2}\left(\xi(t)\left[A,X\right]+\left[A,X\right]\xi(t)\right)\\\
&-\frac{i}{2}\left[\int_{t_{0}}^{t}K(t-t^{\prime})\partial_{t^{\prime}}X(t^{\prime})dt^{\prime}\,,\,\left[A,X\right]\right]_{+}.\end{split}$
(55)

Here $[\cdot,\cdot]_{+}$ is the anti-commutator, and

$\begin{split}\xi(t)&=\sum_{k}c_{k}\left(Y_{k}^{homg.}(t)-\frac{c_{k}}{m_{k}\omega_{k}^{2}}X(t_{0})\cos(\omega_{k}(t-t_{0}))\right)\\\
K(t)&=\sum_{k}\frac{c_{k}^{2}}{m_{k}\omega_{k}^{2}}\cos(\omega_{k}t),\end{split}$
(56)

are, respectively, the _quantum noise term_ and the _dissipative kernel_. We
notice that the last term in (51) is exactly cancelled by the term
proportional to $C(t)$ coming out by the integration by part. From the
fluctuations-dissipation theorem we obtain the specific value of the quantum
noise correlator Breuer  _et al._ (2002). In the high-temperature limit it
reads

$\frac{1}{2}\langle{\left[\xi(t),\xi(t^{\prime})\right]_{+}\rangle}\simeq
2k_{b}TK(t-t^{\prime}).$ (57)

Using Eq. (53) and also the retarded solution, and considering
$C(\pm\infty)\simeq 0$, we get the in/out relation

$Y^{out}=Y^{in}-\int_{-\infty}^{+\infty}K(t-t^{\prime})\dot{X}(t^{\prime})dt^{\prime},$
(58)

where $Y^{in}=\sum_{k}c_{k}Y_{k}^{homg.}(-\infty)$,
$Y^{out}=\sum_{k}c_{k}Y_{k}^{homg.}(+\infty)$.

A useful way to treat the dissipation without having all the details of the
bath is to introduce the bath spectral density

$J(\omega)=\frac{\pi}{2}\sum_{k}\frac{c_{k}^{2}}{m_{k}\omega_{k}}\delta(\omega-\omega_{k}),$
(59)

and recast the dissipator in the form

$K(t)=\int_{0}^{\infty}\frac{d\omega}{\pi}\frac{J(\omega)}{\omega}\cos(\omega
t).$ (60)

Now all bath properties are encoded in the spectral density $J(\omega)$.

##  Appendix B Thermalizing master equation

We consider here the master equation suitable to study relaxation and
thermalization processes in cavity QED under the ultra-strong coupling regime.
For this purpose we consider the treatment used in Beaudoin  _et al._ (2011).
We do not repeat the derivation here, but we stress that the physical
assumptions are almost the same as used in deriving the Langevin equation in
Appendix A, with the further assumption that the coupling between the system
and the bath is very small. In particular this latter one ensures that we can
implement the rotating wave approximation between the system and the bath,
proceeding with the standard textbook derivation.

The crucial step here is to isolate the components of the system coupling
operator $X$ that rotates with positive (negative) frequencies. This can be
done as follows: given an Hamiltonian $H_{\rm sys}$ and a (or multiple) system
operator(s) $X$, we express them on the system eigenbasis
$X=\sum_{n,m}\braket{n}{X}{m}|n\rangle\langle m|$. The jump operators are then
given by the set $\\{{c_{nm}=\braket{n}{X}{m}|n\rangle\langle m|\,,~{}{\rm
such~{}that}~{}n<m\\}}$. In the Heinseberg picture these jump operators evolve
with positive frequencies. This allows to implement the rotating wave
approximation in the standard system-bath linear Hamiltonian in Eq. (50).

The master equation is then given by

$\partial_{t}\rho=\mathcal{L}_{H}(\rho)+\mathcal{L}_{D}(\rho),$ (61)

where the conservative time evolution is generated by

$\mathcal{L}_{H}(\rho)=-i\left[H_{\rm sys},\rho\right],$ (62)

while dissipations are given by

$\begin{split}&\mathcal{L}_{D}(\rho)=\sum_{n<m}\left[1+N_{T}(\omega_{mn})\right]\Gamma_{nm}D\left[|n\rangle\langle
m|,\rho\right]+\\\
&+\sum_{n<m}N_{T}(\omega_{mn})\Gamma_{nm}D\left[|m\rangle\langle
n|,\rho\right].\end{split}$ (63)

Here

$D\left[c,\rho\right]=c\,\rho\,c^{{\dagger}}-\frac{1}{2}\left[c^{{\dagger}}c\,,\,\rho\right]_{+}$
(64)

is the usual dissipator super-operator Gardiner and Zoller (2015), and

$N_{T}(\omega)=\frac{1}{e^{\omega/(k_{B}T)}-1}$ (65)

is the bosonic thermal population, where $k_{B}$ is the Boltzmann constant.
The thermalization rates are given by

$\Gamma_{nm}=J(|\omega_{mn}|)\left|\braket{n}{X}{m}\right|^{2}.$ (66)

Considering the thermal density matrix $\rho_{\rm ss}=e^{-H_{\rm
sys}/(k_{B}T)}/\mathcal{Z}$, where $\mathcal{Z}={\rm Tr}[\rho_{\rm ss}]$, one
can easily prove that is the steady state of the system.

##  Appendix C Diagonalization of the resonant-symmetric Rabi Hamiltonian

In this section we perform the approximated diagonalization of the Rabi model
in the regime where

$\begin{split}&\omega_{d}\simeq\omega_{c}\\\ &\epsilon=0.\end{split}$ (67)

We first transform the original Rabi Hamiltonian through the unitary
transformation $U_{\rm
pol}=\exp\left[g/\omega_{c}(a-a^{{\dagger}})s_{x}\right]$, obtaining the Rabi
polaron Hamiltonian in the form

$\begin{split}\tilde{H}_{\rm
Rabi}&=\omega_{c}a^{{\dagger}}a+\omega_{d}\Big{[}\cosh\left(g/\omega_{c}(a-a^{{\dagger}})\right)s_{z}\\\
&+i\sinh\left(g/\omega_{c}(a-a^{{\dagger}})\right)s_{y}\Big{]}.\end{split}$
(68)

Here $\cosh$, and $\sinh$ operators can be expressed in terms of the
displacement operator $\mathcal{D}(x)=\exp\left[x(a-a^{{\dagger}})\right]$.

In this frame one can then perform a generalized rotating wave approximation
Irish (2007), following from the fact that the Hamiltonian in this basis is
approximately block diagonal. Each block is spanned by the states
$\\{{|\uparrow,n-1\rangle,|\downarrow,n\rangle\\}}$, where $n=1,2\ldots$, and
the ground-state is given by the polaron vacuum state $|\downarrow,0\rangle$.

This block-diagonal structure is ultimately linked to the matrix elements of
the displacements operators in Eq. (68), which are known to be exponentially
suppressed as $\braket{n}{\mathcal{D}(g/\omega_{c})}{m}\sim
e^{-g^{2}/(2\omega_{c}^{2})}$ Cahill and Glauber (1969). Moreover because of
parity selection rule of the $\cosh(g/\omega_{c}(a-a^{{\dagger}}))s_{z}$,
$\sinh(g/\omega_{c}(a-a^{{\dagger}}))s_{y}$ operators only second nearest
neighbour blocks are coupled. Combining this two observations we have that the
most relevant transitions are within each block, for which we need only two
matrix element of the displacement operator per block:

$\begin{split}&\mathcal{D}_{n\,n-1}=\frac{g}{\omega_{c}}\sqrt{\frac{(n-1)!}{n!}}e^{-\frac{g^{2}}{2\omega_{c}^{2}}}L_{n-1}^{(1)}\left(g^{2}/\omega_{c}^{2}\right),\\\
&\mathcal{D}_{n\,n}=e^{-\frac{g^{2}}{2\omega_{c}^{2}}}L_{n}^{(0)}\left(g^{2}/\omega_{c}^{2}\right),\end{split}$
(69)

where $L_{n}^{(\alpha)}(x)$ are the special Laguerre polynomials. Notice that,
since $L_{n}^{(\alpha)}(0)=(n+\alpha)!/(n!)$, we have that
$L_{n-1}^{(1)}(0)=n$, recovering the usual Jaynes-Cummings picture at weak-
coupling.

We can then rewrite the polaron Rabi Hamiltonian as a block-diagonal matrix,
$\tilde{H}_{\rm Rabi}\approx\sum_{n}\tilde{H}_{\rm Rabi}^{n}$, where each
block reads

$\tilde{H}_{\rm Rabi}^{n}=\begin{pmatrix}A_{n}&C_{n}\\\ C_{n}&B_{n}\\\
\end{pmatrix}$ (70)

where

$\begin{split}&A_{n}=\omega_{c}n-\frac{\omega_{d}e^{-g^{2}/(2\omega_{c}^{2})}}{2}L_{n}^{(0)}\left(g^{2}/\omega_{c}^{2}\right),\\\
&B_{n}=\omega_{c}(n-1)+\frac{\omega_{d}e^{-g^{2}/(2\omega_{c}^{2})}}{2}L_{n-1}^{(0)}\left(g^{2}/\omega_{c}^{2}\right),\\\
&C_{n}=\frac{g}{\omega_{c}}\frac{\omega_{d}e^{-g^{2}/(2\omega_{c}^{2})}}{2}\sqrt{\frac{(n-1)!}{n!}}L_{n-1}^{(1)}\left(g^{2}/\omega_{c}^{2}\right).\end{split}$
(71)

The spectrum is

$\omega_{\pm,n}=\frac{A_{n}+B_{n}}{2}\pm\sqrt{\frac{\left(A_{n}+B_{n}\right)^{2}}{4}+C_{n}^{2}-A_{n}B_{n}}$
(72)

and the eigenstates are

$\begin{split}&|+,n\rangle=\cos\frac{\theta_{n}}{2}|\downarrow,n\rangle+\sin\frac{\theta_{n}}{2}|\uparrow,n-1\rangle,\\\
&|-,n\rangle=-\sin\frac{\theta_{n}}{2}|\downarrow,n\rangle+\cos\frac{\theta_{n}}{2}|\uparrow,n-1\rangle,\end{split}$
(73)

where

$\begin{split}&\cos\frac{\theta_{n}}{2}=\pm\sqrt{\frac{1}{2}\left(1+\frac{A_{n}-B_{n}}{\sqrt{\left(A_{n}-B_{n}\right)^{2}+4C_{n}^{2}}}\right)}\\\
&\sin\frac{\theta_{n}}{2}=\pm\sqrt{\frac{1}{2}\left(1-\frac{A_{n}-B_{n}}{\sqrt{\left(A_{n}-B_{n}\right)^{2}+4C_{n}^{2}}}\right)}\end{split}$
(74)

This approximate solution of the symmetric Rabi model is valid in all coupling
regimes for each value of $g$. However its validity is restricted to the cases
when $\omega_{d}\lesssim\omega_{c}$ Irish (2007).

##  Appendix D Matrix element and transition rates of the symmetric Rabi model

As for standard dressed states, the allowed transitions are only between
states of neighbouring blocks, with $(n,n\pm 1)$-excitations, and the only
relevant matrix elements contributing to the transition rates of the cavity
are

$\begin{split}\langle+,n|\left(a^{{\dagger}}-a\right)|-,n-1\rangle&=\sqrt{n-1}\cos\frac{\theta_{n-1}}{2}\sin\frac{\theta_{n}}{2}+\\\
&-\sqrt{n}\cos\frac{\theta_{n}}{2}\sin\frac{\theta_{n-1}}{2},\end{split}$ (75)
$\begin{split}\langle-,n|\left(a^{{\dagger}}-a\right)|+,n-1\rangle&=\sqrt{n-1}\sin\frac{\theta_{n-1}}{2}\cos\frac{\theta_{n}}{2}+\\\
&-\sqrt{n}\sin\frac{\theta_{n}}{2}\cos\frac{\theta_{n-1}}{2},\end{split}$ (76)
$\begin{split}\langle+,n|\left(a^{{\dagger}}-a\right)|+,n-1\rangle&=\sqrt{n-1}\sin\frac{\theta_{n-1}}{2}\sin\frac{\theta_{n}}{2}+\\\
&+\sqrt{n}\cos\frac{\theta_{n}}{2}\cos\frac{\theta_{n-1}}{2},\end{split}$ (77)
$\begin{split}\langle-,n|\left(a^{{\dagger}}-a\right)|-,n-1\rangle&=\sqrt{n-1}\cos\frac{\theta_{n-1}}{2}\cos\frac{\theta_{n}}{2}+\\\
&+\sqrt{n}\sin\frac{\theta_{n}}{2}\sin\frac{\theta_{n-1}}{2},\end{split}$ (78)

and for the ground-state

$\begin{split}&\langle\downarrow,0|\left(a+a^{{\dagger}}\right)|+,1\rangle=\cos\frac{\theta_{1}}{2},\\\
&\langle\downarrow,0|\left(a+a^{{\dagger}}\right)|-,1\rangle=-\sin\frac{\theta_{1}}{2}.\end{split}$
(79)

For the dipole we have a complementary situation

$\begin{split}&\langle+,n|s_{x}|-,n-1\rangle=-\frac{1}{2}\sin\frac{\theta_{n-1}}{2}\sin\frac{\theta_{n}}{2},\end{split}$
(80)

$\begin{split}&\langle-,n|s_{x}|+,n-1\rangle=\frac{1}{2}\cos\frac{\theta_{n-1}}{2}\cos\frac{\theta_{n}}{2},\end{split}$
(81)

$\begin{split}&\langle+,n|s_{x}|+,n-1\rangle=\frac{1}{2}\cos\frac{\theta_{n-1}}{2}\sin\frac{\theta_{n}}{2},\end{split}$
(82)

$\begin{split}&\langle-,n|s_{x}|-,n-1\rangle=-\frac{1}{2}\sin\frac{\theta_{n-1}}{2}\cos\frac{\theta_{n}}{2},\end{split}$
(83)

and the ground-state

$\begin{split}&\langle\downarrow,0|s_{x}|+,1\rangle=\sin\frac{\theta_{1}}{2},\\\
&\langle\downarrow,0|s_{x}|-,1\rangle=\cos\frac{\theta_{1}}{2}.\end{split}$
(84)

##  Appendix E Linear response and absorption spectra

Exciting the cavity corresponds to insert a current in the circuit, which can
be interpreted as a parallel LC filter. We can then define a system circuit
impedance $Z_{\rm sys}(\omega)$ for the LC circuit, which takes into account
the presence of the dipole in the capacitor. By using the standard composition
of circuit impedance, we can derive the response input/output relation from
the current flowing in the resistance (the transmission line)

$\frac{I_{out}}{I_{in}}=\frac{Z_{\rm sys}(\omega)}{R+Z_{\rm sys}(\omega)},$
(85)

where $R=Z_{\rm LC}Q$ is the Ohmic resistance of the transmission line coupled
to the LC cavity, and $Q=\omega_{c}/\gamma$ is the LC cavity quality factor.
In Fig. 1(c) it is shown the general scheme of our circuit approach.

The system impedance is defined by considering the relation between voltage
and current flowing through the circuit, $V=ZI$, which gives

$Z_{\rm sys}(\omega)=\frac{\braket{\dot{\Phi}}(\omega)}{I_{\rm in}(\omega)}.$
(86)

When the input current is very small we can invoke linear response theory
Clerk  _et al._ (2010), for which we have

$\chi_{VI}=\lim_{I_{\rm in}\rightarrow
0}\frac{\braket{\dot{\Phi}}(\omega)}{I_{\rm in}(\omega)},$ (87)

where $\chi_{VI}$ is the voltage-current linear response function Clerk  _et
al._ (2010). From here it follows an operative definition of the system
impedance as

$Z_{\rm sys}(\omega)=-i\omega\chi_{II}(\omega),$ (88)

where

$\chi_{II}=\lim_{I_{\rm in}\rightarrow 0}\frac{\braket{\Phi}(\omega)}{I_{\rm
in}(\omega)}$ (89)

is the current-current linear response function.

The current-current linear response function can be calculated in many ways,
but the simplest one is to use the cavity structure factor

$\mathcal{S}_{c}(\omega)=\sum_{n,m}\frac{e^{-\hbar\omega_{n}/(k_{b}T)}}{\mathcal{Z}}\left|\braket{n}{\Phi}{m}\right|^{2}\delta(\omega-\omega_{mn}).$
(90)

The system impedance is then given by $Z_{\rm
sys}(\omega)=-i\omega\mathcal{S}_{c}(\omega).$

## References

  * Landau and Lifshitz (2013) L. Landau and E. Lifshitz,  _Statistical Physics: Volume 5_, Bd. 5 (Elsevier Science, 2013).
  * Piskulich  _et al._ (2019) Z. A. Piskulich, O. O. Mesele, and W. H. Thompson,  _Activation Energies and Beyond,_  _The Journal of Physical Chemistry A_ , The Journal of Physical Chemistry A 123, 7185 (2019) .
  * Merzbacher (2002) E. Merzbacher,  _The Early History of Quantum Tunneling,_ Physics Today 55, 44 (2002) .
  * Caldeira and Leggett (1981) A. O. Caldeira and A. J. Leggett,  _Influence of Dissipation on Quantum Tunneling in Macroscopic Systems,_ Phys. Rev. Lett. 46, 211 (1981) .
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger,  _Dynamics of the dissipative two-state system,_ Rev. Mod. Phys. 59, 1 (1987) .
  * Milonni (1994) P. Milonni,  _The Quantum Vacuum: An Introduction to Quantum Electrodynamics_ (Elsevier Science, 1994).
  * Sáez-Blázquez  _et al._ (2022) R. Sáez-Blázquez, D. de Bernardis, J. Feist, and P. Rabl,  _Can we observe non-perturbative vacuum shifts in cavity QED?_ arXiv:2212.08675 [quant-ph] (2022).
  * Devoret  _et al._ (2007) M. Devoret, S. Girvin, and R. Schoelkopf,  _Circuit-QED: How strong can the coupling between a Josephson junction atom and a transmission line resonator be?*,_ Annalen der Physik 519, 767 (2007) .
  * De Bernardis  _et al._ (2018a) D. De Bernardis, T. Jaako, and P. Rabl,  _Cavity quantum electrodynamics in the nonperturbative regime,_ Phys. Rev. A 97, 043820 (2018a) .
  * Frisk Kockum  _et al._ (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori,  _Ultrastrong coupling between light and matter,_ Nature Reviews Physics 1, 19 (2019) .
  * Forn-Díaz  _et al._ (2019a) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano,  _Ultrastrong coupling regimes of light-matter interaction,_ Rev. Mod. Phys. 91, 025005 (2019a) .
  * Ciuti  _et al._ (2005) C. Ciuti, G. Bastard, and I. Carusotto,  _Quantum vacuum properties of the intersubband cavity polariton field,_ Phys. Rev. B 72, 115303 (2005) .
  * Hutchison  _et al._ (2012) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen,  _Modifying Chemical Landscapes by Coupling to Vacuum Fields,_ Angewandte Chemie International Edition 51, 1592 (2012) .
  * Paravicini-Bagliani  _et al._ (2019a) G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist,  _Magneto-transport controlled by Landau polariton states,_ Nature Physics 15, 186 (2019a) .
  * Valmorra  _et al._ (2021) F. Valmorra, K. Yoshida, L. C. Contamin, S. Messelot, S. Massabeau, M. R. Delbecq, M. C. Dartiailh, M. M. Desjardins, T. Cubaynes, Z. Leghtas, K. Hirakawa, J. Tignon, S. Dhillon, S. Balibar, J. Mangeney, A. Cottet, and T. Kontos,  _Vacuum-field-induced THz transport gap in a carbon nanotube quantum dot,_ Nature Communications 12, 5490 (2021) .
  * Appugliese  _et al._ (2022) F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist,  _Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect,_ Science 375, 1030 (2022) .
  * Flick  _et al._ (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio,  _Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry,_ Proceedings of the National Academy of Sciences 114, 3026 (2017) .
  * Fregoni  _et al._ (2022) J. Fregoni, F. J. Garcia-Vidal, and J. Feist,  _Theoretical Challenges in Polaritonic Chemistry,_  _ACS Photonics_ , ACS Photonics 9, 1096 (2022) .
  * Schäfer  _et al._ (2022) C. Schäfer, J. Flick, E. Ronca, P. Narang, and A. Rubio,  _Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity,_ Nature Communications 13, 7817 (2022) .
  * Paravicini-Bagliani  _et al._ (2019b) G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist,  _Magneto-transport controlled by Landau polariton states,_ Nature Physics 15, 186 (2019b) .
  * Arwas and Ciuti (2023) G. Arwas and C. Ciuti,  _Quantum electron transport controlled by cavity vacuum fields,_ Phys. Rev. B 107, 045425 (2023) .
  * Schuler  _et al._ (2020) M. Schuler, D. D. Bernardis, A. M. Läuchli, and P. Rabl,  _The vacua of dipolar cavity quantum electrodynamics,_ SciPost Phys. 9, 066 (2020) .
  * Ashida  _et al._ (2020) Y. Ashida, A. m. c. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler,  _Quantum Electrodynamic Control of Matter: Cavity-Enhanced Ferroelectric Phase Transition,_ Phys. Rev. X 10, 041027 (2020) .
  * Andolina  _et al._ (2020) G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini,  _Theory of photon condensation in a spatially varying electromagnetic field,_ Phys. Rev. B 102, 125137 (2020) .
  * Román-Roche  _et al._ (2021) J. Román-Roche, F. Luis, and D. Zueco,  _Photon Condensation and Enhanced Magnetism in Cavity QED,_ Phys. Rev. Lett. 127, 167201 (2021) .
  * Schlawin  _et al._ (2019) F. Schlawin, A. Cavalleri, and D. Jaksch,  _Cavity-Mediated Electron-Photon Superconductivity,_ Phys. Rev. Lett. 122, 133602 (2019) .
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef,  _Cavity quantum materials,_ Applied Physics Reviews 9, 011312 (2022) .
  * Bloch  _et al._ (2022) J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio,  _Strongly correlated electron–photon systems,_ Nature 606, 41 (2022) .
  * Garcia-Vidal  _et al._ (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen,  _Manipulating matter by strong coupling to vacuum fields,_ Science 373, eabd0336 (2021) .
  * Hausinger and Grifoni (2008) J. Hausinger and M. Grifoni,  _Dissipative dynamics of a biased qubit coupled to a harmonic oscillator: analytical results beyond the rotating wave approximation,_ New Journal of Physics 10, 115015 (2008) .
  * Mercurio  _et al._ (2022) A. Mercurio, V. Macrì, C. Gustin, S. Hughes, S. Savasta, and F. Nori,  _Regimes of cavity QED under incoherent excitation: From weak to deep strong coupling,_ Phys. Rev. Res. 4, 023048 (2022) .
  * Settineri  _et al._ (2018) A. Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta,  _Dissipation and thermal noise in hybrid quantum systems in the ultrastrong-coupling regime,_ Phys. Rev. A 98, 053834 (2018) .
  * Schaeverbeke  _et al._ (2019) Q. Schaeverbeke, R. Avriller, T. Frederiksen, and F. Pistolesi,  _Single-Photon Emission Mediated by Single-Electron Tunneling in Plasmonic Nanojunctions,_ Phys. Rev. Lett. 123, 246601 (2019) .
  * Magazzù  _et al._ (2021) L. Magazzù, P. Forn-Díaz, and M. Grifoni,  _Transmission spectra of the driven, dissipative Rabi model in the ultrastrong-coupling regime,_ Phys. Rev. A 104, 053711 (2021) .
  * Chen  _et al._ (2022) Z.-H. Chen, H.-X. Che, Z.-K. Chen, C. Wang, and J. Ren,  _Tuning nonequilibrium heat current and two-photon statistics via composite qubit-resonator interaction,_ Phys. Rev. Res. 4, 013152 (2022) .
  * Koch  _et al._ (2006) J. Koch, F. von Oppen, and A. V. Andreev,  _Theory of the Franck-Condon blockade regime,_ Phys. Rev. B 74, 205438 (2006) .
  * Leturcq  _et al._ (2009) R. Leturcq, C. Stampfer, K. Inderbitzin, L. Durrer, C. Hierold, E. Mariani, M. G. Schultz, F. von Oppen, and K. Ensslin,  _Franck–Condon blockade in suspended carbon nanotube quantum dots,_ Nature Physics 5, 327 (2009) .
  * Cui  _et al._ (2015) Y. Cui, S. Tosoni, W.-D. Schneider, G. Pacchioni, N. Nilius, and H.-J. Freund,  _Phonon-Mediated Electron Transport through CaO Thin Films,_ Phys. Rev. Lett. 114, 016804 (2015) .
  * Vdovin  _et al._ (2016) E. E. Vdovin, A. Mishchenko, M. T. Greenaway, M. J. Zhu, D. Ghazaryan, A. Misra, Y. Cao, S. V. Morozov, O. Makarovsky, T. M. Fromhold, A. Patanè, G. J. Slotman, M. I. Katsnelson, A. K. Geim, K. S. Novoselov, and L. Eaves,  _Phonon-Assisted Resonant Tunneling of Electrons in Graphene–Boron Nitride Transistors,_ Phys. Rev. Lett. 116, 186603 (2016) .
  * Jaako  _et al._ (2016) T. Jaako, Z.-L. Xiang, J. J. Garcia-Ripoll, and P. Rabl,  _Ultrastrong-coupling phenomena beyond the Dicke model,_ Phys. Rev. A 94, 033850 (2016) .
  * Forn-Díaz  _et al._ (2019b) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano,  _Ultrastrong coupling regimes of light-matter interaction,_ Rev. Mod. Phys. 91, 025005 (2019b) .
  * De Bernardis  _et al._ (2018b) D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl,  _Breakdown of gauge invariance in ultrastrong-coupling cavity QED,_ Phys. Rev. A 98, 053819 (2018b) .
  * Di Stefano  _et al._ (2019) O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori,  _Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics,_ Nature Physics 15, 803 (2019) .
  * Yoshihara  _et al._ (2017) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba,  _Superconducting qubit–oscillator circuit beyond the ultrastrong-coupling regime,_ Nature Physics 13, 44 (2017) .
  * Gardiner and Zoller (2015) C. Gardiner and P. Zoller,  _The Quantum World of Ultra-Cold Atoms and Light Book II: The Physics of Quantum-Optical Devices_ (IMPERIAL COLLEGE PRESS, 2015) https://www.worldscientific.com/doi/pdf/10.1142/p983.
  * Dalibard  _et al._ (1982) J. Dalibard, J. Dupont-Roc, and C. Cohen-Tannoudji,  _Vacuum fluctuations and radiation reaction : identification of their respective contributions,_ Journal de Physique 43, 1617 (1982) .
  * Kessler  _et al._ (2012) E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. Lukin, and J. I. Cirac,  _Dissipative phase transition in a central spin system,_ Phys. Rev. A 86, 012116 (2012) .
  * Minganti  _et al._ (2018) F. Minganti, A. Biella, N. Bartolo, and C. Ciuti,  _Spectral theory of Liouvillians for dissipative phase transitions,_ Phys. Rev. A 98, 042118 (2018) .
  * Johansson  _et al._ (2013) J. Johansson, P. Nation, and F. Nori,  _QuTiP 2: A Python framework for the dynamics of open quantum systems,_ Computer Physics Communications 184, 1234 (2013) .
  * Pilar  _et al._ (2020) P. Pilar, D. De Bernardis, and P. Rabl,  _Thermodynamics of ultrastrongly coupled light-matter systems,_ Quantum 4, 335 (2020) .
  * Rossatto  _et al._ (2017) D. Z. Rossatto, C. J. Villas-Bôas, M. Sanz, and E. Solano,  _Spectral classification of coupling regimes in the quantum Rabi model,_ Phys. Rev. A 96, 013849 (2017) .
  * Beaudoin  _et al._ (2011) F. Beaudoin, J. M. Gambetta, and A. Blais,  _Dissipation and ultrastrong coupling in circuit QED,_ Phys. Rev. A 84, 043832 (2011) .
  * Irish (2007) E. K. Irish,  _Generalized Rotating-Wave Approximation for Arbitrarily Large Coupling,_ Phys. Rev. Lett. 99, 173601 (2007) .
  * Claude Cohen-Tannoudji (1998) G. G. Claude Cohen-Tannoudji, Jacques Dupont-Roc,  _Atom-photon interactions: basic processes and applications_, wiley ed., Wiley Science Paperback Series (Wiley-VCH, 1998).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber,  _Ordered Expansions in Boson Amplitude Operators,_ Phys. Rev. 177, 1857 (1969) .
  * Garziano  _et al._ (2015) L. Garziano, R. Stassi, V. Macrì, A. F. Kockum, S. Savasta, and F. Nori,  _Multiphoton quantum Rabi oscillations in ultrastrong cavity QED,_ Phys. Rev. A 92, 063830 (2015) .
  * Ma and Law (2015) K. K. W. Ma and C. K. Law,  _Three-photon resonance and adiabatic passage in the large-detuning Rabi model,_ Phys. Rev. A 92, 023842 (2015) .
  * Niemczyk  _et al._ (2010) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross,  _Circuit quantum electrodynamics in the ultrastrong-coupling regime,_ Nature Physics 6, 772 (2010) .
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff,  _Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet,_ Phys. Rev. 58, 1098 (1940) .
  * Rabl (2011) P. Rabl,  _Photon Blockade Effect in Optomechanical Systems,_ Phys. Rev. Lett. 107, 063601 (2011) .
  * Minoguchi  _et al._ (2019) Y. Minoguchi, P. Kirton, and P. Rabl,  _Environment-Induced Rabi Oscillations in the Optomechanical Boson-Boson Model,_ arXiv:1904.02164 [quant-ph] (2019).
  * Rabl (2010) P. Rabl,  _Cooling of mechanical motion with a two-level system: The high-temperature regime,_ Phys. Rev. B 82, 165320 (2010) .
  * Breuer  _et al._ (2002) H. Breuer, F. Petruccione, and S. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Clerk  _et al._ (2010) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf,  _Introduction to quantum noise, measurement, and amplification,_ Rev. Mod. Phys. 82, 1155 (2010) .