# Nonperturbative Waveguide Quantum Electrodynamics

Yuto Ashida<EMAIL_ADDRESS>Department of Physics, University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Institute for Physics of
Intelligence, University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan  Ataç
$\dot{\mathrm{I}}$mamo$\breve{\mathrm{g}}$lu  Institute of Quantum
Electronics, ETH Zurich, CH-8093 Zürich, Switzerland  Eugene Demler
Department of Physics, Harvard University, Cambridge, MA 02138, USA

###### Abstract

Understanding physics of quantum emitters strongly interacting with quantized
electromagnetic continuum is one of the primary goals in the emergent field of
waveguide quantum electrodynamics (QED). When the light-matter coupling
strength is comparable to or even exceeds energies of elementary excitations,
conventional approaches based on perturbative treatment of light-matter
interactions, two-level description of matter excitations, and photon-number
truncation are no longer sufficient. Here we study in and out of equilibrium
properties of waveguide QED in such nonperturbative regimes by developing a
comprehensive and rigorous theoretical approach using an asymptotic decoupling
unitary transformation. In the transformed frame, factorization between light
and matter becomes exact in the strong-coupling limit and an accurate
effective model can be derived for arbitrary interaction strengths. We uncover
several surprising features ranging from symmetry-protected many-body bound
states in the continuum to strong renormalization of the effective mass and
potential; the latter may explain recent experiments demonstrating cavity-
induced changes in chemical reactivity as well as enhancements of
ferromagnetism or superconductivity. We demonstrate these results by applying
our general formalism to a model of coupled cavity arrays, which is relevant
to experiments in superconducting qubits interacting with microwave resonators
or atoms coupled to photonic crystals. We examine the relation between our
results and delocalization-localization transition in the spin-boson model;
notably, we point out that one can find a quantum phase transition akin to the
superradiant transition in multi-emitter waveguide QED systems with
superlinear photonic dispersion. Besides waveguide resonators, we discuss
possible applications of our framework to other light-matter systems relevant
to quantum optics, condensed matter physics, and quantum chemistry.

##  I Introduction

###  I.1 Background

Quantum many-body states arising from coherent interactions between light and
matter are of fundamental interest in modern physics and lie at the heart of
quantum technologies. Historically, analysis of light-matter systems focused
on the perturbative regime Rabi (1937); Cohen-Tannoudji  _et al._ (1989),
since interaction of atomic dipoles with vacuum electromagnetic fields is weak
due to the smallness of the fine structure constant $\alpha\\!=\\!1/137$.
Recent progress has led to the realization of new platforms in which
electromagnetic continuum is modified to reach stronger light-matter coupling.
In particular, (artificial) atoms coupled to one-dimensional continuum of
photons at microwave Bishop  _et al._ (2009); Astafiev  _et al._ (2010); Hoi
_et al._ (2011, 2012); van Loo  _et al._ (2013); Mlynek  _et al._ (2014);
Mirhosseini  _et al._ (2019); Kannan  _et al._ (2020) or optical Reitz  _et
al._ (2013); Thompson  _et al._ (2013); Arcari  _et al._ (2014); Yalla  _et
al._ (2014); Goban  _et al._ (2014); Lodahl  _et al._ (2015) frequencies
achieve enhancement of light-matter interaction through strong spatial
confinement of electromagnetic modes. This rapidly growing field of research
has been dubbed waveguide quantum electrodynamics (QED).

There exist many conceptual similarities between the questions addressed in
waveguide QED and the problems analyzed in the context of quantum dissipative
systems Leggett  _et al._ (1987); Schmid (1983); Guinea  _et al._ (1985);
Weiss (2012); in the latter, bosonic modes represent phonons or other
collective excitations of condensed matter systems. More recently, light-
matter interaction has also been the subject of intense research in the fields
of polaritonic chemistry Hutchison  _et al._ (2012); Galego  _et al._ (2015);
Ebbesen (2016); Herrera and Spano (2016); Flick  _et al._ (2017); Feist  _et
al._ (2017); Hiura  _et al._ (2018); Hiura and Shalabney (2019); Thomas  _et
al._ (2019); Haugland  _et al._ (2021) and nanostructured plasmonics Barnes
_et al._ (2003); Pitarke  _et al._ (2006); Tame  _et al._ (2013); Baumberg
_et al._ (2019); Mueller  _et al._ (2020). In light of such broad relevance,
models of quantum emitters interacting with a continuum of bosonic excitations
have played a crucial role not only in quantum information science, but also
in condensed matter physics and quantum chemistry. An outstanding challenge
here is to uncover the novel physical phenomena in nonperturbative regimes,
where strong interaction leads to the formation of quantum many-body states
with large entanglement between emitters and bosonic excitations of the
continuum.

Despite recent remarkable advances, our understanding of waveguide QED at
strong couplings is far from complete. Due to virtual excitation of many
photons, the problem becomes intrinsically nonperturbative and standard
approximations of quantum optics fail in many crucial aspects. First and
foremost, it is known that the usual rotating wave approximation becomes no
longer valid Cohen-Tannoudji  _et al._ (1989) due to the processes that create
or annihilate pairs of excitations. Moreover, the inclusion of the diamagnetic
$\hat{A}^{2}$ term and the multilevel structure of emitters becomes more
essential at larger coupling strengths. Importantly, the latter indicates that
the comprehensive understanding of strong coupling physics cannot be achieved
unless one goes beyond the standard two-level descriptions. Such multilevel
structure of quantum emitters is also of current technological importance. For
instance, superconducting transmon qubits are rarely operated as perfect two-
level systems Blais  _et al._ (2020) and such multilevel structures are
potentially useful for the purpose of performing certain quantum information
operations Rosenblum  _et al._ (2018); Elder  _et al._ (2020).

Figure 1:  (a) Schematic illustration summarizing the analysis performed in
the present manuscript. (Left) In the original Coulomb gauge, a single or
multiple emitters interact with common electromagnetic continuum in arbitrary
geometry via light-matter coupling $g$. A quantum emitter is modeled by a
charged quantum particle of mass $m$ and position $Q$ that is trapped in a
potential $V$. The potential is typically assumed to be a double-well
potential as appropriate for an effectively two-level emitter though our
theory is equally applicable to a generic potential profile. (Right) We use
the newly introduced unitary transformation to asymptotically decouple emitter
and photon degrees of freedom in the strong-coupling limit. After the
transformation, emitters and photons interact with each other via vanishingly
weak effective coupling that scales as $g_{\rm eff}\propto g^{-1/2}$ at large
$g$. In contrast, the renormalized mass is enhanced as $m_{\rm
eff}\\!\propto\\!g^{2}$, leading to the tight localization of the emitter at
the potential minima as well as the $1/g$ energy spacing. The potential is
renormalized to $V_{\rm eff}$ with lower potential barrier due to the dressing
by the vacuum electromagnetic fluctuations. (b) Formation of a ladder of the
bound states (BS) and the bound states in the continuum (BIC) on top of the
ground state (GS) in the nonperturbative regimes. The energy spacing and
excitation energies decrease as $\propto\\!g^{-1}$ and thus, these states
become increasingly degenerate at strong couplings.

While significant efforts have been devoted to elucidating the strong and
ultrastrong coupling regimes in this decade Zheng  _et al._ (2010); Koshino
and Nakamura (2012); Grießer and Ritsch (2013); Chang  _et al._ (2013);
González-Tudela and Porras (2013); Peropadre  _et al._ (2013); Ringel  _et
al._ (2014); Sanchez-Burillo  _et al._ (2014); Pichler  _et al._ (2015);
Calajó  _et al._ (2016); Shi  _et al._ (2016); Forn-Díaz  _et al._ (2017); Shi
_et al._ (2018); Gheeraert  _et al._ (2018); Martínez  _et al._ (2019);
Sánchez-Burillo  _et al._ (2019); Léger  _et al._ (2019); Román-Roche  _et
al._ (2020); Mahmoodian  _et al._ (2020); Forn-Díaz  _et al._ (2019); Sheremet
_et al._ (2021), the physics of waveguide QED in the realms of even stronger
light-matter interactions, namely, the deep Casanova  _et al._ (2010) and
extremely strong Ashida  _et al._ (2021) coupling regimes, remains largely
unexplored. There, the coupling strength becomes comparable to elementary
excitation energies or exceeds them, and qualitatively different phenomena are
expected to occur since vacuum fluctuations alone can lead to large
populations of photons in every coupled mode. However, due to the
aforementioned difficulties, a reliable theoretical approach for unveiling
those potentially new phenomena is currently lacking. One of the key goals of
this paper is to present a comprehensive theoretical framework that solves
this issue, allowing one to study strongly interacting light-matter systems in
the previously unexplored regimes.

On another front, the spin-boson model, a supposedly effective description of
waveguide QED systems (e.g., Ref. Peropadre  _et al._ (2013)), has long been
known to exhibit the delocalization-localization transition at strong
couplings Blume  _et al._ (1970); Caldeira and Leggett (1981). Nevertheless,
the breakdown of the usual two-level description Ashida  _et al._ (2021); De
Bernardis  _et al._ (2018) and the relevance of the diamagnetic term De
Bernardis  _et al._ (2018); Rzażewski  _et al._ (1975); Nataf and Ciuti
(2010); Viehmann  _et al._ (2011); Andolina  _et al._ (2019); Stokes and Nazir
(2020); De Liberato (2014); Garcia-Ripoll  _et al._ (2015); Ashida  _et al._
(2021) make it unclear how these known results for quantum dissipative systems
should be interpreted in the context of waveguide QED. More specifically, the
conditions under which a counterpart of the delocalization-localization
transition exists should be carefully examined by using the full-fledged QED
Hamiltonian. One intriguing possibility is that such a quantum phase
transition can be extended to multi-emitter systems and provide a new route
toward realizing a superradiant transition without external driving Ashida
_et al._ (2020a); Pilar  _et al._ (2020); Latini  _et al._ (2021).

In view of recent experimental developments in realizing stronger light-matter
interactions Forn-Díaz  _et al._ (2017); Yoshihara  _et al._ (2017); Mueller
_et al._ (2020), the time is ripe to explore in and out of equilibrium physics
of nonperturbative waveguide QED in a comprehensive manner. Specifically, we
will ask the following questions:

  * (A)

What are the defining physical features of waveguide QED in the previously
unexplored nonperturbative regimes?

  * (B)

How can one construct a proper effective model of waveguide QED at strong
couplings, where existing theoretical descriptions are expected to fail?

  * (C)

Starting from a fully microscopic theory, is it possible to identify a quantum
phase transition akin to the superradiant transition in waveguide QED setups?

The aim of this paper is to develop a consistent and rigorous framework to
address these questions from a unified perspective (see Fig. 1(a)). Below we
summarize the main results at a nontechnical level before presenting a
detailed theoretical formulation in subsequent sections.

###  I.2 Summary of the main results

The first part of this paper develops a general theory of quantum emitters
interacting with arbitrary quantum electromagnetic continuum. This is done by
extending the asymptotic light-matter decoupling unitary transformation that
we introduced earlier in the context of single-mode cavity QED Ashida  _et
al._ (2021) to the present waveguide QED setups. Previously, most of the
studies in waveguide QED approximate an emitter as a simplified two-level
system, which, however, is not necessarily appropriate for many experimentally
relevant systems, including superconducting qubits as mentioned above
Hutchings  _et al._ (2017); Elder  _et al._ (2020); Ma  _et al._ (2020);
Reinhold  _et al._ (2020). The validity of the two-level approximation becomes
particularly questionable in the nonperturbative regimes due to significant
renormalization of both the effective mass and potential as we demonstrate
later (see e.g., Fig. 5 in Sec. IV.3). To provide an adequate model of the
multilevel structure in realistic physical systems, in the present work we
model a quantum emitter as a charged particle moving in a potential with two
degenerate minima (see Fig. 1(a)).

Our first main result is the appearance of a ladder of many-body bound states
(BS) and the many-body bound states in the continuum (BIC) in nonperturbative
regimes (see Fig. 1(b)). We point out formation of increasingly many low-lying
bound states whose energies decrease as $\propto\\!g^{-1}$ in the limit of
strong light-matter coupling $g$. Interestingly, we find that the exact BICs
emerge as a consequence of the $\mathbb{Z}_{2}$-symmetry protection, which is
naturally satisfied in the standard QED Hamiltonian, without artificial fine-
tuning of microscopic parameters such as emitter positions/wavelength as
previously pursued Tufarelli  _et al._ (2013); Gonzalez-Ballestero  _et al._
(2013); Facchi  _et al._ (2016); Calajó  _et al._ (2019); Dinc  _et al._
(2019); Barkemeyer  _et al._ (2021). Even when the exact symmetry is broken,
the lifetime of these states diverges as $\propto\\!g^{3/2}$ in the strong-
coupling limit, and thus they still behave as so-called quasi BIC (see Fig. 3
in Sec. IV). We note that the BICs have recently attracted significant
attention in light of their potential applications to realize quantum memory
Lvovsky  _et al._ (2009) and nondissipative emitter interactions González-
Tudela  _et al._ (2015); Douglas  _et al._ (2015). It is in general
challenging to detect BICs in standard photon scattering experiments, since
bound states are orthogonal to delocalized states in the continuum. Instead,
we propose and numerically demonstrate an experimentally feasible quench
protocol to excite those states in a concrete model of cavity array, which is
fairly standard in quantum optics Hartmann  _et al._ (2006); Zhou  _et al._
(2008); Lombardo  _et al._ (2014). These results establish one of the defining
features in the nonperturbative regimes of waveguide QED and thus address
question (A).

We remark that the present work should be contrasted to earlier studies of
atom-field dressed bound states John and Wang (1990) in several crucial
aspects. In Refs. Longo  _et al._ (2010, 2011); Calajó  _et al._ (2016); Shi
_et al._ (2016); Kocabas (2016); Schneider  _et al._ (2016); Bello  _et al._
(2019); Mahmoodian  _et al._ (2020); Leonforte  _et al._ (2021); Kim  _et al._
(2021), the existence of bound states was predicted in perturbative regimes on
the basis of the rotating wave approximation. However, it turns out that these
bound states in general become resonances with finite lifetimes once the
counter rotating terms are included Sanchez-Burillo  _et al._ (2014); Román-
Roche  _et al._ (2020). In contrast, our analysis does not rely on those
simplifying approximations and rigorously establishes the presence of bound
states at arbitrary coupling strengths for general photonic dispersions. In
particular, a ladder of bound states revealed by our analysis is appreciable
only after multilevel structure of emitters is consistently included in
theory. Said differently, this feature cannot be addressed within any of
previous frameworks that rely on the two-level descriptions. More importantly,
to our knowledge, the existence of (quasi) BICs in nonperturbative waveguide
QED systems has been pointed out in this paper for the first time. We
emphasize that the BICs found here are genuine quantum many-body states that
are distinct from one-body wave phenomena, which have been the focus of
earlier studies Hsu  _et al._ (2016).

The second important result of our work is construction of proper effective
models for waveguide QED that remain valid at arbitrary coupling strengths.
This is made possible through the use of a unitary transformation that
achieves asymptotic decoupling of emitter and photon degrees of freedom in the
limit where light-matter interaction becomes the dominant energy scale (see
Fig. 1(a)). We point out that conventional descriptions become inapplicable in
the nonperturbative regimes because of uncontrolled level truncations in the
Coulomb or Power-Zienau-Woolley (i.e., dipole) gauges. In contrast, following
the unitary transformation used in the present work, such truncations are
well-justified owing to vanishingly small light-matter entanglement at strong
couplings, ensuring the validity of effective models constructed in this new
frame of reference. The obtained effective models take the same standard forms
as the Jaynes-Cummings-type Hamiltonian for a single-emitter case (see Eq.
(IV.3) in Sec. IV) and the inhomogeneous transverse-field Ising Hamiltonian
for a multi-emitter case (see Eq. (84) in Sec. V), but with suitably
renormalized parameters. More generally, our theory features a fully
microscopic treatment of waveguide QED and allows one to seamlessly analyze
the entire coupling region at both qualitative and quantitative levels. These
results thus address question (B).

Finally, building on these analyses, we answer question (C) in the affirmative
way. Specifically, we show that the infrared divergence of the renormalized
emitter mass occurs for a certain gapless photonic dispersion. This in turn
implies the exact two-fold degeneracy of the ground state and thus leads to
the transition to the symmetry-broken (i.e., localized) phase in the
thermodynamic limit. Our results are qualitatively consistent with the ones
previously obtained for the spin-boson model, except for the reentrance into
the delocalized phase in the extremely strong coupling regimes which
originates from the diamagnetic term (see Fig. 8 in Sec. VI). We then discuss
the extension of these results to multi-emitter cases and point out that a
waveguide QED system with superlinear photonic dispersion (corresponding to
sub-Ohmic bath) provides a promising candidate to realize the superradiant
transition (at least) for sufficiently small emitter frequencies.

Overall, it is notable that the key features revealed by this paper do not
rely on fine-tuning of parameters, but should appear generally in strongly
coupled light-matter systems. To obtain these results, it is crucial to
accurately perform analysis without resorting to uncontrolled approximations
that cannot be justified in the nonperturbative regimes. While the emphasis of
our discussion is on the waveguide setups, the present formalism is equally
applicable to other electromagnetic environments in arbitrary geometries.
Examples include cavity QED systems in 2D materials or polaritonic chemistry,
in which the inclusion of multiple photonic modes becomes crucial depending on
the cavity geometry and the coupling strength. Our work thus establishes a
foundation for studying strongly coupled light-matter systems lying at the
intersection of quantum optics, condensed matter physics, and quantum
chemistry in genuinely nonperturbative regimes.

The remainder of the paper is organized as follows. In Sec. II, we present a
general theoretical framework for a quantum emitter coupled to electromagnetic
continuum on the basis of the asymptotically decoupling unitary
transformation. In Sec. III, we unravel key physical features emerging in
nonperturbative regimes of waveguide QED. In Sec. IV, we illustrate the
general properties by providing an explicit numerical solution of a concrete
model of coupled cavity arrays. In Sec. V, we present the extension of the
theoretical formalism to multi-emitter systems. In Sec. VI, we consider the
ground-state properties of waveguide QED systems with a gapless photonic
dispersion and discuss their relation to the superradiant transition. In Sec.
VII, we give a summary of results and suggest several interesting directions
for future investigations.

##  II Asymptotic Light-Matter Decoupling: general formalism

We begin by developing a general theory of a single quantum emitter coupled to
arbitrary quantized electromagnetic environment. To analyze nonperturbative
regimes, we use a disentangling unitary transformation that can asymptotically
decouple light and matter degrees of freedom in the strong-coupling limit.
This asymptotically decoupled (AD) frame significantly simplifies the analysis
of strongly interacting light-matter systems, which allows us to explore the
entire coupling region, even beyond the ultrastrong coupling regimes. We will
later apply the framework to a concrete model of coupled cavity array in Sec.
IV. While a single-emitter setup is considered in this section, we will
generalize the whole formalism to multi-emitter cases in Sec. V.

###  II.1 QED Hamiltonian in the Coulomb gauge

We consider a quantum emitter that is locally coupled to quantized
electromagnetic modes in arbitrary geometries. The emitter is modeled as a
quantum particle of mass $m$ and charge $q$ trapped by a potential $V$, while
the electromagnetic environment is represented as a sum of harmonic
oscillators with frequencies $\omega_{k}$. The corresponding QED Hamiltonian
in the Coulomb gauge is given by

$\displaystyle\hat{H}_{{\rm
C}}=\frac{\bigl{(}\hat{P}-q\hat{A}\bigr{)}^{2}}{2m}+V(\hat{Q})+\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k},$
(1)

where $\hat{Q}$ ($\hat{P}$) is the position (momentum) operator of the emitter
and $\hat{a}_{k}$ ($\hat{a}_{k}^{\dagger}$) is the annihilation (creation)
operator of photons in mode $k$, which satisfy the commutation relations

$\displaystyle[\hat{Q},\hat{P}]=i\hbar,\;\;[\hat{a}_{k},\hat{a}_{k^{\prime}}^{\dagger}]=\delta_{kk^{\prime}}.$
(2)

We denote the vector potential operator as

$\displaystyle\hat{A}=\sum_{k}f_{k}(\hat{a}_{k}+\hat{a}_{k}^{\dagger}),$ (3)

where $f_{k}$ characterizes the electromagnetic amplitude of mode $k$.

It is useful to diagonalize the quadratic photon part of $\hat{H}_{\rm C}$ as
follows (see Appendix A for details):

$\displaystyle\hat{H}_{{\rm
C}}\\!=\\!\frac{\hat{P}^{2}}{2m}+V(\hat{Q})\\!-\\!\hat{P}\sum_{n}\zeta_{n}\left(\hat{b}_{n}\\!+\\!\hat{b}_{n}^{\dagger}\right)\\!+\\!\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n},$

where we perform the canonical transformation to introduce a squeezed photon
operator $\hat{b}_{n}$ labeled by $n\in\mathbb{Z}$ via

$\displaystyle\hat{a}_{k}=\sum_{n}\left(O\right)_{kn}\bigl{[}\cosh\left(r_{nk}\right)\hat{b}_{n}-\sinh\left(r_{nk}\right)\hat{b}_{n}^{\dagger}\bigr{]},$
(5)

and $\zeta_{n}$ is given as

$\displaystyle\zeta_{n}=\sqrt{\frac{\hbar}{m\Omega_{n}}}\sum_{k}g_{k}O_{kn}.$
(6)

Here, $O_{kn}$ is an orthogonal matrix that satisfies

$\displaystyle\sum_{kk^{\prime}}\left(O^{{\rm
T}}\right)_{nk}\left(\delta_{kk^{\prime}}\omega_{k}^{2}+2g_{k}g_{k^{\prime}}\right)O_{k^{\prime}m}=\delta_{nm}\Omega_{n}^{2},$
(7)

where $\Omega_{n}$ is an eigenfrequency of mode $n$, $r_{nk}$ is a squeezing
parameter defined by $e^{r_{nk}}\equiv\sqrt{\Omega_{n}/\omega_{k}}$, and
$g_{k}$ characterizes a coupling strength to mode $k$:

$\displaystyle g_{k}\equiv qf_{k}\sqrt{\frac{\omega_{k}}{m\hbar}}.$ (8)

We note that the magnitudes of $g_{k}$ depend on the system size $L$ (i.e.,
the total number $L$ of electromagnetic modes) via $g_{k}\propto f_{k}\propto
L^{-1/2}$.

Before proceeding further, we make two remarks. First, while we follow the
standard notation in atomic QED to write down the Hamiltonian (1), the present
formulation is equally applicable to circuit QED setups regardless of the
physical nature of each variable. In superconducting circuits, artificial
atoms are locally coupled to the continuum of microwave electromagnetic fields
in a transmission line. There is a well-established analogy between circuit
and atomic QED systems; the charge number operator of a transmon qubit and its
conjugate phase operator precisely correspond to $\hat{P}$ and $\hat{Q}$ in
Eq. (1), respectively, and the charge bias induced by the electromagnetic
fields of microwave resonator plays the role of the vector potential
$\hat{A}$. The same analogy holds true also for a flux qubit, where $\hat{Q}$
is coupled to photons through the dipole-type coupling $\hat{Q}\cdot\hat{E}$
with $\hat{E}$ being the electric field; one can use the Power-Zienau-Woolley
(PZW) transformation Power and Zienau (1959); Woolley (1971) to change this
circuit Hamiltonian back to the standard form as in Eq. (1) (see e.g., Ref. De
Bernardis  _et al._ (2018) or Eq. (III.4) below). In practice, coefficients of
the $\hat{A}^{2}$ term in circuit setups may have to be modified depending on
resonator geometries. Our formalism below can readily be generalized to
include such specifics.

Table 1: Summary of the scaling analysis for each of the renormalized parameters at different coupling strengths $g$ in Eq. (18) normalized by the characteristic photon frequency $\omega$ in Eq. (17). The second, third, and fourth columns represent the results in the ultrastrong coupling (USC), deep strong coupling (DSC), and extremely strong coupling (ESC) regimes, respectively. The renormalized frequencies $\Omega_{n}$ are determined by the eigenvalue problem (7). The label $n\\!=\\!0$ indicates the dominant electromagnetic mode with the largest eigenfrequency. The length scales $\xi_{n}$ in Eq. (11) are normalized by $x_{\omega}$ in Eq. (20) and characterize the effective light-matter coupling strengths for the dominant $n\\!=\\!0$ and the other modes $n\\!\neq\\!0$ in the transformed frame. The effective mass $m_{\rm eff}$ is defined by Eq. (15) with the renormalization factor (16). The two lowest rows correspond to the expectation values of the total photon numbers with respect to the low-energy eigenstates in the transformed frame (denoted by $U$) or in the Coulomb gauge (denoted by ${\rm C}$); see Eqs. (30) and (38). Parameter | USC | DSC | ESC
---|---|---|---
$g/\omega\sim 0.1$ | $g/\omega\sim 1$ | $g/\omega>1$
$\Omega_{0}/\omega$ |  $\simeq 1$ |  $\simeq\sqrt{1+2g^{2}/\omega^{2}}$ |  $\propto g$
$\Omega_{n\neq 0}/\omega$ |  $\simeq 1$ |  $\simeq 1$ |  $\propto g^{0}$
$\xi_{0}/x_{\omega}$ |  $\simeq g^{2}/\omega^{2}$ |  $\simeq\sqrt{\frac{g^{2}/\omega^{2}}{\left(1+2g^{2}/\omega^{2}\right)^{3/2}}}$ |  $\propto g^{-1/2}$
$\xi_{n\neq 0}/x_{\omega}$ |  $\simeq g/(\omega\sqrt{L})$ |  $={\rm O}(\delta^{2}/(g\omega\sqrt{L}))$ |  $\propto g^{-1}$
$m_{\rm eff}/m$ |  $\simeq 1$ |  $\simeq 1+2g^{2}/\omega^{2}$ |  $\propto g^{2}$
$\langle\sum_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}\rangle_{U}$ |  $={\rm O}(g^{2}/\omega^{2})$ |  $\simeq\frac{g^{2}/\omega^{2}}{\left(1+2g^{2}/\omega^{2}\right)^{5/2}}$ |  $\propto g^{-3}$, $g^{-2}$
$\langle\sum_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\rm C}$ |  $={\rm O}(g^{2}/\omega^{2})$ |  $={\rm o}(g^{2}/\omega^{2})$ |  $\propto g$

Second, we invoke neither the two-level approximation of an emitter nor the
rotating wave approximation (RWA), which are often used in the literature but
will break down when the light-matter interaction becomes sufficiently strong.
In particular, it will be crucial to take into account the multilevel
structure of an emitter to unveil the key physics in nonperturbative regimes
as we demonstrate later. We also note that the $\hat{A}^{2}$ term must be
incorporated to retain the gauge invariance of the theory, and its inclusion
becomes particularly essential when one goes beyond the ultrastrong coupling
regime. Meanwhile, we assume that the length scale of a quantum emitter is
much smaller than the photon wavelength in such a way that the $\hat{Q}$
dependence of the vector potential $\hat{A}$ can be neglected. This long-wave
assumption ultimately puts an upper limit on the light-matter coupling when
the confinement length scale of the emerging localized mode analyzed below
becomes comparable to the emitter size.

###  II.2 Asymptotic decoupling transformation

We now introduce a unitary transformation to asymptotically decouple light and
matter degrees of freedom Ashida  _et al._ (2021):

$\displaystyle\hat{U}=\exp\left(-\frac{i}{\hbar}\hat{P}\,\hat{\Xi}\right),$
(9)

where $\hat{\Xi}$ is given as

$\displaystyle\hat{\Xi}$ $\displaystyle\equiv$
$\displaystyle\sum_{n}i\xi_{n}(\hat{b}^{\dagger}_{n}-\hat{b}_{n}),$ (10)
$\displaystyle\xi_{n}$ $\displaystyle\equiv$
$\displaystyle\frac{\zeta_{n}}{\Omega_{n}}.$ (11)

This transformation acts on individual operators via

$\displaystyle\hat{U}^{\dagger}\hat{Q}\hat{U}$ $\displaystyle=$
$\displaystyle\hat{Q}+\hat{\Xi},$ (12)
$\displaystyle\hat{U}^{\dagger}\hat{b}_{n}\hat{U}$ $\displaystyle=$
$\displaystyle\hat{b}_{n}+\frac{\xi_{n}\hat{P}}{\hbar},$ (13)

where the emitter position is shifted by the gauge-field-dependent
displacement $\hat{\Xi}$ while each photon mode is subject to the momentum-
dependent shift $\xi_{n}\hat{P}/\hbar$. We note that the displacement
variables $\xi_{n}$ in Eq. (11) are chosen in such a way that the
$\hat{P}\\!\cdot\\!(\hat{b}\\!+\\!\hat{b}^{\dagger})$ term in Eq. (II.1) will
be precisely cancelled by the contributions arising from the displacement of
the $\hat{b}^{\dagger}\hat{b}$ term via Eq. (13).

The resulting Hamiltonian in the asymptotically decoupled frame is

$\displaystyle\hat{H}_{U}$ $\displaystyle=$
$\displaystyle\hat{U}^{\dagger}\hat{H}_{{\rm C}}\hat{U}$ (14) $\displaystyle=$
$\displaystyle\frac{\hat{P}^{2}}{2m_{{\rm
eff}}}+V(\hat{Q}+\hat{\Xi})+\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}.$

Here the effective mass is defined as

$\displaystyle m_{{\rm eff}}$ $\displaystyle\equiv$ $\displaystyle
m(1+2\Theta),$ (15) $\displaystyle\Theta$ $\displaystyle\equiv$
$\displaystyle\sum_{k}\left(\frac{g_{k}}{\omega_{k}}\right)^{2},$ (16)

where the mass enhancement is characterized by the dimensionless quantity
$\Theta$ whose expression (16) follows from Eq. (7). This renormalization
comes from the $\hat{P}^{2}$ terms arising from the residual contributions
generated by displacing the
$\hat{P}\\!\cdot\\!(\hat{b}\\!+\\!\hat{b}^{\dagger})$ and
$\hat{b}^{\dagger}\hat{b}$ terms. After the transformation, the light-matter
interaction is incorporated in the external potential $V$ in the form of the
gauge-field-dependent shift of the emitter, and its effective coupling
strength is characterized by $\xi_{n}$ instead of the bare coupling $g_{k}$.

Hereafter we first focus on the case of a gapped dispersion with frequencies
$\omega_{k}\\!>\\!0$ $\forall k$, for which $\Theta$ remains finite. This
includes experimentally relevant systems such as coupled cavity arrays and
open microwave transmission lines. The case of a gapless dispersion should be
analyzed separately, since one can find an infrared divergence of $\Theta$ in
that case; we will revisit this issue in Sec. VI.

###  II.3 Scaling analysis of the effective parameters

To demonstrate the asymptotic light-matter decoupling, we perform the scaling
analysis of the renormalized parameters with respect to the interaction
strength. To this end, we introduce the characteristic photonic frequency
$\omega$ and the coupling strength $g$ as follows:

$\displaystyle\omega^{2}$ $\displaystyle\equiv$
$\displaystyle\sum_{k}\omega_{k}^{2}/L,$ (17) $\displaystyle g^{2}$
$\displaystyle\equiv$ $\displaystyle\sum_{k}g_{k}^{2},$ (18)

where we note $g\\!=\\!{\rm O}(L^{0})$ since $g_{k}\\!\propto\\!L^{-1/2}$. We
begin by considering the regime $g/\omega\\!>\\!1$ in which the coupling
strength is dominant over other energy scales; we shall refer to it as the
extremely strong coupling (ESC) regime Ashida  _et al._ (2021). There, the
eigenvalue problem (7) has a single dominant mode (which we label $n\\!=\\!0$)
with the largest eigenfrequency $\Omega_{0}\\!\propto\\!g$ and
$\sum_{k}g_{k}O_{0k}\\!\simeq\\!g$, while the other frequencies remain
$\Omega_{n\neq 0}\\!=\\!{\rm O}(\omega)$ with $\sum_{k}g_{k}O_{n\neq
0k}\\!\propto\\!{\rm O}(g^{-1})$. This leads to the scalings $\xi_{0}\propto
g^{-1/2}$ and $\xi_{n\neq 0}\propto g^{-1}$, i.e., the light-matter
interaction in $\hat{H}_{U}$ asymptotically vanishes in the strong-coupling
limit.

In the deep strong coupling (DSC) regime $g/\omega\\!\sim\\!1$ Casanova  _et
al._ (2010), one can continue the scaling analysis in the similar manner and
obtain slightly refined expressions as summarized in Table 1. There, we denote
the variance of a photonic dispersion (or the effective bandwidth) as

$\displaystyle\delta^{2}=\sum_{k}(\omega_{k}-\omega)^{2}/L,$ (19)

and normalize the length scale by the characteristic one

$\displaystyle x_{\omega}=\sqrt{\frac{\hbar}{m\omega}}.$ (20)

In Table 1, we also summarize the scaling relations in the ultrastrong
coupling (USC) regime $g/\omega\sim 0.1$, which can readily be obtained by the
perturbative analysis. All these scalings will later be demonstrated in a case
study of coupled cavity array (see Fig. 2 below).

Besides the asymptotic decoupling in the strong-coupling limit, one notable
result of this scaling analysis is that the displacement parameters $\xi_{n}$
and consequently the effective light-matter couplings in the AD frame remain
small over the entire region of $g$. As detailed below, this fact allows us to
significantly simplify the analysis in a broad range of coupling strengths,
including the realms beyond the USC regime which are otherwise challenging to
investigate in any previous theoretical approaches.

###  II.4 Vacuum-dressed potential and decoupled excitations

To analyze low-energy eigenstates of $\hat{H}_{U}$, it is useful to rewrite it
in the following manner:

$\displaystyle\hat{H}_{U}=\hat{H}_{{\rm matter}}+\hat{H}_{{\rm
int}}+\hat{H}_{{\rm light}},$ (21)

where we define the matter Hamiltonian by

$\displaystyle\hat{H}_{{\rm matter}}=\frac{\hat{P}^{2}}{2m_{{\rm
eff}}}+V_{{\rm eff}}(\hat{Q})$ (22)

with $m_{\rm eff}$ being the renormalized mass (15) and $V_{\rm eff}$ being
the dressed potential given as

$\displaystyle V_{{\rm eff}}(Q)$ $\displaystyle\equiv$ $\displaystyle
V(Q)+\sum_{l=1}\frac{\xi^{2l}}{(2l)!!}V^{(2l)}(Q),$ (23)
$\displaystyle\xi^{2}$ $\displaystyle\equiv$
$\displaystyle\sum_{n}\xi_{n}^{2},$ (24)

where $V^{(l)}$ is the $l$-th derivative of $V$. The interaction Hamiltonian
is given by

$\displaystyle\hat{H}_{{\rm
int}}=\sum_{l=1}\frac{:\hat{\Xi}^{l}:}{l!}V^{(l)}(\hat{Q}),$ (25)

where $:\hat{O}:\equiv\hat{O}-\langle 0|\hat{O}|0\rangle$ represents the
normal ordering of photonic operators with $|0\rangle$ being the vacuum state
in the AD frame:

$\displaystyle\hat{b}_{n}|0\rangle=0\;\;\;\forall{n}.$ (26)

We emphasize that this vacuum state is distinct from the original vacuum of
$\hat{a}$ operators in the Coulomb gauge due to the squeezing (cf. Eq. (5)).
Finally, we denote the photon Hamiltonian as

$\displaystyle\hat{H}_{{\rm
light}}=\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}.$ (27)

Equations (23) and (25) can be obtained by expanding the interaction term in
Eq. (14) with respect to $\hat{\Xi}$ and using the relations $\langle
0|\hat{\Xi}^{2l}|0\rangle\\!=\\!(2l-1)!!\xi^{2l}$ and $\langle
0|\hat{\Xi}^{2l-1}|0\rangle\\!=\\!0$.

Arguably, the most celebrated signature of quantum fluctuations of the
electromagnetic field for an isolated ground-state atom is the Lamb shift.
Incorporation of light-matter coupling exclusively through a modification of
the external potential $V(\hat{Q}\\!+\\!\hat{\Xi})$ renders the origin of the
Lamb shift explicit; fluctuations in $\hat{\Xi}$ add to the intrinsic
fluctuations of $\hat{Q}$ to enhance the variance of the effective particle
position. This feature manifests itself as the nonvanishing dressing in the
effective potential (23) without any externally excited photons; it originates
from the zero-point fluctuations of the electromagnetic fields. Here, the
dominant contribution to the dressing strength $\xi$ comes from the mode $n=0$
and thus $\xi$ basically obeys the same scaling relation satisfied by
$\xi_{0}$ in Table 1. If necessary, the summation over $l$ can in practice be
truncated at a certain order that scales inversely with the coupling strength
$g$ owing to the asymptotic vanishment of $\xi$. In particular, in the limit
of weak light-matter coupling, the energy shift of the ground state is
captured by the $l\\!=\\!1$ term of Eq. (23). We emphasize however that the
effective-mass enhancement also modifies the ground-state energy.

Because of the decoupling $\xi_{0}\\!\propto\\!g^{-1/2}$ and enhancement of
the effective photon frequency $\Omega_{0}\\!\propto\\!g$, low-energy
eigenstates of $\hat{H}_{U}$ in the strong-coupling limit can be written as a
product of the emitter eigenstates and the photon vacuum as follows:

$\displaystyle|\Psi_{\alpha}\rangle_{U}=|\psi_{\alpha}\rangle|0\rangle,$ (28)

where $|\psi_{\alpha}\rangle$ with $\alpha=1,2,\ldots$ are single-particle
eigenstates of $\hat{H}_{\rm matter}$ in Eq. (22); we then represent it as

$\displaystyle\hat{H}_{\rm
matter}=\sum_{\alpha}E_{\alpha}(g)|\psi_{\alpha}\rangle\langle\psi_{\alpha}|$
(29)

with $E_{1}\leq E_{2}\leq\cdots$ being the corresponding eigenenergies. These
single-particle energies $E_{\alpha}$ provide asymptotically exact excitation
energies of the total Hamiltonian $\hat{H}_{U}$, which in the original Coulomb
gauge corresponds to an intrinsically many-body problem with highly entangled
light-matter degrees of freedom (cf. Eq. (1)). We note that the mass
enhancement $m_{\rm eff}\propto g^{2}$ leads to the tight localization of
$|\psi_{\alpha}\rangle$ around the bottom of $V_{\rm eff}$; accordingly, the
excitation energies $\delta E_{\alpha}\equiv E_{\alpha}-E_{1}$ decrease as
$\delta E_{\alpha}\propto g^{-1}$ in the nonperturbative regimes as long as
$V_{\rm eff}$ has well-defined minima.

###  II.5 Few-photon many-body ansatz at arbitrary coupling strengths

The scaling analysis indicates that the total number of photons in the AD
frame remains small over the entire coupling region (cf. Table 1). In
particular, in the ESC regime, the standard perturbation theory predicts that
the photon number in the ground state is on the order of
$(\xi_{0}/\Omega_{0})^{2}$ as long as the bandwidth $\delta$ is narrow,
resulting in the scaling

$\displaystyle\bigl{\langle}\sum_{n}\hat{b}^{\dagger}_{n}\hat{b}_{n}\bigr{\rangle}_{U}\propto
g^{-3},$ (30)

where $\langle\cdots\rangle_{U}$ represents an expectation value with respect
to a low-energy eigenstate in the AD frame. As the bandwidth becomes broad,
the contributions from $n\neq 0$ modes eventually dominate the $n\\!=\\!0$
contribution above, and the perturbation theory leads to $\sum_{n\neq
0}(\xi_{n}/\Omega_{n})^{2}\\!\propto\\!g^{-2}$, where we used $\xi_{n\neq
0}\\!\propto\\!\delta^{2}/(g\omega\sqrt{L})$ and $\Omega_{n\neq 0}\sim\omega$
(see Table 1). This crossover occurs when the bandwidth reaches a value around
$g\delta^{4}/\omega^{5}={\rm O}(1)$, at which the contributions from
$n\\!\neq\\!0$ and $n\\!=\\!0$ modes become comparable. In any cases, the
total photon number in the transformed frame asymptotically vanishes in the
strong-coupling limit.

This fact motivates us to introduce the few-photon many-body ansatz by
projecting the whole Hilbert space onto the following subspace:

$\displaystyle{\cal H}_{N_{c}}\equiv{\rm
span}\left\\{|\psi_{\alpha}\rangle\otimes|\psi_{{\rm
photon},i}\rangle\right\\},$ (31)

where we recall that $|\psi_{\alpha}\rangle$ are single-particle eigenstates
of $\hat{H}_{\rm matter}$ in Eq. (22) while we denote $|\psi_{{\rm
photon},i}\rangle$ as an arbitrary bosonic many-body state that satisfies

$\displaystyle\langle\psi_{{\rm
photon},i}|\sum_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}|\psi_{{\rm
photon},i}\rangle\leq N_{c}$ (32)

with a photon-number cutoff $N_{c}$; note that this cutoff is imposed on the
total photon number summed over all the modes, but not on each individual
electromagnetic mode. With this definition, the decoupled excitations (28)
discussed above correspond to the simplest subspace ${\cal H}_{0}$ with no
photons.

We here emphasize that the complexity is no longer exponential, but it is
polynomial with respect to the system size $L$; the Hilbert-space dimension of
the few-photon manifold ${\cal H}_{N_{c}}$ scales as $\propto L^{N_{c}}$. This
allows us to study the (exact) waveguide QED Hamiltonian (1) at arbitrary
coupling strengths in a very efficient and accurate manner. Indeed, our exact
diagonalization analysis shows that the results converge within (at most)
$\sim\\!1$% deviation already at a small total photon-number cutoff
$N_{c}\\!=\\!2$-$4$ (see Appendix B for further details).

This should be contrasted to previous approaches, where one has to include
more excitations for each electromagnetic mode at a greater coupling $g$. The
corresponding Hilbert-space dimension in general exponentially increases with
$L$, which severely limits their applicabilities in the strong-coupling
regions. Some variational states, such as the displaced-oscillator states
Silbey and Harris (1984); Bera  _et al._ (2014); Gheeraert  _et al._ (2018);
Román-Roche  _et al._ (2020), can provide useful approximative methods up to a
rather modest coupling regime. However, they should also ultimately become
inaccurate especially beyond the USC regime because of the breakdown of the
polaron picture. More importantly, the usual two-level description of an
emitter, on which most of the previous studies rely, becomes invalid once one
enters into the DSC and ESC regimes as detailed below. We will show that it is
actually such multilevel structure that leads to a defining feature of the
waveguide QED in genuinely nonperturbative regimes. Our approach gets around
these difficulties encountered in previous studies by employing the
(asymptotically exact) disentangling unitary transformation, after which the
whole low-energy eigenstates are well restricted into the few-photon manifold
(31) at any coupling strengths.

##  III Generic features of nonperturbative waveguide QED

We next delineate key physical features of waveguide QED that emerge when one
enters into the nonperturbative regimes on the basis of the theoretical
formalism developed in Sec. II. To understand the qualitative physics, it is
sufficient, as a first step, to consider the decoupled excitations (28) that
belong to the simplest, zero-photon subspace ${\cal H}_{0}$. The results
discussed here establish universal nonperturvative features which hold true
independent of specific choices of microscopic parameters. We will make these
predictions quantitatively accurate in the next section by extending the
analysis to the few-photon many-body ansatz in the subspace ${\cal
H}_{N_{c}>0}$.

###  III.1 Ladder of many-body bound states

One of the most surprising results is the appearance of a ladder of many-body
bound states. To see this, we recall that the excitation energies of the
decoupled states (28) decrease as $\delta E_{\alpha}\\!\propto\\!1/g$ due to
the mass enhancement, and eventually lie outside of the photon continuum,

$\displaystyle\delta E_{\alpha}\notin[\hbar\omega_{\rm L},\hbar\omega_{\rm
U}],$ (33)

where $\omega_{{\rm L}({\rm U})}$ represents the lower (upper) frequency limit
of the photon dispersion. These states are energetically separated from
scattering states and thus form bound states (BS), i.e., the excitation
energies are localized to the emitter degree of freedom and cannot decay to
the continuum at all. This emergence of multiple low-lying BS is inaccessible
by the commonly used two-level treatments that can be valid only up to the USC
regime. For this reason, the appearance of BS ladder can be considered as one
of the defining features of the DSC and ESC regimes of waveguide QED.

Interestingly, these bound states appear with equal energy spacing that
narrows as $\delta E\propto 1/g$. This results from the increase of the
emitter mass $m_{\rm eff}\\!\propto\\!g^{2}$ and the ensuing tight
localization of the wavefunction, which can be best understood in the AD
frame. The low-energy spectrum then reduces to the harmonic one as long as the
potential is well-behaved and can be expanded quadratically around the minima.
Importantly, in the original frame, these states behave as the many-body BS,
which are strongly entangled states including high-momentum emitter states and
exponentially localized (virtual) photons. It is worthwhile to note that
photon localization in these bound states becomes increasingly tight at
greater $g$ and can be much smaller than the (bare) emitter-transition
wavelength; this feature should be contrasted to usual atom-field dressed
bound states John and Wang (1990).

###  III.2 Many-body bound states in the continuum

Even when the decoupled excitations (28) lie within the photon continuum, they
can behave as either symmetry-protected bound states in the continuum (BIC) or
quasi BIC with lifetime that diverges in the strong-coupling limit. First, to
understand the origin of such symmetry-protected BIC, suppose that the
potential respects the inversion symmetry, $V(Q)\\!=\\!V(-Q)$. The total
light-matter Hamiltonian then satisfies the following $\mathbb{Z}_{2}$
symmetry:

$\displaystyle\hat{{\cal P}}^{-1}\hat{H}_{U}\hat{{\cal P}}$ $\displaystyle=$
$\displaystyle\hat{H}_{U},$ (34) $\displaystyle\hat{{\cal P}}^{2}$
$\displaystyle=$ $\displaystyle 1,$ (35)

where $\hat{\cal P}$ acts on the emitter operators as $\hat{{\cal
P}}^{-1}\hat{P}\hat{{\cal P}}=-\hat{P}$, $\hat{{\cal P}}^{-1}\hat{Q}\hat{{\cal
P}}=-\hat{Q},$ and also transforms the photon field via $\hat{{\cal
P}}^{-1}\hat{b}_{n}\hat{{\cal P}}=-\hat{b}_{n}$. It is clear that, if a
decoupled excitation lying in the photon continuum has a different parity from
that of scattering states, it leads to the exact BIC protected by the above
$\mathbb{Z}_{2}$ symmetry. For instance, the lowest photon continuum has the
odd parity, while a half of the decoupled excitations (28) have the even
parity and thus can behave as the BIC when the excitation energies lie within
the continuum.

Interestingly, even if a decoupled excitation has the same parity as
scattering states, it can still behave as a long-lifetime resonance, which is
often called quasi BIC. Indeed, the scaling analysis of its decay rate given
by the Fermi’s golden rule results in

$\displaystyle\Gamma_{\rm qBIC}\propto g^{-3/2},$ (36)

which vanishes in the strong-coupling limit. The same argument also applies to
the case when the $\mathbb{Z}_{2}$ symmetry is not exact due to, e.g., the
broken inversion symmetry, $V(Q)\\!\neq\\!V(-Q)$; the symmetry-protected BIC
discussed above then become resonances in this case, however, their lifetimes
still diverge in the strong-coupling limit. Physically, these (quasi) BIC
originate from the strong light-matter interaction containing the diamagnetic
effect, which tends to prevent scattering photons from interacting with the
many-body BS consisting of virtual photons localized around the emitter. We
emphasize that the physics of (quasi) BIC discussed here qualitatively remains
the same also in the case of a gapless dispersion unless $\Theta$ diverges
(see Sec. V).

It is worthwhile to note that fine-tuning of the coupling strength is not
necessary to observe the (quasi) BIC here. Specifically, there always exists a
nonzero-measure parameter regime of $g$ such that a certain excitation lies in
the photon continuum,

$\displaystyle\delta E_{\alpha}(g)\in[\hbar\omega_{\rm L},\hbar\omega_{\rm
U}].$ (37)

The reason for this is as follows. At zero coupling, one can find an emitter
state lying above the continuum, i.e., $\delta
E_{\alpha}(g\\!=\\!0)\\!>\\!\hbar\omega_{\rm U}$. In the ESC regime, this
excitation energy asymptotically decreases as $\delta
E_{\alpha}(g)\\!\propto\\!g^{-1}$ and ultimately converges to zero. Thus,
between these two limits, there must exist an intermediate coupling regime
such that the relation (37) is satisfied. The metastability of these states
stems from the fact that radiation field modes that are resonant with them
have small amplitudes at the emitter position.

###  III.3 Vacuum-induced suppression in potential barrier

Yet another common feature in the nonperturbative regimes is the vacuum-
induced suppression of potential barrier in $V_{\rm eff}$. This can readily be
understood from Eq. (23), where the vacuum fluctuations decrease (increase)
the energies at local maxima (minima), thus lowering the potential barrier
felt by the particle when tunneling to a different local minimum (see Fig. 4
below for an illustrative example of the double-well potential). The amount of
this suppression nonmonotonically depends on the coupling strength, since it
is solely determined by the displacement parameter $\xi$ that exhibits the
nonmonotonic $g$ dependence (see Eq. (24) and Table 1).

When one considers quantum tunneling, the mass renormalization eventually
dominates the barrier suppression and thus the tunneling rate is ultimately
exponentially suppressed in the strong coupling limit. In contrast, if thermal
activation over the barrier, i.e., thermal escape, is of interest, the escape
rate is basically characterized by the ratio of the potential barrier to the
temperature, but independent of the mass. Thus, it is a universal feature that
a thermal escape should be enhanced by strong light-matter couplings owing to
the vacuum-induced suppression of the barrier. This may provide a physical
explanation of recent experimental observations in polaritonic chemistry Hiura
_et al._ (2018); Hiura and Shalabney (2019); Thomas  _et al._ (2019), where
the thermally activated chemical reaction was found to be enhanced due to
cavity confinement.

###  III.4 Breakdown of level truncations in the Coulomb and PZW gauges

We finally point out that an analysis relying on the Coulomb or PZW gauges
must in general become invalid at a sufficiently strong light-matter coupling.
This difficulty arises from the breakdown of level truncations in photon and
emitter degrees of freedom in the nonperturbative regimes. Specifically, we
first note that the mean photon number in the Coulomb gauge grows as (cf.
Table 1)

$\displaystyle\bigl{\langle}\sum_{k}\hat{a}^{\dagger}_{k}\hat{a}_{k}\bigr{\rangle}_{{\rm
C}}\propto g.$ (38)

The same scaling also applies to the photon-number fluctuations. The number of
photon basis required to analyze the Coulomb-gauge Hamiltonian (1) thus
exponentially diverges as $g$ is increased. This eventually invalidates
truncation of photon levels, which is actually inevitable in almost any
analysis of bosonic many-body systems.

Similarly, the truncation of matter levels also becomes ill-justified at a
sufficiently large $g$ in the conventional gauges; in particular, this fact
indicates the breakdown of the usual two-level descriptions of quantum
emitters in the nonperturbative regimes. To see this, we use the unitary
transformation (9) to express the decoupled states (28) in the Coulomb gauge

$\displaystyle|\Psi_{\alpha}\rangle_{\rm C}$ $\displaystyle=$
$\displaystyle\hat{U}|\psi_{\alpha}\rangle|0\rangle$ (39) $\displaystyle=$
$\displaystyle\int dP\,\psi_{\alpha}(P)|P\rangle e^{-iP\hat{\Xi}}|0\rangle,$

where we expand an emitter state $|\psi_{\alpha}\rangle$ in terms of the
momentum eigenstates $|P\rangle$. In the strong-coupling limit, the variance
of the momentum distribution $|\psi_{\alpha}(P)|^{2}$ diverges with
$\sigma_{P}\\!\propto\\!g$ due to the mass renormalization $m_{\rm
eff}\\!\propto\\!g^{2}$. Thus, an energy eigenstate in the Coulomb gauge is a
strongly entangled emitter-photon state consisting of a superposition of
higher momentum states at larger $g$. This fact eventually invalidates the
common analyses that rely on either two-level approximation or a fixed
momentum cutoff for an emitter.

Note that these difficulties are carried over in the PZW gauge (also known as
the dipole gauge). To see this, we recall that the corresponding Hamiltonian
in the long-wavelength limit is given by $\hat{H}_{\rm
PZW}\\!=\\!\hat{U}^{\dagger}_{\rm PZW}\hat{H}_{\rm C}\hat{U}_{\rm PZW}$ with
the PZW transformation $\hat{U}_{\rm PZW}\\!=\\!\exp(iq\hat{Q}\hat{A}/\hbar)$:

$\displaystyle\hat{H}_{{\rm PZW}}$ $\displaystyle=$
$\displaystyle\frac{\hat{P}^{2}}{2m}+V(\hat{Q})+mg^{2}\hat{Q}^{2}+q\hat{Q}\hat{E}+\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k},$
$\displaystyle\hat{E}$ $\displaystyle\equiv$
$\displaystyle\sum_{k}if_{k}\omega_{k}(\hat{a}_{k}^{\dagger}-\hat{a}_{k}).$
(40)

As is the case with the Coulomb gauge, the mean/fluctuation of the photon
number in this gauge rapidly grows at strong couplings, while the mass remains
at the bare value which leads to eventual breakdown of matter-level
truncation. In contrast, the AD frame makes both photon- and emitter-level
truncations well-justified and allows us to reveal the key features in the
nonperturbative regimes as outlined above.

Figure 2:  (a) Renormalized photon frequencies $\Omega_{n}$ in Eq. (7), (b)
displacement parameters $\xi_{n}$ in Eq. (11), and (c) expectation values of
the total photon number plotted against the light-matter coupling strength $g$
in Eq. (45). In (a,b), the red solid (blue dashed) curve shows the values
corresponding to the dominant mode $n\\!=\\!0$ (the other modes $n\neq 0$),
where the dominant mode is characterized by the highest frequency
$\Omega_{0}\\!\simeq\\!\sqrt{\omega_{c}^{2}+2g^{2}}$. Note that only a part of
$n\\!\neq\\!0$ modes is plotted for the sake of visibility. In (c), the red
solid curves show the total number of dressed photons for the two-lowest
eigenstates in the asymptotically decoupled (AD) frame (cf. Eq. (30)), while
the blue dashed ones show the corresponding number of bare photons in the
Coulomb gauge (cf. Eq. (38)). Parameters are $J=0.2$ in (a,b), and $J=0.2$,
$v=0.5$, and $d=1.2$ in (c). All the numerical values are shown in the unit
$\omega_{c}\\!=\\!\hbar\\!=\\!m\\!=\\!1$ throughout this paper.

##  IV Application to coupled cavity array

Figure 3:  (Left) Low-energy excitation spectrum obtained by the exact
diagonalization (ED) of the AD-frame Hamiltonian (14) with the double-well
potential (46) at different coupling strengths $g$. The red dashed curves show
the energies corresponding to the decoupled excitations (28), which become
asymptotically exact in the strong-coupling limit. Panels (a) and (b) are
closeups of the main panel, where the exact BIC shows no anticrossings in (a),
while the quasi BIC exhibits tiny anticrossings in (b). Parameters are
$J=0.1$, $v=0.5$, $d=0.87$, and $L=19$. The emitter parameters correspond to
the on-resonant condition $\Delta/\omega_{c}\simeq 1$ with $\Delta$ being the
bare emitter frequency (see Eqs. (54) and (58) for the definition of
$\Delta$).

We here demonstrate all the generic features discussed in Sec. III by
analyzing a concrete model of coupled cavity arrays. Extending the above
analysis to the few-photon many-body ansatz (31), we provide experimentally
testable predictions of bound states, excitation energies, and quench
dynamics, which are quantitatively accurate over the entire coupling region.

We consider a waveguide realized by a one-dimensional array of coupled
cavities with nearest-neighbor coupling

$\displaystyle\hat{H}_{{\rm
light}}\\!=\\!-\frac{J}{2}\sum_{x}\left(\hat{a}_{x+1}^{\dagger}\hat{a}_{x}\\!+\\!{\rm
H.c.}\right)\\!+\\!\hbar\omega_{c}\\!\sum_{x}\hat{a}_{x}^{\dagger}\hat{a}_{x},$
(41)

where $J$ is a hopping parameter, $\omega_{c}$ is a resonator frequency, and
$\hat{a}_{x}\equiv\frac{1}{\sqrt{L}}\sum_{k}\hat{a}_{k}e^{-ikx}$ is a photonic
annihilation operator of the resonator mode at site $x$. The corresponding
dispersion is

$\displaystyle\hbar\omega_{k}=\hbar\omega_{c}-J\cos k$ (42)

with wavevector $k\in[-\pi,\pi)$. This specific choice of the waveguide is not
essential to the qualitative physics we discuss below, but is amenable to
numerical calculations and experimental implementations. The emitter is
coupled to the waveguide at $x\\!=\\!0$ and the vector potential in the
Coulomb-gauge Hamiltonian (1) is given by

$\displaystyle\hat{A}={\cal A}(\hat{a}_{x=0}+\hat{a}_{x=0}^{\dagger}),$ (43)

which corresponds to electromagnetic amplitudes $f_{k}$ (see Eq. (3))

$\displaystyle f_{k}=\frac{{\cal A}}{\sqrt{L}}.$ (44)

As the amplitudes are independent of $k$, it is useful to define the
characteristic light-matter coupling strength by

$\displaystyle g=q{\cal A}\sqrt{\frac{\omega_{c}}{m\hbar}}.$ (45)

We model a quantum emitter as a charged particle trapped in the standard
double-well potential,

$\displaystyle V(Q)=v\left(1-\frac{Q^{2}}{d^{2}}\right)^{2},$ (46)

where $v$ is a potential depth and $d$ characterizes a position of the
potential minima. While we here focus on this minimal model for a quantum
emitter, our theoretical formalism is equally applicable to a general
potential landscape that may be more complex depending on each specific system
or problem, such as transmon/flux qubits or chemical reactions.

Figure 2(a,b) shows the renormalized parameters in the AD frame at different
coupling strengths, which are fully consistent with the scaling analysis
summarized in Table 1. Specifically, beyond the USC regime, a single mode
labeled by $n\\!=\\!0$ turns out to have a large renormalized eigenfrequency
and dominantly couples to the emitter, while the other modes with
$n\\!\neq\\!0$ basically remain at the bare frequencies and are almost
decoupled from the emitter. As shown in Fig. 2(c), the total photon number in
the AD frame vanishes as $\propto\\!g^{-3}$ as consistent with the scalings
$\Omega_{0}\\!\propto\\!g$ and $\xi_{0}\\!\propto\\!g^{-1/2}$, while in the
Coulomb gauge the photon number increases as $\propto\\!g$.

###  IV.1 Bound states, symmetry-protected BIC, quasi BIC

Figure 3 shows the low-energy excitation spectrum of a quantum emitter coupled
to the cavity array in a broad range of the coupling strength $g$. This
spectrum is obtained by the exact diagonalization of the QED Hamiltonian in
the AD frame (14) within the few-photon many-body ansatz (31) (see Appendix B
for details about the method). We note that the emitter parameters $v$ and $d$
in Fig. 3 are chosen in such a way that the bare emitter frequency is resonant
to the cavity frequency, i.e., $(E_{2}-E_{1})/\hbar\simeq\omega_{c}$ at $g=0$.

The eigenstates insensitive to $g$ and staying in the photonic band,

$\displaystyle\delta E_{\rm sca}\in[\hbar\omega_{c}-J,\hbar\omega_{c}+J],$
(47)

correspond to the single-photon scattering states, which are extended over the
waveguide and constitute the energy continuum in the thermodynamic limit. In
contrast, the eigenstates lying out of the band continuum,

$\displaystyle\delta E_{\rm BS}\notin[\hbar\omega_{c}-J,\hbar\omega_{c}+J],$
(48)

behave as the bound states and are accompanied by virtual photons localized
around the emitter. In the nonperturbative regimes, the energies of these
bound states decrease as $\delta E_{\rm BS}\\!\propto\\!g^{-1}$ and are
asymptotically determined by the excitation energies $\delta E_{\alpha}$ of
the decoupled states (28). This point is confirmed in the left panel of Fig.
3, where the exact spectrum (blue solid curves) eventually agrees with the
asymptotic values (red dashed curves).

As discussed earlier, when these bound states lie in the band continuum,

$\displaystyle\delta E_{\rm(q)BIC}\in[\hbar\omega_{c}-J,\hbar\omega_{c}+J],$
(49)

they behave either as the $\mathbb{Z}_{2}$-symmetry-protected BIC or as the
quasi BIC. These features manifest themselves as the absence of anticrossings
in the finite-size spectrum (Fig. 3(a)) or as the presence of tiny
anticrossings with scattering states (Fig. 3(b)), respectively. We note that
this tiny anticrossing of the quasi BIC originates from its vanishingly small
decay rate, which can be estimated as (cf. Eq. (36) and the related
discussions in Sec. III.2)

$\displaystyle\Gamma_{\rm
qBIC}\sim\frac{(J/\hbar)^{2}}{g\sqrt{m\omega_{c}^{3}d^{3}}}\left(\frac{v^{3}}{m_{\rm
eff}}\right)^{1/4}\propto g^{-3/2}.$ (50)

In the ESC regime, the origin of these bound states can also be understood
from the fact that the cavity mode spatially overlapping with the emitter is
shifted in frequency outside the photonic continuum and thereby hopping to
neighboring cavities is strongly suppressed. All the excited light-emitter
states within the photonic continuum will then become (quasi) BIC because of
this suprression.

We also remark that, in Fig. 3, one can also find several continuum spectra
that connect the two-photon continuum with the single-photon one as $g$ is
increased. Physically, these states consist of the single-photon scattering
states on top of the first, second, and third excited bound states. The $g$
dependence of those energies can be understood as follows. They first rapidly
decrease with increasing $g$ until $g/\omega_{c}\\!\sim\\!5$. There, bound-
state energies are initially above the height of the potential barrier at
$Q\\!=\\!0$ and hence there are no double degeneracies. As we increase $g$
further and bound-state energies go well below the potential barrier, these
states become nearly degenerate because they now form a pair of symmetric and
antisymmetric combinations of excitations localized around each of the two
minima in the double-well potential. For instance, in Fig. 3 the energy of the
third excited state eventually approaches that of the second excited state and
they begin to overlap and become doubly degenerate from
$g/\omega_{c}\\!>\\!5$.

Figure 4:  Potential barrier $v_{\rm eff}$ of the vacuum-dressed effective
potential $V_{\rm eff}(Q)$ at different coupling strengths $g$ (see Eq. (52)).
The suppression and eventual restoration of the barrier arises from the
nonmonotonic $g$ dependence of $\xi$ in Eq. (24) (cf. Fig. 2(b)). Insets show
typical spatial profiles of the effective potential in each regime. Parameters
are $J=0.1$, $v=0.5$, and $d=1$ in (a), and $J=0.1$, $v=0.5$, and $d=0.6$ in
(b).

###  IV.2 Dressed potential

The effective emitter potential (23) in the AD frame is dressed by vacuum
electromagnetic fields. In the present case of the double-well potential, the
corresponding dressed potential is given by

$\displaystyle V_{{\rm eff}}(Q)=v_{{\rm eff}}\left(1-\frac{Q^{2}}{d_{{\rm
eff}}^{2}}\right)^{2},$ (51)

where we neglect an irrelevant constant, introduce the renormalized potential
barrier $v_{\rm eff}$ as

$\displaystyle v_{{\rm
eff}}=\begin{cases}v\left(1-\frac{3\xi^{2}}{d^{2}}\right)^{2}&\xi\leq\frac{d}{\sqrt{3}}\\\
0&\xi>\frac{d}{\sqrt{3}}\end{cases},$ (52)

and define the effective dipole length by $d_{\rm eff}/d\equiv(v_{\rm
eff}/v)^{1/4}$. As expected from the general argument in Sec. III.3, the
barrier $v_{\rm eff}$ is always suppressed compared to the bare value $v$ and
the suppression is most significant when $\xi$ becomes maximum, which occurs
around the DSC regime (see Fig. 4(a)). Interestingly, when the displacement
parameter $\xi$ becomes sufficiently large such that $\xi\\!>\\!d/\sqrt{3}$,
even the full suppression of the potential barrier, i.e., $v_{\rm eff}=0$, is
possible. Nevertheless, this does not mean that the entire potential is
suppressed because the dipole length $d_{\rm eff}$ in Eq. (51) also converges
to zero in this case; the resulting potential is the quartic one with the
single minimum (see Fig. 4(b) for an illustrative example).

###  IV.3 Two-level effective model and its breakdown in the Coulomb gauge

Construction of the Jaynes-Cummings-type effective model is often useful to
simplify the analysis of the original QED Hamiltonian, especially when low-
energy excitations are of interest. This can be done by performing the two-
level truncation of the emitter and assuming the rotating wave approximation.
In the AD frame, such procedure leads to the standard Jaynes-Cummings
Hamiltonian, but with the suitably renormalized parameters,

$\displaystyle\hat{H}_{U}^{\text{JC}}\\!=\\!\frac{\hbar\Delta_{g}}{2}\hat{\sigma}^{z}\\!+\\!\Bigl{(}\hat{\sigma}^{-}\sum_{n}\hbar\tilde{g}_{n}\hat{b}_{n}^{\dagger}\\!+\\!{\rm
H.c.}\Bigr{)}\\!+\sum_{n}\hbar\Omega_{n}\hat{b}^{\dagger}_{n}\hat{b}_{n},$

where the renormalized emitter frequency $\Delta_{g}$ and the effective
coupling strengths $\tilde{g}_{n}$ are defined by

$\displaystyle\Delta_{g}$ $\displaystyle\equiv$
$\displaystyle\frac{E_{2}-E_{1}}{\hbar}>0,$ (54) $\displaystyle\tilde{g}_{n}$
$\displaystyle\equiv$
$\displaystyle\frac{\xi_{n}}{\hbar}\langle\psi_{1}|\frac{dV}{dQ}|\psi_{2}\rangle.$
(55)

We recall that $E_{1,2}$ ($|\psi_{1,2}\rangle$) represent the two-lowest
eigenenergies (eigenstates) of the renormalized emitter Hamiltonian (22), and
thus depend on the coupling strength $g$ through $m_{\rm eff}$ and $V_{\rm
eff}$.

Figure 5:  Comparisons of the first-excited energy between the exact result
and the two-level effective models. The blue solid curve is obtained by the
exact diagonalizion (ED) of the AD-frame Hamiltonian (14). The green dots
represent the results of the two-level model with the rotating wave
approximation (RWA) in the AD frame (cf. Eq. (IV.3)). The red dashed curve
shows the corresponding results in the Coulomb gauge (cf. Eq. (IV.3)).
Parameters are $J=0.1$, $v=0.5$, and $d=0.87$.

Importantly, since the effective spin-bath couplings $\tilde{g}_{n}$ remain
small over the entire region (cf. Fig. 2(b)), the rotating wave approximation
in Eq. (IV.3) can be performed even beyond the USC regime. Also, the two-level
truncation for the emitter in the AD frame here should remain meaningful as
discussed in Sec. III.4. We thus expect the effective model
$\hat{H}_{U}^{\text{JC}}$ to be valid not only at weak $g$, but also in the
nonperturbative regimes.

Figure 6:  (a-c) Spatiotemporal dynamics of photon occupancy
$n_{x}=\langle\hat{a}^{\dagger}_{x}\hat{a}_{x}\rangle_{\rm C}$ after the
quench at different coupling strengths $g$. (a,b) Excitations of propagating
photons in the deep strong coupling regimes, and (c) photon confinement around
the emitter at $x\\!=\\!0$ in the extremely strong coupling regime. (d-f)
Corresponding dynamics of the photon occupancy $n_{x=0}$ at the emitter
position, and (g-i) the initial weights in terms of eigenstates of the post-
quench Hamiltonian (cf. Eq. (60)). In (d,e), the excitations of the (quasi)
BIC lead to the slow, long-lasting oscillatory dynamics whose period is
characterized by the bound-state energies. In (f), a ladder of bound states
manifests itself as the oscillatory dynamics with a long period $T_{\rm
osc}\\!=\\!2\pi/\omega_{\rm osc}$ that diverges in the strong-coupling limit
(cf. Eq. (62)). Parameters are $g=1.2$, $J=0.1$, $v=0.5$, $d_{i}=0.6$, and
$d_{f}=0.87$ in (a,d,g), $g=1.3$, $J=0.2$, $v=1$, $d_{i}=0.9$, and $d_{f}=2.7$
in (b,e,h), and $g=5$, $J=0.1$, $v=1$, $d_{i}=2$, and $d_{f}=2.5$ in (c,f,i).

To demonstrate this, we plot in Fig. 5 the lowest excitation energy $E_{\rm
ex}$ of $\hat{H}_{U}^{\text{JC}}$ which is obtained from the following
analytic relation for the single-excitation subspace,

$\displaystyle E_{{\rm
ex}}-\Delta_{g}=\sum_{n}\frac{\tilde{g}_{n}^{2}}{E_{{\rm ex}}-\Omega_{n}}.$
(56)

Indeed, it agrees well with the exact values even in the DSC regime. In
contrast, the conventional two-level model constructed from the Coulomb-gauge
Hamiltonian is given by

$\displaystyle\hat{H}^{\text{JC}}_{{\rm
C}}\\!=\\!\frac{\hbar\Delta}{2}\hat{\sigma}^{z}\\!+\\!\Bigl{(}\hat{\sigma}^{-}\sum_{k}\hbar\tilde{g}_{k}\hat{a}_{k}^{\dagger}\\!+\\!{\rm
H.c.}\Bigr{)}\\!+\\!\sum_{k}\hbar\omega_{k}\hat{a}^{\dagger}_{k}\hat{a}_{k},$

where $\Delta$ and $\tilde{g}_{k}$ are the bare parameters defined by

$\displaystyle\Delta$ $\displaystyle\equiv$ $\displaystyle\Delta_{g=0},$ (58)
$\displaystyle\tilde{g}_{k}$ $\displaystyle\equiv$
$\displaystyle\frac{g}{\sqrt{L}}x_{\omega_{c}}\langle\psi_{1}^{g=0}|\partial_{Q}|\psi_{2}^{g=0}\rangle$
(59)

with $x_{\omega_{c}}=\sqrt{\hbar/m\omega_{c}}$. This simplified Hamiltonian
$\hat{H}^{\text{JC}}_{{\rm C}}$ takes exactly the same form as
$\hat{H}^{\text{JC}}_{{U}}$ in Eq. (IV.3), but with unrenormalized parameters.
While this construction is valid at weak $g$, it is well-known that such a
naïve two-level description in the Coulomb gauge breaks down once one enters
into the USC regime in which nonresonant processes become relevant and the
two-level truncation becomes ill-justified (see the red dashed curve in Fig.
5).

In this respect, the AD frame significantly expands the applicability of the
Jaynes-Cummings description beyond the weak coupling regimes, and thus allows
one to use the standard techniques valid within the rotating wave
approximation, such as the Wigner-Weisskopf theory, in a broad range of $g$.
Nevertheless, we remark that the effective Hamiltonian
$\hat{H}^{\text{JC}}_{{U}}$ constructed in the AD frame should also ultimately
become invalid when $\Delta_{g}\\!<\\!\tilde{g}_{n}\\!\ll\\!\omega_{c}$, where
the counter rotating terms turn out to be equally important as rotating ones;
this typically occurs in $g\gtrsim 5$. Instead, a more accurate description
including the counter rotating terms, i.e., the quantum Rabi model, still
remains valid even in such ESC regime.

###  IV.4 Quench dynamics

The many-body bound states and the BIC lead to rich nonequilibrium dynamics in
the nonperturbative regimes. To be concrete, we consider the quench protocol
with the emitter parameter $d$ in the double-well potential (46) being
abruptly changed as $d_{i}\\!\to\\!d_{f}$ at time $t\\!=\\!0$ while keeping
all the other parameters constant. This effectively changes the positions of
the minima of $V_{\rm eff}$ and also modifies the qubit frequency. The initial
state $|\Psi(0)\rangle$ is chosen to be the ground state of the QED
Hamiltonian at $d\\!=\\!d_{i}$ and large fixed $g$. We emphasize that, in the
Coulomb gauge, this initial state is already a strongly entangled light-matter
state consisting of virtual photons localized around the emitter. The quench
protocol discussed here should be realized in the current experimental
techniques of, e.g., circuit QED that deals with photons in microwave regime.
As detailed below, this procedure provides a feasible way to experimentally
detect the signature of the predicted many-body BIC, which cannot be excited
by a single-photon scattering by definition.

We calculate the real-time dynamics by transforming to the AD frame, since the
analysis in the Coulomb gauge becomes exponentially hard at large $g$ as
discussed earlier. Specifically, we exactly diagonalize the post-quench
Hamiltonian $\hat{H}_{U}$ at $d\\!=\\!d_{f}$ (see Appendix B for details) and
obtain the time evolution via

$\displaystyle|\Psi(t)\rangle_{U}$ $\displaystyle=$ $\displaystyle
e^{-i\hat{H}_{U}t}|\Psi(0)\rangle_{U}$ (60) $\displaystyle=$
$\displaystyle\sum_{i}c_{i}e^{-iE_{i}t}|\Psi_{i}\rangle_{U},$

where $E_{i}$ ($|\Psi_{i}\rangle_{U}$) are the corresponding energies
(eigenstates), and $c_{i}$ are expansion coefficients of the initial state. We
then calculate the evolution of an observable $\hat{O}$ in the original
Coulomb gauge through the unitary transformation

$\displaystyle\langle\hat{O}\rangle_{\rm
C}\\!=\\!\langle\hat{U}^{\dagger}\hat{O}\hat{U}\rangle_{U}.$ (61)

Figure 6(a,b) shows the typical spatiotemporal dynamics of photon occupancy
$n_{x}\\!=\\!\langle\hat{a}_{x}^{\dagger}\hat{a}_{x}\rangle_{\rm C}$ in the
DSC regimes. One can find the nondecaying oscillatory dynamics that is most
pronounced around the emitter position $x\\!=\\!0$ as well as the emission of
propagating photons. The former originates from the existence of the many-body
bound states and the (quasi) BIC, while the latter can be considered as the
analogue of the dynamical Casimir effect in which physical photons are
generated by quenching the vacuum Yablonovitch (1989); Schwinger (1992).

To gain further insights into the oscillatory dynamics, we plot the time
evolution of the photon occupancy at the origin $x\\!=\\!0$ in Fig. 6(d,e). We
also show the corresponding initial weights $|c_{i}|$ in Fig. 6(g,h), where
the blue shaded regions represent the energy continuum. One can see that the
quench protocol excites the (quasi) BIC and the bound states with substantial
weights and that their frequencies characterize the long-period oscillation in
the dynamics, which typically has a period $T\\!=\\!{\rm O}(10)$. Thus, those
bound states manifest themselves as the long-lasting oscillatory behavior that
associates with photons bouncing back and forth around the emitter.

In the ESC regime, photons are so strongly bound by the emitter that they
cannot propagate into the waveguide (see Fig. 6(c)). Besides such photon
confinement, the dynamics exhibits the increasingly slow coherent oscillation
at larger $g$ (see Fig. 6(f)). This oscillatory behavior can be understood as
follows. In the AD frame, the emitter and photons are asymptotically
disentangled and the low-energy dynamics is solely governed by the
renormalized emitter Hamiltonian (22) with no photon excitations. The present
quench protocol then corresponds to the sudden shift of the potential minima
of $V_{\rm eff}$. Because the mass is enhanced as $m_{\rm
eff}\\!\propto\\!g^{2}$ and the wavepacket is tightly localized, this quench
initiates the oscillatory dynamics where the wavepacket (initially localized
at $d\\!\simeq\\!d_{i}$) oscillates around the new minima at
$d\\!\simeq\\!d_{f}$. Such oscillation frequency can be estimated as

$\displaystyle\omega_{{\rm osc}}=\sqrt{\frac{8v}{d_{f}^{2}m_{{\rm
eff}}}\left(1-\frac{3\xi^{2}}{d_{f}^{2}}\right)}\propto g^{-1},$ (62)

which vanishes in the strong-coupling limit. The estimated period $T_{\rm
osc}\\!=\\!2\pi/\omega_{\rm osc}$ agrees well with the observed long-period
oscillation in Fig. 6(f). In the energy basis, this slow coherent dynamics
manifests itself as excitations of a ladder of bound states corresponding to
the decoupled eigenstates (28) (see Fig. 6(i)).

##  V Extension to multiple quantum emitters

We now extend our theoretical formalism to the case of multiple quantum
emitters. We discuss several limiting cases and construct the effective two-
level model that is valid in a broad range of the light-matter coupling
strength.

###  V.1 General formalism

We consider $N$ emitters of mass $m_{j}$ that are subject to potential $V_{j}$
and locally interact with the common electromagnetic modes at positions
$x_{j}$ with $j\\!=\\!1,2,\ldots,N$. We thus start from the multi-emitter QED
Hamiltonian in the Coulomb gauge,

$\displaystyle\hat{H}_{{\rm
C}}\\!=\\!\\!\sum_{j=1}^{N}\Biggl{[}\frac{\left(\hat{P}_{j}\\!-\\!q\hat{A}_{x_{j}}\right)^{2}}{2m_{j}}\\!+\\!V_{j}(\hat{Q}_{j})\Biggr{]}\\!+\\!\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k},$
(63)

where the position and momentum operators of the emitters satisfy

$\displaystyle[\hat{Q}_{i},\hat{P}_{j}]=i\hbar\delta_{ij},$ (64)

and the vector potential is given by

$\displaystyle\hat{A}_{x_{j}}=\sum_{k}f_{kj}\left(\hat{a}_{k}e^{ikx_{j}}+\hat{a}_{k}^{\dagger}e^{-ikx_{j}}\right)$
(65)

with $f_{kj}$ characterizing an electromagnetic coupling between photonic mode
$k$ and emitter $j$; for the sake of simplicity, we assume $f_{kj}=f_{-kj}$.

Generalizing the unitary transformation (9) to such multi-emitter case, we
obtain the following Hamiltonian (see Appendix C for details):

$\displaystyle\hat{H}_{U}=\hat{H}_{{\rm matter}}+\hat{H}_{{\rm
int}}+\hat{H}_{{\rm light}},$ (66)

where the emitter part is given by

$\displaystyle\hat{H}_{{\rm
matter}}\\!=\\!\sum_{j}\Bigl{[}\frac{\hat{P}_{j}^{2}}{2m_{{\rm
eff},j}}\\!+\\!V_{{\rm
eff},j}(\hat{Q}_{j})\Bigr{]}\\!+\\!\\!\sum_{i>j}\mu_{ij}\hat{P}_{i}\hat{P}_{j}$
(67)

with

$\displaystyle m_{{\rm eff},j}$ $\displaystyle=$ $\displaystyle
m_{j}/\left[(1+2G)^{-1}\right]_{jj},$ (68) $\displaystyle\mu_{ij}$
$\displaystyle=$
$\displaystyle\left[(1+2G)^{-1}\right]_{ij}/\sqrt{m_{i}m_{j}}.$ (69)

Here, $G$ is the $N\\!\times\\!N$ matrix defined by

$\displaystyle G_{ij}$ $\displaystyle\equiv$
$\displaystyle\sum_{k}\frac{g_{ki}g_{kj}}{\omega_{k}^{2}}\cos(k(x_{i}-x_{j})),$
(70) $\displaystyle g_{kj}$ $\displaystyle\equiv$ $\displaystyle
qf_{kj}\sqrt{\frac{\omega_{k}}{m_{j}\hbar}}.$ (71)

Physically, $m_{{\rm eff},j}$ in Eq. (68) represents the renormalized mass
similar to $m_{\rm eff}$ in the single-emitter case considered before except
for its dependence on emitter positions through $G$. The coupling $\mu_{ij}$
in Eq. (69) represents the waveguide-mediated interaction between emitters; in
the original Coulomb gauge, this corresponds to the nondissipative coupling
mediated by virtual photons in the waveguide.

In the renormalized multi-emitter Hamiltonian (67), the vacuum-dressed
effective potentials $V_{{\rm eff},j}$ are given by

$\displaystyle V_{{\rm eff},j}(Q)$ $\displaystyle\equiv$ $\displaystyle
V_{j}(Q)+\sum_{l=1}\frac{\xi_{j}^{2l}}{(2l)!!}V_{j}^{(2l)}(Q),$ (72)
$\displaystyle\xi_{j}^{2}$ $\displaystyle\equiv$
$\displaystyle\sum_{n}|\xi_{nj}|^{2}.$ (73)

The interaction Hamiltonian is

$\displaystyle\hat{H}_{{\rm int}}=\sum_{j}:V_{j}(\hat{Q}_{j}+\hat{\Xi}_{j}):$
(74)

with

$\displaystyle\hat{\Xi}_{j}\equiv\sum_{n}i(\xi_{nj}\hat{b}_{n}^{\dagger}-\xi_{nj}^{*}\hat{b}_{n}),$
(75)

where the displacement parameters $\xi_{nj}$ characterize the effective
coupling strengths between dressed photon mode $n$ and emitter $j$. The photon
part $\hat{H}_{\rm light}$ takes the same form as the single-emitter case in
Eq. (27) 111We emphasize that dressed photon modes discussed here are in
general distinct from those in the single-emitter system, as inferred from the
sensitivity to emitter positions in multi-emitter cases (see Appendix C).
Nevertheless, we shall use the same subscript $n$ to label electromagnetic
mode for the sake of notational simplicity..

To be concrete, from now on we consider the case of identical emitters with
$m_{j}=m$, $V_{j}=V$, and $f_{kj}=f_{k}\;\forall j$. This simplification leads
to the identical bare light-matter couplings, $g_{kj}=g_{k}\;\forall j$. In
contrast, we note that the effective masses $m_{{\rm eff},j}$, the
displacement parameters $\xi_{nj}$, and the dressed potentials $V_{{\rm
eff},j}$ still depend on the emitter positions, and thus we need subscript $j$
to distinguish them. While we only illustrate key aspects of several limiting
cases below, the multi-emitter formulation developed here establishes a solid
starting point to accurately analyze arbitrary multi-emitter waveguide QED
systems in a broad range of the light-matter coupling strength.

###  V.2 Two emitters

Figure 7:  Renormalized parameters in the two-emitter case. (a) Waveguide-
mediated emitter-emitter interaction strength $\mu_{21}$ in Eq. (69) and (b)
effective mass $m_{\rm eff}$ in Eq. (68) are plotted against the emitter
separation at different light-matter couplings $g$. The waveguide is assumed
to be the coupled cavity array as in Sec. IV and we set $J=0.2$.

We begin by discussing the two-emitter case. Figure 7 plots the waveguide-
mediated coupling $\mu_{21}$ and the effective mass $m_{\rm eff}$ against the
emitter separation at different coupling strengths $g$. We here assume the
waveguide to be the same cavity array as considered in Sec. IV. In the USC
regime, the coupling $\mu_{21}$ exhibits the oscillatory behavior and can be
long-ranged. As $g$ is further increased, it becomes increasingly short-ranged
with oscillations being damped (see Fig. 7(a)). Such suppression can be
interpreted as a nonperturbative effect originating from the tighter
confinement of virtual photons around the emitters.

Figure 7(b) shows that the effective mass monotonically increases at larger
$g$. As noted earlier, the effective mass in multi-emitter systems is
sensitive to the emitter separation; for the two-emitter case, it starts from
$m\frac{1+4\Theta}{1+2\Theta}$ at zero separation and eventually saturates to
the single-emitter limit $m(1\\!+\\!2\Theta)$ when the separation surpasses
the cavity length $x_{\omega_{c}}\\!=\\!1$ (cf. Eq. (15) and (16)).

###  V.3 Localized $N$ emitters

We next consider the case in which all the emitters are coupled to the
waveguide at the same position $x_{1}\\!=\\!\cdots\\!=\\!x_{N}\\!=\\!0$. This
corresponds to the case of $N$ (artificial) atoms collectively coupled to the
common electromagnetic fields, which is relevant to various experimental
setups. It is useful to introduce the collective momentum and coordinate by

$\displaystyle\hat{P}_{{\rm CM}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{j}\hat{P}_{j},\;\;\;\;\hat{Q}_{{\rm
CM}}\equiv\frac{1}{\sqrt{N}}\sum_{j}\hat{Q}_{j},$ (76)

as well as the relative ones via

$\displaystyle\hat{p}_{j}$ $\displaystyle\equiv$
$\displaystyle\hat{P}_{j}-\frac{\hat{P}_{{\rm
CM}}}{\sqrt{N}},\;\;\;\;\hat{q}_{j}\equiv\hat{Q}_{j}-\frac{\hat{Q}_{{\rm
CM}}}{\sqrt{N}}.$ (77)

The Hamiltonian can be rewritten as

$\displaystyle\hat{H}_{U}=\hat{H}_{{\rm CM}}+\hat{H}_{{\rm int}}+\hat{H}_{{\rm
light}}+\hat{H}_{{\rm rel}},$ (78)

where

$\displaystyle\hat{H}_{{\rm CM}}=\frac{\hat{P}_{{\rm CM}}^{2}}{2M_{{\rm
eff}}}+N\,V_{{\rm eff}}\Bigl{(}\frac{\hat{Q}_{{\rm CM}}}{\sqrt{N}}\Bigr{)}$
(79)

governs the dynamics of the collective mode with the effective mass

$\displaystyle M_{{\rm eff}}=m\left(1+2N\Theta\right).$ (80)

The electromagnetic interaction between the collective mode and photons is
given by

$\displaystyle\hat{H}_{{\rm int}}=\;:N\,V\Bigl{(}\frac{\hat{Q}_{{\rm
CM}}}{\sqrt{N}}+\hat{\Xi}\Bigr{)}:.$ (81)

The relative motion of emitters is governed by

$\displaystyle\hat{H}_{{\rm
rel}}\\!=\\!\sum_{j}\Bigl{[}\frac{\hat{p}_{j}^{2}}{2m}\\!+\\!V\Bigl{(}\hat{q}_{j}\\!+\\!\frac{\hat{Q}_{{\rm
CM}}}{\sqrt{N}}\\!+\\!\hat{\Xi}\Bigr{)}\\!-\\!V\Bigl{(}\frac{\hat{Q}_{{\rm
CM}}}{\sqrt{N}}\\!+\\!\hat{\Xi}\Bigr{)}\Bigr{]}.$

Here we recall that the relative variables satisfy the constraints
$\sum_{j}\hat{p}_{j}\\!=\\!\sum_{j}\hat{q}_{j}\\!=\\!0$ and thus they contain
$N\\!-\\!1$ degrees of freedom.

Since the collective mode dominantly couples to the electromagnetic fields in
the present case, one may neglect dynamics of the relative degrees of freedom.
The total QED Hamiltonian (78) then becomes equivalent to the single-emitter
one (14) upon the replacements $\hat{Q}\\!\to\\!\hat{Q}_{{\rm CM}},$
$\hat{P}\\!\to\\!\hat{P}_{{\rm CM}}$, $g_{k}\\!\to\\!\sqrt{N}g_{k}$,
$d\\!\to\\!\sqrt{N}d$, and $v\\!\to\\!Nv$. This equivalence is nothing but the
well-known $\sqrt{N}$ collective enhancement of the dipole and the coupling
strength. Nevertheless, we emphasize that our analysis puts no limitations on
the coupling strength and rigorously shows that the $\sqrt{N}$ enhancement
(which was previously discussed up to the USC regime) is present at arbitrary
coupling strengths, even in the DSC/ESC regimes where the multilevel nature of
emitters becomes crucial.

###  V.4 Two-level effective model

We next extend the construction of the two-level effective model discussed in
Sec. IV.3 to arbitrary multi-emitter systems. The projection onto the two-
lowest dressed emitter states in the AD frame results in the effective model

$\displaystyle\hat{H}_{U}=\hat{{H}}_{\text{Ising}}+\hat{H}^{\text{JC}}_{\rm
int}+\hat{H}_{\rm light},$ (83)

where the matter part corresponds to the (inhomogeneous) transverse-field
Ising model,

$\displaystyle\hat{{H}}_{\text{Ising}}=\sum_{j}\frac{\hbar\Delta_{g,j}}{2}\,\hat{\sigma}_{j}^{z}+\sum_{i>j}J_{ij}\,\hat{\sigma}_{i}^{x}\hat{\sigma}_{j}^{x}$
(84)

with $J_{ij}$ being the waveguide-mediated qubit-qubit interaction given by

$\displaystyle
J_{ij}=\hbar^{2}\mu_{ij}\langle\psi_{1i}|{\partial}_{Q_{i}}|\psi_{2i}\rangle\langle\psi_{1j}|{\partial}_{Q_{j}}|\psi_{2j}\rangle.$
(85)

Here, $|\psi_{1,2j}\rangle$ represent the two-lowest eigenstates of the
renormalized emitter Hamiltonian $\hat{P}_{j}^{2}/2m_{{\rm
eff},j}\\!\\!+\\!\\!V_{{\rm eff},j}$, and $\Delta_{g,j}\\!>\\!0$ is the
corresponding excitation energy. We recall that $\Delta_{g,j}$ depends on
emitter positions through $m_{{\rm eff},j}$ and $V_{{\rm eff},j}$. The Jaynes-
Cummings-type light-matter interaction is

$\displaystyle\hat{H}^{\rm JC}_{\rm
int}=\sum_{j}\hat{\sigma}^{-}_{j}\sum_{n}\hbar\tilde{g}_{nj}\hat{b}_{n}^{\dagger}\\!+\\!{\rm
H.c.},$ (86)

where $\tilde{g}_{nj}$ are the effective qubit-boson couplings given by

$\displaystyle\tilde{g}_{nj}=\frac{\xi_{nj}}{\hbar}\langle\psi_{1j}|\frac{dV_{j}}{dQ}|\psi_{2j}\rangle.$
(87)

As discussed before, in contrast to the Coulomb/PZW gauges, our construction
in the AD frame remains valid even at large $g$ because the level truncations
can increasingly be well-justified in the strong-coupling limit. The two-level
effective model (83) thus should accurately capture low-energy physics of the
original multi-emitter QED Hamiltonian in a broad range of the coupling
strength.

In the single-emitter case, we recall that the asymptotic decoupling and the
enhanced photon frequency led to the decoupled eigenstates (28). Similarly, in
the present multi-emitter case, one can set the photon state to be the vacuum
and reduce the whole problem to the Ising Hamiltonian (84), from which the
ground-state properties and elementary excitations can be extracted. This
treatment becomes exact in the strong-coupling limit and should capture the
qualitative physics in the nonperturbative regimes at large, but finite $g$.
To make the results quantitatively accurate, one can extend the analysis to
the few-photon sector in the same manner as done in Sec. IV when necessary.

It is worthwhile to discuss a simple case of homogeneous configuration, in
which the emitters are periodically placed with the same separation. In this
case, the renormalized frequency and the spin-boson couplings are independent
of an emitter, $\Delta_{g,j}=\Delta_{g}$, $\tilde{g}_{nj}=\tilde{g}_{n}$
$\forall j$, and the spin-spin interaction becomes translationally invariant
$J_{ij}=J_{|i-j|}$. Then, the multi-emitter Hamiltonian (84) reduces to the
standard homogeneous transverse-field Ising model with couplings $J_{|i-j|}$
that are in general long-ranged. In particular, in the limit of zero emitter
separation, the effective Hamiltonian reduces to the Lipkin-Meshkov-Glick
(LMG) model:

$\displaystyle\hat{H}_{\rm
LMG}=\frac{\hbar\Delta_{g}}{2}\hat{S}^{z}+J^{\prime}(\hat{S}^{x})^{2},$ (88)

where $J^{\prime}\\!<\\!0$ is the all-to-all ferromagnetic coupling,
$\hat{S}^{\gamma}\\!\equiv\\!\sum_{j}\hat{\sigma}_{j}^{\gamma}$ are the
collective spin operators with $\gamma\\!\in\\!\\{x,y,z\\}$, and we neglect
the irrelevant constant.

These simplifications in the AD frame allow us to export the insights and
techniques originally developed in studies of the transverse-field Ising
models to the analysis of the multi-emitter QED Hamiltonian in the
nonperturbative regimes. Indeed, it is known that, in the many-emitter limit
$N\\!\to\\!\infty$, such model exhibits rich phase diagrams depending on
dimension, lattice geometry, or decay exponent of the long-range coupling
$J_{|i-j|}$ Suzuki  _et al._ (2012); Humeniuk (2016); Saadatmand  _et al._
(2018); Fey  _et al._ (2019). Moreover, a disordered transverse-field Ising
model is argued to realize many-body localization Kjäll  _et al._ (2014);
Hauke and Heyl (2015), and such disorder is fairly ubiquitous in the multi-
emitter Hamiltonian (84) where disorder comes into play through emitter
positions. While we leave the full understanding of such multi-emitter physics
at strong light-matter couplings to future investigations, our analysis
provides a reliable starting point for this and shows promise for realizing
the above exotic phases in the waveguide QED.

##  VI Quantum phase transitions with gapless dispersions

We finally turn our attention to the case of gapless dispersions.
Specifically, we consider the photon frequencies that scale as

$\displaystyle\omega_{k}\propto k^{l},$ (89)

which correspond to the spectral function

$\displaystyle
J(\omega)\equiv\pi\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k})\propto\omega^{1/l}$
(90)

with $l\\!>\\!0$ being an exponent characterizing the gapless dispersion. This
type of dispersions can be realized in waveguide QED systems by using
transmission lines or by designing mode frequencies with fabricated
resonators. Note that the couplings scale as $g_{k}\\!\propto\\!L^{-1/2}$ in
our notation (cf. Eq. (8)).

One of the key questions here is whether or not a waveguide QED system
governed by the Hamiltonian (1) undergoes a quantum phase transition as the
light-matter coupling is increased. For instance, in the single-emitter case
$N\\!=\\!1$, such a class of models constitutes a prototypical system for
investigating quantum dissipation Leggett  _et al._ (1987); Schmid (1983);
Guinea  _et al._ (1985); Weiss (2012), and it has been known that (depending
on $l$) a quantum phase transition can occur in the thermodynamic limit of
bosonic modes $L\\!\to\\!\infty$.

Below we revisit this fundamental problem and discuss the related issues from
a new perspective of the AD frame. We show that the presence or absence of
transition can be understood in a very simple and unified manner in terms of
the mass renormalization after the unitary transformation. We then extend the
analysis to multi-emitter cases and point out that the superradiant phase,
which is known to be absent in cavity QED setups with a single photonic mode,
can in general occur in waveguide QED systems with a superlinear gapless
dispersion.

###  VI.1 Delocalization-localization transition in the single-emitter case

The ground state of a single-emitter system displays either delocalized or
localized phase that is characterized by the following order parameter

$\displaystyle q$ $\displaystyle\equiv$
$\displaystyle\lim_{h\to+0}\lim_{L\to\infty}\langle\hat{Q}\rangle_{h,{\rm C}}$
(91) $\displaystyle=$
$\displaystyle\lim_{h\to+0}\lim_{L\to\infty}\langle\hat{Q}+\hat{\Xi}\rangle_{h,U},$
(92)

where $\langle\hat{Q}\rangle_{h,{{\rm C}}}$ represents an emitter displacement
in the ground state of the QED Hamiltonian in the Coulomb gauge (1) with a
bias potential $-h\hat{Q}$ being added to $V(\hat{Q})$; we here assume $V$ to
be the standard double-well potential (46) though our arguments can be applied
to a generic potential profile. In the AD frame, this order parameter
corresponds to an expectation value of
$\langle\hat{Q}\\!+\\!\hat{\Xi}\rangle_{U}$ (cf. Eq. (61)). The delocalized
phase is characterized by the vanishing order parameter $q\\!=\\!0$ and the
unique, nondegenerate ground state. In contrast, in the localized phase, the
ground state exhibits the two-fold degeneracy in the thermodynamic limit
$L\\!\to\\!\infty$ corresponding to localization to each of the two minima in
the potential. In this case, the order parameter takes a nonzero value
$q\\!>\\!0$, which indicates the broken $\mathbb{Z}_{2}$ symmetry in Eq. (34).

To understand whether or not the ground state exhibits the two-fold degeneracy
in the thermodynamic limit, the renormalization factor $\Theta$ in the
effective mass $m_{\rm eff}$ plays a crucial role (cf. Eq. (16)). Depending on
the exponent $l$, it scales with the system size $L$ as

$\displaystyle\Theta\propto\begin{cases}{\rm
O}(1)&l<1\;\;({\text{sublinear/super-Ohmic}})\\\
L^{l-1}&l>1\;\;(\text{superlinear/sub-Ohmic})\\\
\ln(L)&l=1\;\;(\text{linear/Ohmic})\end{cases},$ (93)

where the divergence in the superlinear and linear regimes originates from the
infrared one. It is worthwhile to note that recent experiments have indeed
reached such regime that the infrared divergence becomes relevant, the so-
called superstrong coupling regime Kuzmin  _et al._ (2019).

At sufficiently large coupling $g$, the emitter and photons are decoupled in
the AD frame, and the first excitation energy $\Delta_{g}$ can be estimated
from the tunneling rate in the dressed potential (23), resulting in Jona-
Lasinio  _et al._ (1981)

$\displaystyle\hbar\Delta_{g}\simeq\frac{\hbar^{2}}{m_{\rm eff}d_{\rm
eff}^{2}}\exp\left[-\frac{4}{3}\sqrt{\frac{2m_{\rm eff}d_{\rm eff}^{2}v_{\rm
eff}}{\hbar^{2}}}\right].$ (94)

Thus, the divergent $\Theta$ leads to the gap closing
$\lim_{L\to\infty}\Delta_{g}\\!=\\!0$, indicating a possible ground-state
degeneracy, i.e., transition to the localized phase. In contrast, as long as
$m_{\rm eff}$ remains finite, there should be no such exact two-fold
degeneracy of the ground state even in the thermodynamic limit; this implies
the absence of transition. Below we delineate general properties in each case
of dispersions on the basis of this observation.

#### Sublinear

When the dispersion is sublinear with $l\\!<\\!1$ (also known as the super-
Ohmic bath), the integral in $\Theta$ converges and thus the effective mass
always remains finite $m_{\rm eff}\\!<\\!\infty$. Thus, the whole results
discussed in Secs. II-IV for a gapped dispersion qualitatively remain the
same, except for the point that all the bound states now behave as the (quasi)
BIC in the present gapless case. Importantly, the ground state can thus be
well-approximated by the lowest decoupled eigenstate
$|\psi_{0}\rangle|0\rangle$ (cf. Eq. (28)), which has $q\\!=\\!0$ and provides
the unique ground state due to the nonvanishing excitation energy
$\Delta_{g}\\!>\\!0$ (see Eq. (94)). Note that, while $\Delta_{g}$ can be
exponentially small as $g$ is increased, it still remains nonzero in
$L\\!\to\\!\infty$ at any finite $g$. Thus, the system cannot exhibit the
exact two-fold degeneracy, and we conclude that it resides in the delocalized
phase at arbitrary coupling strengths in a sublinear case $l\\!<\\!1$. It is
worthwhile to note that the same argument also rules out the possibility of a
phase transition for general gapped dispersions, which include some
experimentally relevant situations, such as cavity array and open transmission
lines.

#### Superlinear

Figure 8:  Schematic figure illustrating the ground-state phase diagram of the
full QED Hamiltonian (1) with a superlinear photonic dispersion and the
double-well emitter potential (46). It is qualitatively consistent with what
has been found in the sub-Ohmic spin-boson model, except that the QED system
here should ultimately reenter into the delocalized phase at a sufficiently
large light-matter coupling due to the diamagnetic term (red dashed horizontal
arrow).

The integral in $\Theta$ exhibits the infrared divergence for superlinear
dispersion with $l\\!>\\!1$ (often called the sub-Ohmic case), where the
effective mass grows as $m_{\rm eff}\\!\propto\\!L^{l-1}$. One can also check
that such superlinear dispersion leads to the polynomial divergence of $\xi$
in the dressed emitter potential (23); the resulting effective potential
$V_{\rm eff}$ then has the unique minimum at $Q\\!=\\!0$ (cf. Eq. (52)) and
becomes infinitely tight in $L\\!\to\\!\infty$ 222While the infrared
divergence also occurs for linear dispersion, it associates with (slow)
logarithmic size dependence and one should discuss this case separately from
the present superlinear case.. Thus, in the thermodynamic limit, the emitter
wavefunction in the AD frame is completely localized at $Q\\!=\\!0$, and the
total system is solely governed by the photonic part

$\displaystyle\hat{H}_{U}=\;:V(\hat{\Xi}):+\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}.$
(95)

An order of its ground state is characterized by the expectation value of
$\langle\hat{\Xi}\rangle_{U}$ (see Eq. (92)). Building up on this observation,
we now argue that both delocalized and localized phases are possible depending
on the potential depth $v$.

First, in the limit of a deep potential with large $v$, the first term in Eq.
(95) is dominant, and the ground state should exhibit the two-fold degeneracy
at $\langle\hat{\Xi}\rangle_{U}\\!=\\!\pm d$ in $L\\!\to\\!\infty$. This leads
to the localized phase with $q>0$. Second, in the opposite limit of a shallow
potential, the second (harmonic) term in Eq. (95) dominates over the potential
term, and the vacuum state gives the unique ground state, which has
$\langle\hat{\Xi}\rangle_{U}\\!=\\!0$ and leads to the delocalized phase with
$q\\!=\\!0$. Finally, in between these two limits, the potential landscape
effectively changes from the double-well profile to the harmonic one, and
accordingly, the order parameter $q$ continuously decreases from $d$ and
vanishes at certain $v^{*}\\!>\\!0$. We recall that, in general, a deep
(shallow) potential depth $v$ corresponds to a small (large) bare qubit
frequency $\Delta$. Altogether, one expects a continuous quantum phase
transition between the localized phase at strong $v$ (resp. small $\Delta$)
and the delocalized phase at weak $v$ (resp. large $\Delta$); see the blue
dashed vertical arrow in Fig. 8. We remark that this conclusion is fully
consistent with the numerical renormalization group results of the sub-Ohmic
spin-boson model (see e.g., Ref. Bulla  _et al._ (2003)).

The present consideration of the full QED Hamiltonian also reveals an
intriguing new possibility beyond the usual spin-boson model. Namely, due to
the nonmonotonic $g$ dependence of $\xi_{n}$ (cf. Fig. 2(b)), the system
should again transition into the delocalized phase at a sufficiently large
light-matter coupling (see the red dashed horizontal arrow in Fig. 8). This is
expected to occur in the ESC regime, where the coupling strength dominates the
other energy scales. Physically, the origin of this favoring of the
delocalized phase can be traced back to the diamagnetic $\hat{A}^{2}$ term
that suppresses the displacements of the bosonic modes from the vacuum. The
latter prohibits populating a macroscopic number of low-momentum photons,
which is necessary to induce the transition to the localized phase Spohn and
Dümcke (1985); this is reminiscent of what has been discussed in the context
of the no-go theorems of the superradiant transition Rzażewski  _et al._
(1975); Nataf and Ciuti (2010); Viehmann  _et al._ (2011); De Bernardis  _et
al._ (2018); Andolina  _et al._ (2019); Stokes and Nazir (2020). Nevertheless,
we emphasize that our analysis here does not disprove the well-established
results in the standard spin-boson Hamiltonian, but rather point out the
crucial role of the diamagnetic term in the QED systems, which was not
manifestly taken into account in the original spin-boson model.

#### Linear

Finally, the special case of $l\\!=\\!1$ corresponds to the so-called Ohmic
bath, which is a rather marginal case due to the infrared logarithmic
divergence. On one hand, one may argue that the well-known mapping between the
Ohmic spin-boson model and the anisotropic Kondo model Blume  _et al._ (1970);
Leggett  _et al._ (1987) can be used in this case, from which the Kosterlitz-
Thouless-type transition from the delocalized to localized phase is expected
at a finite light-matter coupling. On the other hand, however, the validity of
such argument needs to be carefully reexamined in the present context of the
nonperturbative waveguide QED, because the relevance of the diamagnetic term
and the possible breakdown of the two-level description may make an important
difference between the usual spin-boson model and the full QED Hamiltonian.
Our preliminary analysis indicates that both $m_{\rm eff}$ and $\xi$ exhibit
the logarithmic divergence. Thus, (at least for the double-well potential) we
expect the situation to be analogous to what has been discussed in the sub-
Ohmic case above, but with the suppressed parameter region for the localized
phase. Nevertheless, the issue of how the known results of the Ohmic spin-
boson model can be related to the full-fledged QED Hamiltonian with a linear
dispersion is an ongoing subject especially in the context of the resistively
shunted Josephson junction Murani  _et al._ (2020); Kaur  _et al._ (2020); we
leave its full analysis to a future publication.

###  VI.2 Superradiant phase transition in multi-emitter cases

Our arguments above can be generalized to the case of multi-emitter systems
with large, but finite $N$. First, when all the emitters are collectively
coupled to the waveguide at the same position as discussed in Sec. V.3, the
analysis of the ground-state physics can be reduced to that of the collective-
mode Hamiltonian (79). One can then follow the same arguments above by
replacing the single-emitter effective mass $m_{\rm eff}$ by the multi-emitter
one $M_{\rm eff}$ in Eq. (80). Namely, as long as the renormalized mass
remains finite $\lim_{L\to\infty}M_{{\rm eff}}\\!<\\!\infty$, the ground state
should be unique with vanishing order parameter, which leads to the disordered
phase. This occurs when the photonic dispersion is either gapped or gapless
but with the sublinear exponent $l\\!<\\!1$. In particular, this conclusion is
consistent with the known no-go theorems of the superradiant phase transition
in cavity QED systems that accompany only a single (gapped) photonic mode. In
this context, our analysis greatly extends the applicability of those no-go
theorems to the realm of a multimode waveguide QED system with a
gapped/sublinear photonic dispersion.

Remarkably, in contrast, when the dispersion is gapless with superlinear
exponent $l\\!>\\!1$, transition to the superradiant phase with nonvanishing
order parameter is expected to occur. The reason is that, as discussed above,
the effective mass $M_{\rm eff}$ diverges polynomially in the thermodynamic
limit $L\\!\to\\!\infty$ in this case and, following the same arguments below
Eq. (95), we conclude that the ground state should be two-fold degenerate at a
sufficiently large $v$, i.e., at small qubit frequencies $\Delta$ (see Fig.
8). To our knowledge, this provides yet another pathway to realize the genuine
superradiant transition originally pursued along the line of the Dicke model.
We expect that our conclusion remains valid also when emitters are placed at
different positions, since the mass divergence originates from the infrared
divergence in the infinite-wavelength limit, where emitter separation should
not play a significant role.

Meanwhile, we note that taking the limit of many emitters $N\\!\to\\!\infty$
provides alternative way to realize a quantum phase transition. Importantly,
this is the case even when the dispersion is gapped/sublinear for which
transition is not expected at finite $N$ as discussed above. As detailed in
Sec. V.3, one can use the two-level effective model (84) to determine the
ground-state properties in these cases. This mapping to the transverse-field
Ising model immediately indicates the possibilities of inducing a transition
between the disordered (i.e., delocalized) phase and the ferromagnetic ordered
(i.e., localized) phase or realizing exotic phases such as the many-body
localized phase.

##  VII Summary and Discussions

We analyzed equilibrium and dynamical properties of light-matter systems
consisting of quantum emitters strongly interacting with quantized
electromagnetic continuum in the nonperturbative regimes, including the
previously unexplored deep and extremely strong coupling regimes. There,
traditional theoretical approaches utilizing the Coulomb or PZW gauges are no
longer sufficient, since substantial light-matter entanglement invalidates
truncations of emitter/photon levels in these gauges. We resolved this problem
by using the unitary transformation (9) that asymptotically disentangles
emitter and photon degrees of freedom in the strong-coupling limit; this new
frame of reference then enabled us to construct an accurate theoretical
framework at any finite interaction strengths. Our formulation puts no
limitations on the light-matter coupling and allows for consistently analyzing
the entire coupling region of waveguide QED systems at both qualitative and
quantitative levels.

We first applied this new formalism to the single-emitter system (see Eq.
(14)), and elucidated the essential features in the nonperturbative regimes on
the basis of general arguments. In particular, we demonstrated the emergence
of a ladder of many-body bound states and the (quasi) BIC, the vacuum
fluctuations induced suppression of potential barrier, and the strong
renormalization of the effective mass. We then analyzed these nonperturbative
features in a concrete model of cavity-array waveguide. All of these results
have relevance to ongoing experiments in superconducting qubits interacting
with microwave resonators or atoms coupled to photonic crystals. We proposed
that the BIC can experimentally be observed by analyzing nonequilibrium
dynamics induced by the quench of a parameter in the waveguide QED
Hamiltonian. This protocol should be implemented in circuit QED systems using
currently available experimental techniques. The parameter regimes we studied
are either directly relevant to state-of-the-art experimental systems in,
e.g., superconducting devices Yoshihara  _et al._ (2017); Forn-Díaz  _et al._
(2017) and plasmonic crystals Mueller  _et al._ (2020) or (at least) expected
to be accessible in the near future in view of rapid developments in achieving
stronger light-matter coupling regimes Forn-Díaz  _et al._ (2019); Sheremet
_et al._ (2021). To explore those nonperturbative regimes in setups of atoms
coupled to photonic crystals, one can utilize Rydberg atoms Agarwal (1984);
Bienias  _et al._ (2014); Liang  _et al._ (2018) and/or the collective
$\sqrt{N}$ enhancement of the light-matter coupling by assembling a large
number of weakly coupled components Roy  _et al._ (2017) as discussed in Sec.
V.3. We envision that the predicted tightly localized bound states in the
waveguides can be used as a photon storage in quantum information
applications.

We next extended the analysis to the multi-emitter case and established a
general framework for studying multi-emitter QED systems without relying on
uncontrolled approximations or assumptions. Building on this formalism, we
showed that the ground-state physics can be understood from the perspective of
the standard transverse-field Ising model (84) but with suitably renormalized
parameters. Finally, we analyzed the case of gapless photonic dispersions and
argued that a quantum phase transition can in general occur if the
renormalized mass diverges in the thermodynamic limit, while transition is not
expected when the mass remains finite. In certain special limits, this
analysis led to the conclusions that are consistent with previous results of
the delocalization-localization transition in the spin-boson model and the no-
go theorems for the superradiant transition in cavity QED. In this respect,
one of the conceptual advantages in our approach is establishing a unified
perspective on these seemingly unrelated phenomena. Most interestingly, we
pointed out the possibility of realizing the superradiant phase by coupling
(artificial) atoms to electromagnetic modes with superlinear dispersion, which
should provide a new twist to ongoing discussions about the possibility of the
superradiant transition in multimode QED setups (see e.g., Ref. Andolina  _et
al._ (2020)).

Several open questions remain for future studies. First, it merits further
study to figure out how losses and decoherence either for emitters or photons
can affect the present results. While the waveguide coupling efficiency can be
made close to the unit fidelity in, e.g., superconducting devices Mirhosseini
_et al._ (2019), losses are still ubiquitous in atoms coupled to photons in
optical domains Goban  _et al._ (2015). These issues can be addressed by
combining the present nonperturbative QED formalism with the standard
framework of Markovian open systems Ashida  _et al._ (2020b). Second, it will
be intriguing to analyze the multi-emitter QED systems in more detail,
especially in view of its mapping to the inhomogeneous long-range Ising model.
Third, instead of the exact diagonalization performed here, one can apply more
efficient numerical methods, such as the matrix-product-states calculations
Peropadre  _et al._ (2013); Sanchez-Burillo  _et al._ (2014); Shi  _et al._
(2018); Barkemeyer  _et al._ (2021) or a hybrid variational approach Wang  _et
al._ (2020), to analyze the asymptotically decoupled QED Hamiltonian within
the few-photon ansatz (31). This should be particularly useful when one is
interested in a larger system with many emitters being coupled to common
multiple electromagnetic modes.

Finally, while the emphasis was placed on the waveguide QED in this paper, our
theory is equally applicable to cavity QED setups with multiple photonic modes
(see e.g., Ref. Vaidya  _et al._ (2018)) whose inclusion is often important
depending on the cavity geometry and the coupling strength. In this respect,
our study is also relevant to a variety of strongly coupled light-matter
systems recently realized by using both solid state and quantum chemistry
platforms. For instance, the predicted vacuum-induced suppression of the
potential barrier may lie at the heart of the enhanced chemical reactivity
observed in polaritonic chemistry Hiura  _et al._ (2018); Hiura and Shalabney
(2019); Thomas  _et al._ (2019) as discussed in Sec. III.3. More generally,
our study reveals that the mass enhancement is one of the universal features
of strongly interacting light-matter systems. This naturally associates with
the higher density of states, which could lead to enhancements of certain
many-body properties, including superconductivity or ferromagnetism Thomas
_et al._ (2019, 2021); Basov  _et al._ (2016); Curtis  _et al._ (2019);
Schlawin  _et al._ (2019); Sentef  _et al._ (2018); Kiffner  _et al._ (2019);
Li and Eckstein (2020); Cottet  _et al._ (2015).

###### Acknowledgements.

We are grateful to Andrea Cavalleri, Ignacio Cirac, Jerome Faist, Victor
Galitski, Dieter Jaksch, Israel Klich, Angel Rubio, Hossein Sadeghpour, and
Tao Shi for useful discussions. Y.A. acknowledges support from the Japan
Society for the Promotion of Science through Grant No. JP19K23424. E.D.
acknowledges support from Harvard-MIT CUA, AFOSR-MURI Photonic Quantum Matter
(Grant No. FA95501610323), the ARO grant “Control of Many-Body States Using
Strong Coherent Light-Matter Coupling in Terahertz Cavities”, the NSF EAGER-
QAC-QSA Grant No. 2038011 “Quantum Algorithms for Correlated Electron-Phonon
System”, and the hospitality of the physics department at ETH Zurich.

##  Appendix A Diagonalization of the quadratic photon Hamiltonian

Here we provide details about the diagonalization of the quadratic photon
Hamiltonian including the $\hat{A}^{2}$ term in the Coulomb gauge. To this
end, we introduce the conjugate pairs of variables via

$\displaystyle\hat{X}_{k}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\hbar}{2\omega_{k}}}\left(\hat{a}_{k}+\hat{a}_{k}^{\dagger}\right),$
(96) $\displaystyle\hat{P}_{k}$ $\displaystyle=$ $\displaystyle
i\sqrt{\frac{\hbar\omega_{k}}{2}}\left(\hat{a}_{k}^{\dagger}-\hat{a}_{k}\right),$
(97)

and rewrite the quadratic part of $\hat{H}_{\rm C}$ in Eq. (1) as

$\displaystyle\frac{q^{2}\hat{A}^{2}}{2m}+\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}$
(98) $\displaystyle=$
$\displaystyle\sum_{k}\frac{\hat{P}_{k}^{2}}{2}\\!+\\!\frac{1}{2}\sum_{kk^{\prime}}\left(\delta_{kk^{\prime}}\omega_{k}^{2}\\!+\\!2g_{k}g_{k^{\prime}}\right)\hat{X}_{k}\hat{X}_{k^{\prime}}.$

The last term in Eq. (98) can readily be diagonalized by an orthogonal
transformation,

$\hat{P}_{k}=\sum_{n}O_{kn}\hat{\tilde{P}}_{n},\;\;\hat{X}_{k}=\sum_{n}O_{kn}\hat{\tilde{X}}_{n},$
(99)

where an orthogonal matrix $O$ satisfies Eq. (7). The quadratic photon
Hamiltonian then becomes

$\frac{1}{2}\sum_{n}\left(\hat{\tilde{P}}_{n}^{2}+\Omega_{n}^{2}\hat{\tilde{X}}_{n}^{2}\right)=\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n},$
(100)

where we introduce the annihilation operator after the orthogonal
transformation by

$\hat{b}_{n}=\sqrt{\frac{\Omega_{n}}{2\hbar}}\hat{\tilde{X}}_{n}+\frac{i}{\sqrt{2\hbar\Omega_{n}}}\hat{\tilde{P}}_{n},$
(101)

which gives Eq. (5). Using these squeezed photon operators $\hat{b}_{n}$, the
total Hamiltonian in the Coulomb gauge is now given by Eq. (II.1) in the main
text.

##  Appendix B Numerical diagonalization

Figure 9:  Convergence of the excitation energies against the total photon-
number cutoff $N_{c}$ at different coupling strengths $g$. The results are
obtained by the exact diagonalization of the transformed Hamiltonian (107)
within the few-photon many-body ansatz (31) with the total photon-number
cutoff $N_{c}$. Parameters are $J=0.2$, $v=0.5$, and $d=0.87$.

We describe details about the method used in Sec. IV to numerically
diagonalize the QED Hamiltonian in the AD frame. We begin with folding the
electromagnetic modes of the cavity array onto even and odd modes with respect
to the spatial parity. Since only the even modes,

$\hat{a}_{x=0}^{e}=\hat{a}_{x=0},\;\;a_{x>0}^{e}=\frac{1}{\sqrt{2}}(\hat{a}_{x}+\hat{a}_{-x}),$
(102)

interact with the emitter, we neglect contributions from the odd modes. The
photon Hamiltonian in the real-space basis is then written as

$\displaystyle\hat{H}_{{\rm light}}$
$\displaystyle=-\frac{J}{\sqrt{2}}\left(\hat{a}_{x=0}^{\dagger
e}\hat{a}_{x=1}^{e}+{\rm H.c.}\right)$ (103)
$\displaystyle-\frac{J}{2}\sum_{x=1}^{(L-1)/2}\left(\hat{a}_{x+1}^{\dagger
e}\hat{a}_{x}^{e}+{\rm
H.c.}\right)\\!+\\!\hbar\omega_{c}\sum_{x=0}^{(L-1)/2}\hat{a}_{x}^{\dagger
e}\hat{a}_{x}^{e}$
$\displaystyle\equiv\sum_{p=0}^{(L-1)/2}\hbar\omega_{p}\hat{a}_{p}^{\dagger}\hat{a}_{p},$

where we use an orthogonal matrix $M$ to obtain the diagonalized form (103)
with

$\hat{a}_{x}^{e}=\sum_{p=0}^{(L-1)/2}M_{xp}\hat{a}_{p}.$ (104)

The vector field is expressed in this basis as

$\displaystyle\hat{A}_{x=0}$ $\displaystyle=$ $\displaystyle{\cal
A}(\hat{a}_{x=0}+\hat{a}_{x=0}^{\dagger})$ (105) $\displaystyle=$
$\displaystyle\sum_{p=0}^{(L-1)/2}{\cal
A}M_{0p}\left(\hat{a}_{p}+\hat{a}_{p}^{\dagger}\right),$

which results in the electromagnetic amplitudes (cf. Eq. (3))

$f_{p}={\cal A}M_{0p}.$ (106)

We use the frequencies $\omega_{p}$ and the amplitudes $f_{p}$ to diagonalize
the photon Hamiltonian including the $\hat{A}^{2}$ term as explained in
Appendix A. Then, we perform the unitary transformation (9) and arrive at the
Hamiltonian in the AD frame

$\displaystyle\hat{H}_{U}$ $\displaystyle=$
$\displaystyle\frac{\hat{P}^{2}}{2m_{{\rm eff}}}+V_{{\rm eff}}(\hat{Q})$ (107)
$\displaystyle+$
$\displaystyle\sum_{l=1}\frac{:\hat{\Xi}^{l}:}{l!}V^{(l)}(\hat{Q})+\\!\sum_{n=0}^{(L-1)/2}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n}.$

To numerically diagonalize this Hamiltonian efficiently, we employ the few-
photon many-body ansatz (31) and express the matrix elements of $\hat{H}_{U}$
in terms of the basis

$|\psi_{\alpha}\rangle_{\rm emitter}\otimes|n_{0}n_{1}\cdots
n_{(L-1)/2}\rangle_{\rm photon}$ (108)

with level truncations,

$\alpha=1,2,\ldots,\alpha_{c},\;\;\sum_{j=0}^{(L-1)/2}n_{j}\leq N_{c}.$ (109)

We recall that $|\psi_{\alpha}\rangle$ are single-particle eigenstates of the
renormalized emitter Hamiltonian (22), and $|n_{0}n_{1}\cdots
n_{(L-1)/2}\rangle$ is a many-body bosonic Fock state in terms of
$\hat{b}_{n}$ operators with $n_{j}=0,1,\cdots$. The corresponding Hilbert-
space dimension is

$D=\alpha_{c}\sum_{i=1}^{N_{c}}[(L+1)/2]^{i}=\frac{\alpha_{c}([(L+1)/2]^{N_{c}+1}-1)}{(L+1)/2-1},$
(110)

which grows polynomially with the system size $L$. Figure 9 demonstrates that
the numerical results converge very efficiently with $N_{c}$ in a broad range
of the light-matter coupling strength. It typically suffices to set
$N_{c}=2$-$4$ and $\alpha_{c}={\rm O}(10)$ to achieve the accuracy with an
error below $\sim\\!1\%$. When only the low-energy spectrum is of interest,
one can use the Lanczos method to further reduce the computational cost.

##  Appendix C Derivation of the multi-emitter Hamiltonian in the
asymptotically decoupled frame

We here derive the asymptotically decoupled multi-emitter Hamiltonian
discussed in Sec. V. As we have done for the single-emitter case, we first
diagonalize the quadratic photon part of the Coulomb-gauge Hamiltonian
including the $\hat{A}^{2}$ term. To do so, we introduce the position-
dependent multi-emitter coupling strengths by

$\displaystyle g_{kj}^{c}$ $\displaystyle=$ $\displaystyle
qf_{kj}\sqrt{\frac{\omega_{k}}{m_{j}\hbar}}\cos(kx_{j}),$ (111) $\displaystyle
g_{kj}^{s}$ $\displaystyle=$ $\displaystyle
qf_{kj}\sqrt{\frac{\omega_{k}}{m_{j}\hbar}}\sin(kx_{j}),$ (112)

and rewrite the quadratic part as (aside constant)

$\displaystyle\sum_{j}\frac{q^{2}\hat{A}_{x_{j}}^{2}}{2m_{j}}+\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}$
(113) $\displaystyle=$
$\displaystyle\sum_{kk^{\prime}}\left(\delta_{kk^{\prime}}+2\sum_{j}\frac{g_{kj}^{s}g_{k^{\prime}j}^{s}}{\omega_{k}\omega_{k^{\prime}}}\right)\frac{\hat{P}_{k}\hat{P}_{k^{\prime}}}{2}$
$\displaystyle+\sum_{kk^{\prime}}\left(\delta_{kk^{\prime}}\omega_{k}^{2}+2\sum_{j}g_{kj}^{c}g_{k^{\prime}j}^{c}\right)\frac{\hat{X}_{k}\hat{X}_{k^{\prime}}}{2}$
$\displaystyle-\sum_{kk^{\prime}}\sum_{j}\frac{g_{kj}^{c}g_{k^{\prime}j}^{s}}{\omega_{k^{\prime}}}\left(\hat{X}_{k}\hat{P}_{k^{\prime}}+\hat{P}_{k^{\prime}}\hat{X}_{k}\right),$

where we recall that the conjugate operators $\hat{X}_{k}$ and $\hat{P}_{k}$
are defined by Eqs. (96) and (97). Equation (113) can then be diagonalized by
the symplectic transformation,

$\left(\begin{array}[]{c}\hat{\bm{X}}\\\
\hat{\bm{P}}\end{array}\right)=S\left(\begin{array}[]{c}\hat{\tilde{\bm{X}}}\\\
\hat{\tilde{\bm{P}}}\end{array}\right)\equiv\left(\begin{array}[]{cc}S^{XX}&S^{XP}\\\
S^{PX}&S^{PP}\end{array}\right)\left(\begin{array}[]{c}\hat{\tilde{\bm{X}}}\\\
\hat{\tilde{\bm{P}}}\end{array}\right),$ (114)

where the matrix $S$ satisfies

$\displaystyle S\sigma S^{\rm T}=\sigma$ (115)

with $\sigma=i\sigma^{y}\otimes{\rm I}_{L}$ and ${\rm I}_{L}$ being the
$L\times L$ identity matrix. This leads to the diagonalized form

$\displaystyle\sum_{j}\frac{q^{2}\hat{A}_{x_{j}}^{2}}{2m_{j}}+\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n}\left(\hat{\tilde{P}}_{n}^{2}+\Omega_{n}^{2}\hat{\tilde{X}}_{n}^{2}\right)$
(116) $\displaystyle=$
$\displaystyle\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n},$

where we define the squeezed photon operators $\hat{b}_{n}$ in the same manner
as in Eq. (101).

In terms of these new photon operators, the Coulomb-gauge Hamiltonian is
expressed as

$\displaystyle\hat{H}_{{\rm C}}$ $\displaystyle=$
$\displaystyle\sum_{j}\left[\frac{\hat{P}_{j}^{2}}{2m_{j}}+V(\hat{Q}_{j})\right]$
(117)
$\displaystyle-\sum_{jn}\hat{P}_{j}\left(\zeta_{nj}^{*}\hat{b}_{n}+\zeta_{nj}\hat{b}_{n}^{\dagger}\right)+\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n},$

where we define

$\displaystyle\zeta_{nj}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\hbar}{m_{j}\Omega_{n}}}$ $\displaystyle\times$
$\displaystyle\sum_{k}\left[g_{kj}^{c}S_{kn}^{XX}-g_{kj}^{s}S_{kn}^{PX}+i\left(g_{kj}^{c}S_{kn}^{XP}-g_{kj}^{s}S_{kn}^{PP}\right)\right].$

We now introduce the multi-emitter extension of the asymptotically decoupling
unitary transformation by

$\displaystyle\hat{U}$ $\displaystyle=$
$\displaystyle\exp\Biggl{[}\frac{1}{\hbar}\sum_{jn}\hat{P}_{j}\frac{(\zeta_{nj}\hat{b}_{n}^{\dagger}-\zeta_{nj}^{*}\hat{b}_{n})}{\Omega_{n}}\Biggr{]}$
(119) $\displaystyle\equiv$
$\displaystyle\exp\Bigl{(}-\frac{i}{\hbar}\sum_{j}\hat{P}_{j}\hat{\Xi}_{j}\Bigr{)},$

where we recall
$\hat{\Xi}_{j}=\sum_{n}i(\xi_{nj}\hat{b}_{n}^{\dagger}-\xi_{nj}^{*}\hat{b}_{n})$
with the displacement parameters

$\displaystyle\xi_{nj}=\frac{\zeta_{nj}}{\Omega_{n}}.$ (120)

This transformation acts on the photon and emitter operators as

$\displaystyle\hat{U}^{\dagger}\hat{b}_{n}\hat{U}$ $\displaystyle=$
$\displaystyle\hat{b}_{n}+\sum_{j}\frac{\xi_{nj}\hat{P}_{j}}{\hbar},$ (121)
$\displaystyle\hat{U}^{\dagger}\hat{Q}_{j}\hat{U}$ $\displaystyle=$
$\displaystyle\hat{Q}_{j}+\hat{\Xi}_{j},$ (122)

and transforms Eq. (117) to the following form:

$\displaystyle\hat{H}_{U}$ $\displaystyle=$
$\displaystyle\hat{U}^{\dagger}\hat{H}_{{\rm C}}\hat{U}$ (123)
$\displaystyle=$ $\displaystyle\sum_{j}\left[\frac{\hat{P}_{j}^{2}}{2m_{{\rm
eff},j}}+V(\hat{Q}_{j}+\hat{\Xi}_{j})\right]-\sum_{i>j}\mu_{ij}\hat{P}_{i}\hat{P}_{j}$
$\displaystyle+\sum_{n}\hbar\Omega_{n}\hat{b}_{n}^{\dagger}\hat{b}_{n},$

which gives Eq. (66) in the main text. Here, we introduce the effective mass
$m_{{\rm eff},j}$ for each emitter and the emitter-emitter coupling $\mu_{ij}$
as

$\displaystyle m_{{\rm eff},j}$ $\displaystyle\equiv$
$\displaystyle\frac{m_{j}}{1-\sum_{n}\frac{2m_{j}|\zeta_{nj}|^{2}}{\hbar\Omega_{n}}},$
(124) $\displaystyle\mu_{ij}$ $\displaystyle\equiv$
$\displaystyle\sum_{n}\frac{2\zeta_{ni}^{*}\zeta_{nj}}{\hbar\Omega_{n}}.$
(125)

The expressions (68) and (69) in Sec. V follow from the relations (116) and
(C).

## References

  * Rabi (1937) I. I. Rabi,  _Space quantization in a gyrating magnetic field,_ Phys. Rev. 51, 652 (1937).
  * Cohen-Tannoudji  _et al._ (1989) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms_ (Wiley, New York, 1989).
  * Bishop  _et al._ (2009) L. S. Bishop, J. Chow, J. Koch, A. Houck, M. Devoret, E. Thuneberg, S. Girvin, and R. Schoelkopf,  _Nonlinear response of the vacuum rabi resonance,_ Nat. Phys. 5, 105 (2009).
  * Astafiev  _et al._ (2010) O. Astafiev, A. M. Zagoskin, A. A. Abdumalikov, Y. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai,  _Resonance fluorescence of a single artificial atom,_ Science 327, 840 (2010).
  * Hoi  _et al._ (2011) I.-C. Hoi, C. M. Wilson, G. Johansson, T. Palomaki, B. Peropadre, and P. Delsing,  _Demonstration of a single-photon router in the microwave regime,_ Phys. Rev. Lett. 107, 073601 (2011).
  * Hoi  _et al._ (2012) I.-C. Hoi, T. Palomaki, J. Lindkvist, G. Johansson, P. Delsing, and C. M. Wilson,  _Generation of nonclassical microwave states using an artificial atom in 1D open space,_ Phys. Rev. Lett. 108, 263601 (2012).
  * van Loo  _et al._ (2013) A. F. van Loo, A. Fedorov, K. Lalumière, B. C. Sanders, A. Blais, and A. Wallraff,  _Photon-mediated interactions between distant artificial atoms,_ Science 342, 1494 (2013).
  * Mlynek  _et al._ (2014) J. A. Mlynek, A. A. Abdumalikov, C. Eichler, and A. Wallraff,  _Observation of dicke superradiance for two artificial atoms in a cavity with high decay rate,_ Nat. Commun. 5, 1 (2014).
  * Mirhosseini  _et al._ (2019) M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang, and O. Painter,  _Cavity quantum electrodynamics with atom-like mirrors,_ Nature 569, 692 (2019).
  * Kannan  _et al._ (2020) B. Kannan, M. J. Ruckriegel, D. L. Campbell, A. F. Kockum, J. Braumüller, D. K. Kim, M. Kjaergaard, P. Krantz, A. Melville, B. M. Niedzielski, A. Vepsalainen, R. Winik, J. L. Yoder, F. Nori, T. P. Orlando, and W. D. Gustavsson, S.and Oliver,  _Waveguide quantum electrodynamics with superconducting artificial giant atoms,_ Nature 583, 775 (2020).
  * Reitz  _et al._ (2013) D. Reitz, C. Sayrin, R. Mitsch, P. Schneeweiss, and A. Rauschenbeutel,  _Coherence properties of nanofiber-trapped cesium atoms,_ Phys. Rev. Lett. 110, 243603 (2013).
  * Thompson  _et al._ (2013) J. D. Thompson, T. G. Tiecke, N. P. de Leon, J. Feist, A. V. Akimov, M. Gullans, A. S. Zibrov, V. Vuletić, and M. D. Lukin,  _Coupling a single trapped atom to a nanoscale optical cavity,_ Science 340, 1202 (2013).
  * Arcari  _et al._ (2014) M. Arcari, I. Söllner, A. Javadi, S. Lindskov Hansen, S. Mahmoodian, J. Liu, H. Thyrrestrup, E. H. Lee, J. D. Song, S. Stobbe, and P. Lodahl,  _Near-unity coupling efficiency of a quantum emitter to a photonic crystal waveguide,_ Phys. Rev. Lett. 113, 093603 (2014).
  * Yalla  _et al._ (2014) R. Yalla, M. Sadgrove, K. P. Nayak, and K. Hakuta,  _Cavity quantum electrodynamics on a nanofiber using a composite photonic crystal cavity,_ Phys. Rev. Lett. 113, 143601 (2014).
  * Goban  _et al._ (2014) A. Goban, C.-L. Hung, S.-P. Yu, J. Hood, J. Muniz, J. Lee, M. Martin, A. McClung, K. Choi, D. E. Chang, O. Painter, and J. Kimble,  _Atom-light interactions in photonic crystals,_ Nat. Commun. 5, 1 (2014).
  * Lodahl  _et al._ (2015) P. Lodahl, S. Mahmoodian, and S. Stobbe,  _Interfacing single photons and single quantum dots with photonic nanostructures,_ Rev. Mod. Phys. 87, 347 (2015).
  * 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).
  * Schmid (1983) A. Schmid,  _Diffusion and localization in a dissipative quantum system,_ Phys. Rev. Lett. 51, 1506 (1983).
  * Guinea  _et al._ (1985) F. Guinea, V. Hakim, and A. Muramatsu,  _Diffusion and localization of a particle in a periodic potential coupled to a dissipative environment,_ Phys. Rev. Lett. 54, 263 (1985).
  * Weiss (2012) U. Weiss,  _Quantum dissipative systems_ , Vol. 13 (World scientific, 2012).
  * 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,_ Angew. Chem. Int. Ed. 51, 1592 (2012).
  * Galego  _et al._ (2015) J. Galego, F. J. Garcia-Vidal, and J. Feist,  _Cavity-induced modifications of molecular structure in the strong-coupling regime,_ Phys. Rev. X 5, 041022 (2015).
  * Ebbesen (2016) T. W. Ebbesen,  _Hybrid light–matter states in a molecular and material science perspective,_ Acc. Chem. Res. 49, 2403 (2016).
  * Herrera and Spano (2016) F. Herrera and F. C. Spano,  _Cavity-controlled chemistry in molecular ensembles,_ Phys. Rev. Lett. 116, 238301 (2016).
  * 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,_ Proc. Natl. Acad. Sci. U.S.A. 114, 3026 (2017).
  * Feist  _et al._ (2017) J. Feist, J. Galego, and F. J. Garcia-Vidal,  _Polaritonic chemistry with organic molecules,_ ACS Photonics 5, 205 (2017).
  * Hiura  _et al._ (2018) H. Hiura, A. Shalabney, and J. George,  _Cavity Catalysis ‒Accelerating Reactions under Vibrational Strong Coupling,_ ChemRxiv (2018), 10.26434/chemrxiv.7234721.v3.
  * Hiura and Shalabney (2019) H. Hiura and A. Shalabney,  _A Reaction Kinetic Model for Vacuum-Field Catalysis Based on Vibrational Light-Matter Coupling,_ ChemRxiv (2019), 10.26434/chemrxiv.9275777.v1.
  * Thomas  _et al._ (2019) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen,  _Tilting a ground-state reactivity landscape by vibrational strong coupling,_ Science 363, 615 (2019).
  * Haugland  _et al._ (2021) T. S. Haugland, C. Schaefer, E. Ronca, A. Rubio, and H. Koch,  _Intermolecular interactions in optical cavities: An ab initio qed study,_ J. Chem. Phys. 154, 094113 (2021).
  * Barnes  _et al._ (2003) W. L. Barnes, A. Dereux, and T. W. Ebbesen,  _Surface plasmon subwavelength optics,_ Nature 424, 824 (2003).
  * Pitarke  _et al._ (2006) J. M. Pitarke, V. M. Silkin, E. V. Chulkov, and P. M. Echenique,  _Theory of surface plasmons and surface-plasmon polaritons,_ Rep. Prog. Phys. 70, 1 (2006).
  * Tame  _et al._ (2013) M. S. Tame, K. McEnery, Ş. Özdemir, J. Lee, S. Maier, and M. Kim,  _Quantum plasmonics,_ Nat. Phys. 9, 329 (2013).