Nonperturbative Waveguide Quantum Electrodynamics

Understanding physical properties of quantum emitters strongly interacting with quantized electromagnetic modes 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 on the basis of a comprehensive and rigorous theoretical approach using an asymptotic decoupling unitary transformation. 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. To illustrate our general results with concrete examples, we use our formalism to study 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 a reentrant transition can occur in the regimes where the coupling strength becomes the dominant energy scale. We also discuss applications of our results to other problems in different fields, including quantum optics, condensed matter physics, and quantum chemistry.


A. Background
Quantum states arising from strong coherent interaction between light and matter are not only interesting from the perspective of fundamental many-body physics, but also provide promising new platforms for quantum technologies. Historically, analysis of light-matter systems focused on the perturbative regime [1,2], since interaction of atomic dipoles with vacuum electromagnetic fields is weak due to the smallness of the fine structure constant α = 1/137. Recent progress has led to experimental realizations of new systems in which electromagnetic field is modified to reach stronger light-matter coupling. In particular, (artificial) atoms coupled to onedimensional continuum of photons at microwave [3][4][5][6][7][8][9][10][11] or optical [12][13][14][15][16][17] frequencies achieve enhancement of lightmatter 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 [18][19][20][21]; 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 [22][23][24][25][26][27][28][29][30][31] and nanostructured plasmonics [32][33][34][35][36]. In light of such broad relevance, models of quantum emitters interacting with a continuum of bosonic excitations have played a crucial role in quantum information science as well as 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 among 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 [2] due to the processes that create or annihilate pairs of excitations. Moreover, the inclusion of the diamagneticÂ 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 [37] and such multilevel structures are potentially useful for the purpose of performing certain quantum information operations [38,39]. While significant efforts have been devoted to elucidating the strong and ultrastrong coupling regimes in the last decade , the physics of waveguide QED in the realms of even stronger light-matter interactions, namely, the deep [64] and extremely strong [65] 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 these intriguing phenomena is currently lacking. The primary goal of this paper is to reveal physics of strongly interacting light-matter systems in the previously unexplored regimes on the basis of a comprehensive theoretical framework that avoids problems discussed above.
On another front, the spin-boson model, a supposedly effective description of waveguide QED systems (e.g., Ref. [45]), has long been known to exhibit the delocalization-localization transition at strong couplings [66,67]. Nevertheless, the breakdown of the usual two-level description [65,68] and the relevance of the diamagnetic term [65,[68][69][70][71][72][73][74][75] have made it unclear until now 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 [76][77][78].
In view of recent experimental developments in realizing stronger light-matter interactions [36,51,79], 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 delocalization-localization transition in waveguide QED setups, and if so, does there exist a new feature?
The main aim of this paper is to reveal the new physics and develop understanding of strongly interacting light-matter systems by addressing these questions from a unified perspective. Below we summarize the main results at a nontechnical level before presenting a detailed theoretical formulation in subsequent sections.

B. Summary of the main results
Our first main result is the appearance of a ladder of manybody bound states (BS) and the many-body bound states in the continuum (BIC) in nonperturbative regimes (see Fig. 1). We point out formation of increasingly many low-lying bound states whose energies decrease as ∝ g −1 in the limit of strong light-matter coupling g. The exact BICs emerge as a consequence of the Z 2 -symmetry that is linked to microscopic QED Hamiltonians. Previous studies have so far discussed the realizations of one-body BICs, which relied on artificial tuning of either emitter positions or resonator wavelengths/geometry [80][81][82][83][84][85]. In contrast, a new type of BICs found here does not rely on either of them, but emerge from strong light-matter interaction (without artificial fine tunings) and thus have a many-body origin. Even when the symmetry is not exact, the lifetime of these states diverges as ∝ g 3/2 in the strongcoupling 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 for realizing quantum memory [86] and nondissipative emitter interactions [87,88]. 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 the states in a model of cavity array [89][90][91], leading to rich nonequilibrium dynamics in which the bound states and the dynamical Casimir effect are intertwined (see Fig. 6 in Sec. IV). 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 [92] in several crucial aspects. In Refs. [49,50,58,[93][94][95][96][97][98][99], 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 [47,57]. 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 or BICs revealed by our analysis are appreciable only after multilevel structure of emitters is consistently included in theory. We also remark that these bound states are genuine quantum many-body states in contrast to one-body wave phenomena, which have been the main focus of earlier studies [100].
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, fol- 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 eff ∝ g −1/2 at large g. In contrast, the renormalized mass is enhanced as m eff ∝ 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 eff with lower potential barrier due to the dressing by the vacuum electromagnetic fluctuations. Note that, when going back to the Coulomb gauge, Q in the asymptotically decoupled frame contains both matter and light contributions. (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 ∝ g −1 and thus, these states become increasingly degenerate at strong couplings. We note that there also exist the extra degeneracy corresponding to the number of degenerate potential minima, for which the energy spacing closes exponentially as g is increased.
lowing 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. (53) in Sec. IV) and the inhomogeneous transverse-field Ising Hamiltonian for a multi-emitter case (see Eq. (85) in Sec. V), but with suitably renormalized parameters. These results 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, which is reminiscent of the delocalization-localization transition in the spinboson models. Our results also indicate a qualitatively new feature, not present in the simplified spin-boson descriptions, such as the reentrant transition into the delocalized phase in the extremely strong coupling regimes which originates from the mass acquisition in the transformed frame (see Figs. 8 and 9 in Sec. VI). We demonstrate these results by applying the functional renormalization group method to a concrete model of resistively shunted Josephson junctions.
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 can-not be justified in the nonperturbative regimes. Below we thus start by developing a rigorous framework for describing quantum emitters coupled to arbitrary multiple quantum electromagnetic modes, including the case of a continuum spectrum. This is done by extending the asymptotic light-matter decoupling unitary transformation that we introduced earlier in the context of single-mode cavity QED [65] to the present waveguide QED setups. Most of the previous studies approximated an emitter as a simplified two-level system, which, however, is not a valid approximation for many experimentally relevant systems, including superconducting qubits as mentioned before [39,[101][102][103]. 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 C). 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)).
While the emphasis of our discussion is on the waveguide setups, the present formalism can be extended 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 quan-tum 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 delocalization-localization 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 first develop a general theory of a single quantum emitter coupled to arbitrary quantized electromagnetic environment. 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 singleemitter setup is considered in this section, we will generalize the whole formalism to multi-emitter cases in Sec. V.

A. 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 ω k . The corresponding QED Hamiltonian in the Coulomb gauge is given bŷ whereQ (P ) is the position (momentum) operator of the emitter andâ k (â † k ) is the annihilation (creation) operator of photons in mode k, which satisfy the commutation relations We denote the vector potential operator aŝ where f k characterizes the electromagnetic amplitude of mode k.
It is useful to diagonalize the quadratic photon part ofĤ C as follows (see Appendix A for details): where we perform the canonical transformation to introduce a squeezed photon operatorb n labeled by n ∈ Z viâ and ζ n is given as Here, O kn is an orthogonal matrix that satisfies where Ω n is an eigenfrequency of mode n, r nk is a squeezing parameter defined by e r nk ≡ Ω n /ω k , and g k characterizes a coupling strength to mode k: We note that the magnitudes of g k depend on the size of the electromagnetic environment L via g k ∝ f k ∝ L −1/2 . In a concrete model discussed later (cf. Eq. (41)), the environment consists of the coupled cavity arrays and the variable L corresponds to the total number of cavities. 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 wellestablished analogy between circuit and atomic QED systems; the charge number operator of a transmon qubit and its conjugate phase operator precisely correspond toP andQ in Eq. (1), respectively, and the charge bias induced by the electromagnetic fields of microwave resonator plays the role of the vector potentialÂ (see also Sec. VI B). The same analogy holds true also for a flux qubit, whereQ is coupled to photons through the dipole-type couplingQ ·Ê withÊ being the electric field; one can use the Power-Zienau-Woolley (PZW) transformation [104,105] to change this circuit Hamiltonian back to the standard form as in Eq. (1) (see e.g., Ref. [68] or Eq. (40) below). In practice, coefficients of theÂ 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.
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 TABLE I. 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 ω 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 Ωn are determined by the eigenvalue problem (7). The label n = 0 indicates the dominant electromagnetic mode with the largest eigenfrequency. The length scales ξn in Eq. (11) are normalized by xω in Eq. (20) and characterize the effective light-matter coupling strengths for the dominant n = 0 and the other modes n = 0 in the transformed frame. The effective mass m 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 C); see Eqs. (30) and (38).
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 thê 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 theQ dependence of the vector potentialÂ 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.

B. Asymptotic decoupling transformation
We now introduce a unitary transformation to asymptotically decouple light and matter degrees of freedom [65]: whereΞ is given asΞ This transformation acts on individual operators viâ where the emitter position is shifted by the gauge-fielddependent displacementΞ while each photon mode is subject to the momentum-dependent shift ξ nP / [106]. We note that the displacement variables ξ n in Eq. (11) are chosen in such a way that theP ·(b+b † ) term in Eq. (4) will be precisely cancelled by the contributions arising from the displacement of theb †b term via Eq. (13). The resulting Hamiltonian in the asymptotically decoupled frame isĤ Here the effective mass is defined as where the mass enhancement is characterized by the dimensionless quantity Θ whose expression (16) follows from Eq. (7). This renormalization comes from theP 2 terms arising from the residual contributions generated by displacing theP · (b +b † ) andb †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 ξ n instead of the bare coupling g k .
Hereafter we first focus on the case of a gapped dispersion with frequencies ω k > 0 ∀k, for which Θ 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 Θ in that case; we will revisit this issue in Sec. VI.

C. 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 ω and the coupling strength g as follows: where we note g = O(L 0 ) since g k ∝ L −1/2 . We begin by considering the regime g/ω > 1 in which the coupling strength is dominant over other energy scales; we shall refer to it as the extremely strong coupling (ESC) regime [65]. There, the eigenvalue problem (7) has a single dominant mode (which we label n = 0) with the largest eigenfrequency Ω 0 ∝ g and . This leads to the scalings ξ 0 ∝ g −1/2 and ξ n =0 ∝ g −1 , i.e., the light-matter interaction inĤ U asymptotically vanishes in the strong-coupling limit.
In the deep strong coupling (DSC) regime g/ω ∼ 1 [64], one can continue the scaling analysis in the similar manner and obtain slightly refined expressions as summarized in Table I. There, we denote the variance of a photonic dispersion (or the effective bandwidth) as and normalize the length scale by the characteristic one In Table I, we also summarize the scaling relations in the ultrastrong coupling (USC) regime g/ω ∼ 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 ξ n and consequently the effective lightmatter 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.

D. Vacuum-dressed potential and decoupled excitations
To analyze low-energy eigenstates ofĤ U , it is useful to rewrite it in the following manner: where we define the matter Hamiltonian bŷ with m eff being the renormalized mass (15) and V eff being the dressed potential given as where V (l) is the l-th derivative of V . The interaction Hamiltonian is given byĤ where :Ô :≡Ô − 0|Ô|0 represents the normal ordering of photonic operators with |0 being the vacuum state in the AD frame:b We emphasize that this vacuum state is distinct from the original vacuum ofâ operators in the Coulomb gauge due to the squeezing (cf. Eq. (5)). Finally, we denote the photon Hamiltonian asĤ Equations (23) and (25) can be obtained by expanding the interaction term in Eq. (14) with respect toΞ and using the relations 0|Ξ 2l |0 = (2l − 1)!!ξ 2l and 0|Ξ 2l−1 |0 = 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 (Q+Ξ) renders the origin of the Lamb shift explicit; fluctuations inΞ add to the intrinsic fluctuations ofQ 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 ξ comes from the mode n = 0 and thus ξ basically obeys the same scaling relation satisfied by ξ 0 in Table I. 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 ξ. 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) and the effectivemass enhancement in Eq. (15).
Because of the decoupling ξ 0 ∝ g −1/2 and enhancement of the effective photon frequency Ω 0 ∝ g, low-energy eigenstates ofĤ U in the strong-coupling limit can be written as a product of the emitter eigenstates and the photon vacuum as follows: where |ψ α with α = 1, 2, . . . are single-particle eigenstates ofĤ matter in Eq. (22); we then represent it aŝ with E 1 ≤ E 2 ≤ · · · being the corresponding eigenenergies. These single-particle energies E α provide asymptotically exact excitation energies of the total HamiltonianĤ 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)). Said differently, when transforming back to the original Coulomb gauge, the above decoupled emitter states are in general entangled, light-matter correlated states. We note that the mass enhancement m eff ∝ g 2 leads to the tight localization of |ψ α around the bottom of V eff ; accordingly, the excitation energies δE α ≡ E α − E 1 decrease as δE α ∝ g −1 in the nonperturbative regimes as long as V eff has well-defined minima.

E. Few-photon 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 I). In particular, in the ESC regime, the standard perturbation theory predicts that the photon number in the ground state is on the order of (ξ 0 /Ω 0 ) 2 as long as the bandwidth δ is narrow, resulting in the scaling where · · · 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 = 0 modes eventually dominate the n = 0 contribution above, and the perturbation theory leads to and Ω n =0 ∼ ω (see Table I). This crossover occurs when the bandwidth reaches a value around gδ 4 /ω 5 = O(1), at which the contributions from n = 0 and n = 0 modes become comparable. In any case, the total photon number in the transformed frame asymptotically vanishes in the strong-coupling limit.
This fact motivates us to introduce the few-photon ansatz by projecting the whole Hilbert space onto the following subspace: where we recall that |ψ α are single-particle eigenstates of H matter in Eq. (22) while we denote |ψ photon,i as an arbitrary bosonic many-body state that satisfies 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 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 H Nc scales as ∝ L Nc . 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) ∼ 1% deviation already at a small total photon-number cutoff N c = 2-4 (see Appendix B for further details).
This point should be contrasted to previous approaches; eigenstates in the Coulomb gauge can possess large photon occupation numbers (see also Sec. III D below), and one has to include more excitations for each electromagnetic mode at a greater coupling g. Hence, the corresponding Hilbert-space dimension exponentially increases with L, which severely limits their applicabilities in the strong-coupling regions. Some variational states, such as the displaced-oscillator states [53,57,107,108], 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 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 now present 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 H 0 . The results discussed here establish universal nonperturvative features which hold true independent of specific choices or fine-tuning of microscopic parameters. We will make these predictions quantitatively accurate in the next section by extending the analysis to the few-photon ansatz in the subspace H Nc>0 .
A. 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 δE α ∝ 1/g due to the mass enhancement, and eventually lie outside of the photon continuum, where ω L(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 δE ∝ 1/g. This results from the increase of the emitter mass m eff ∝ 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 manybody BS, which are strongly entangled states including highmomentum 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 [92].

B. Many-body bound states in the continuum
Even when the decoupled excitations (28) lie within the photon continuum, they can behave as either symmetryprotected bound states in the continuum (BIC) or quasi BIC with lifetime that diverges in the strong-coupling limit. 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 Z 2 symmetry: whereP acts on the emitter operators asP −1PP = −P , P −1QP = −Q, and also transforms the photon field viâ P −1b nP = −b n . We emphasize that this Z 2 symmetry is intrinsically linked to microscopic light-matter Hamiltonians without making artificial fine tuning. For instance, in a circuit setup, the potential term V (Q) typically consists of the sum of the Josephson energy −E J cos(Q) and the inductive term E L Q 2 /2, both of which clearly satisfy the above symmetry.
At a more fundamental level, since the first-principle QED Hamiltonian in the Coulomb gauge also naturally satisfies this symmetry [2], the present QED Hamiltonian (1) should also respect that in general.
It is then 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 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 which vanishes in the strong-coupling limit. The same argument also applies to the case when the Z 2 symmetry is not exact due to, e.g., the broken inversion symmetry, V (Q) = 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) BICs 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 the mass renormalization factor Θ diverges (see Sec. V).
It is worthwhile to note again 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, The reason for this is as follows. At zero coupling, one can find an emitter state lying above the continuum, i.e., δE α (g = 0) > ω U . In the ESC regime, this excitation energy asymptotically decreases as δE α (g) ∝ 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.

C. Vacuum-induced suppression in potential barrier
Yet another common feature in the nonperturbative regimes is the vacuum-induced suppression of potential barrier in V 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 ξ that exhibits the nonmonotonic g dependence (see Eq. (24) and Table I).
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 [28][29][30], where the thermally activated chemical reaction was found to be enhanced due to cavity confinement.

D. 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 I 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 illjustified at a sufficiently large g in the conventional gauges; in particular, this fact indicates the breakdown of the usual twolevel 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 where we expand an emitter state |ψ α in terms of the momentum eigenstates |P . In the strong-coupling limit, the variance of the momentum distribution |ψ α (P )| 2 diverges with σ P ∝ g due to the mass renormalization m eff ∝ 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Ĥ PZW =Û † PZWĤ CÛPZW with the PZW transformationÛ PZW = exp(iqQÂ/ ): 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. We remark that, when transforming back to the Coulomb gauge, E in the above PZW gauge corresponds to the dielectric displacement field consisting of the electric field and the emitter shift.

IV. APPLICATION TO COUPLED CAVITY ARRAY
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 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 couplinĝ where J is a hopping parameter, ω c is a resonator frequency, andâ x ≡ 1 √ L kâ k e −ikx is a photonic annihilation operator of the resonator mode at site x. The corresponding dispersion is with wavevector k ∈ [−π, π). 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 Coulombgauge Hamiltonian (1) is given bŷ which corresponds to electromagnetic amplitudes f k (see Eq. (3)) As the amplitudes are independent of k, it is useful to define the characteristic light-matter coupling strength by Note that this expression is consistent with the definition (18).
We model a quantum emitter as a charged particle trapped in the standard double-well potential, 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. are shown in the unit ω c = = m = 1 throughout this paper. These results are fully consistent with the scaling analysis summarized in Table I. 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 = 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 ∝ g −3 as consistent with the scalings Ω 0 ∝ g and ξ 0 ∝ g −1/2 , while in the Coulomb gauge the photon number increases as ∝ g.
A. 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 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 )/ ω c at g = 0. The eigenstates insensitive to g and staying in the photonic band, 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, 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 δE BS ∝ g −1 and are asymptotically determined by the excitation energies δE α 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, they behave either as the 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 B) 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 suppression. 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/ω c ∼ 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/ω c > 5. Meanwhile, the apparent absence of two-photon scattering states in this regime is motivated by the clarity of presentation, because of which only a certain number of the lowest eigenstates are presented in Fig. 3; this avoids excessive overlaps of the continuous spectra.

B. 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 where we neglect an irrelevant constant, introduce the renormalized potential barrier v eff as and define the effective dipole length by d eff /d ≡ (v eff /v) 1/4 . As expected from the general argument in Sec. III C, the barrier v eff is always suppressed compared to the bare value v and the suppression is most significant when ξ becomes maximum, which occurs around the DSC regime (see Fig. 4(a)). Interestingly, when the displacement parameter ξ becomes sufficiently large such that ξ > d/ √ 3, even the full suppression of the potential barrier, i.e., v eff = 0, is possible. Nevertheless, this does not mean that the entire potential is suppressed because the dipole length d eff in Eq. (51) also converges to zero in this case. The resulting potential then contains both the quartic and quadratic contributions with the same sign, leading to the single minimum (see Fig. 4(b) for an illustrative example).

C. 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, where the renormalized emitter frequency ∆ g and the effective coupling strengthsg n are defined by We recall that E 1,2 (|ψ 1,2 ) represent the two-lowest eigenenergies (eigenstates) of the renormalized emitter Hamiltonian (22), and thus depend on the coupling strength g through m eff and V eff . Importantly, since the effective spin-bath couplingsg n remain small over the entire region (cf. Fig. 2(b)), the rotating wave approximation in Eq. (53) 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 D. We thus expect the effective modelĤ JC U to be valid not only at weak g, but also in the nonperturbative regimes.
To demonstrate this, we plot in Fig. 5(a) the lowest excitation energy E ex ofĤ JC U which is obtained from the following analytic relation for the single-excitation subspace, 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 bŷ

TBIC
TqBIC-BS 2 Tosc  61)). 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 Tosc = 2π/ωosc that diverges in the strong-coupling limit (cf. Eq. (63) where ∆ andg k are the bare parameters defined by with x ωc = /mω c . This simplified HamiltonianĤ JC C takes exactly the same form asĤ JC U in Eq. (53), 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(a)).
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Ĥ JC U constructed in the AD frame should also ultimately become invalid when ∆ g <g n ω c , where the counter rotating terms turn out to be equally important as rotating ones; this typically occurs in g 5.
Instead, a more accurate description including the counter rotating terms still remains valid even in such ESC regime. Specifically, we can construct the quantum Rabi model in the AD frame, which gives almost the exact results in the ESC regime (see, for instance, the comparison in Fig. 5(b)). There, we note that the lowest excitation energy exponentially vanishes as g is increased (cf. Eq. (93) below), while the higher excitation energies lie well above this two-level manifold with the energy spacing that scales as ∝ 1/g. This is the reason why the quantum Rabi description becomes asymptotically exact in the AD frame.

D. 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 → d f at time t = 0 while keeping all the other parameters constant. This effectively changes the positions of the minima of V eff and also modifies the qubit frequency. The initial state |Ψ(0) 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Ĥ U at d = d f (see Appendix B for details) and obtain the time evolution via where E i (|Ψ i U ) are the corresponding energies (eigenstates), and c i are expansion coefficients of the initial state. We then calculate the evolution of an observableÔ in the original Coulomb gauge through the unitary transformation Figure 6(a,b) shows the typical spatiotemporal dynamics of photon occupancy n x = â † xâx 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 [109,110].
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 = 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 eff . Because the mass is enhanced as m eff ∝ g 2 and the wavepacket is tightly localized, this quench initiates the oscillatory dynamics where the wavepacket (initially localized at d d i ) oscillates around the new minima at d d f . Such oscillation frequency can be estimated as which vanishes in the strong-coupling limit. The estimated period T osc = 2π/ω osc agrees well with the observed longperiod 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.

A. 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, . . . , N . We thus start from the multi-emitter QED Hamiltonian in the Coulomb gauge, where the position and momentum operators of the emitters satisfy and the vector potential is given bŷ 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 multiemitter case, we obtain the following Hamiltonian (see Appendix C for details): where the emitter part is given bŷ Here, G is the N ×N matrix defined by Physically, m eff,j in Eq. (69) represents the renormalized mass similar to m eff in the single-emitter case considered before except for its dependence on emitter positions through G. The coupling µ ij in Eq. (70) 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 (68), the vacuum-dressed effective potentials V eff,j are given by The interaction Hamiltonian iŝ where the displacement parameters ξ nj characterize the effective coupling strengths between dressed photon mode n and emitter j (see Appendix C). The photon partĤ light takes the same form as the single-emitter case in Eq. (27) [111].
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 ∀j. This simplification leads to the identical bare light-matter couplings, g kj = g k ∀j. In contrast, we note that the effective masses m eff,j , the displacement parameters ξ nj , and the dressed potentials V eff,j can still depend on the emitter positions, and thus we need subscript j to distinguish them in general. Below we illustrate aspects of several limiting cases, but leave the full understanding of multi-emitter waveguide QED systems to a future work.

B. Two emitters
We begin by discussing the two-emitter case. Figure 7 plots the waveguide-mediated coupling µ 21 and the effective mass m 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 µ 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 1+4Θ 1+2Θ at zero separation and eventually saturates to the single-emitter limit m(1+2Θ) when the separation surpasses the cavity length x ωc = 1 (cf. Eq. (15) and (16)).

C. Localized N emitters
We next consider the case in which all the emitters are coupled to the waveguide at the same position x 1 = · · · = 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 intro- duce the collective momentum and coordinate bŷ as well as the relative ones viâ The Hamiltonian can be rewritten aŝ whereĤ governs the dynamics of the collective mode with the effective mass The electromagnetic interaction between the collective mode and photons is given bŷ The relative motion of emitters is governed bŷ Here we recall that the relative variables satisfy the constraints jp j = jq j = 0 and thus they contain N − 1 degrees of freedom. When the collective mode dominantly couples to the electromagnetic fields, one may neglect fluctuations and dynamics of the relative degrees of freedom. The total QED Hamiltonian (79) then becomes equivalent to the singleemitter one (14) upon the replacementsQ →Q CM ,P →P CM , g k → √ N g k , d → √ N d, and v → N v. While this equivalence is nothing but the well-known √ N collective enhancement of the dipole and the coupling strength, our analysis indicates that it can remain even in the DSC/ESC regimes where the multilevel nature of emitters becomes crucial, as long as relative motion does not play a substantial role.

D. Two-level effective model
We next extend the construction of the two-level effective model discussed in Sec. IV C to arbitrary multi-emitter systems. The projection onto the two-lowest dressed emitter states in the AD frame results in the effective model where the matter part corresponds to the (inhomogeneous) transverse-field Ising model, with J ij being the waveguide-mediated qubit-qubit interaction given by Here, |ψ 1,2j represent the two-lowest eigenstates of the renormalized emitter HamiltonianP 2 j /2m eff,j +V eff,j , and ∆ g,j ≥ 0 is the corresponding excitation energy. We recall that ∆ g,j depends on emitter positions through m eff,j and V eff,j . The Jaynes-Cummings-type light-matter interaction iŝ whereg nj are the effective qubit-boson couplings given bỹ As discussed before, in contrast to the Coulomb/PZW gauges, our construction in the AD frame should remain valid even at large g because the level truncations can increasingly be well-justified in the strong-coupling limit. Thus, the twolevel effective model (84) can be used to accurately capture low-energy physics of the original multi-emitter QED Hamiltonian in a broad range of the coupling strength. Nevertheless, we note that in the ESC regime the counterrotating terms can be important and, in such a case, the Rabi-type interaction instead of Eq. (87) should give more accurate results. 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 may set the photon state to be the vacuum and reduce the whole problem to the Ising Hamiltonian (85), from which the ground-state properties and elementary excitations can be extracted. 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, ∆ g,j = ∆ g ,g nj =g n ∀j, and the spin-spin interaction becomes translationally invariant J ij = J |i−j| . Then, the multi-emitter Hamiltonian (85) 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: where J > 0 is the all-to-all antiferromagnetic coupling,Ŝ γ ≡ jσ γ j are the collective spin operators with γ ∈ {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 multiemitter QED Hamiltonian in the nonperturbative regimes. Indeed, it is known that, in the many-emitter limit N → ∞, such model exhibits rich phase diagrams depending on dimension, lattice geometry, or sign/decay exponent of the longrange coupling J |i−j| [112][113][114][115]. Moreover, a disordered transverse-field Ising model is argued to realize many-body localization [116,117], and such disorder is fairly ubiquitous in the multi-emitter Hamiltonian (85) where disorder comes into play through emitter positions. While we leave the full understanding of such multi-emitter physics at strong lightmatter 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, in the low-energy limit, scale as where l > 0 is 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. 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. Below we discuss that the presence or absence of transition can be understood in terms of the mass renormalization after the unitary transformation, and demonstrate it by analyzing a concrete model of circuit QED.

A. Delocalization-localization transition and mass renormalization
The ground state of a single-emitter system displays either delocalized or localized phase that is characterized by the following order parameter where Q h,C represents an emitter displacement in the ground state of the QED Hamiltonian in the Coulomb gauge (1) with a bias potential −hQ being added to V (Q); from now on, we 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 Q +Ξ U (cf. Eq. (62)). The delocalized phase is characterized by the vanishing order parameter O = 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 → ∞ corresponding to localization to each of the two minima in the potential. In this case, the order parameter takes a nonzero value O > 0, which indicates the broken Z 2 symmetry in Eq. (34).
At sufficiently large coupling g, the emitter and photons are decoupled in the AD frame, and the first excitation energy ∆ g can be estimated from the tunneling rate in the dressed potential (23), resulting in [118] Thus, the divergent m eff leads to the gap closing lim L→∞ ∆ g = 0, indicating a possible ground-state degeneracy, i.e., transition to the localized phase. In contrast, as long as m eff remains finite, such exact two-fold degeneracy of the ground state is unlikely to happen; this implies the absence of transition. We delineate general properties in each of these cases on the basis of this observation. Firstly, when the effective mass remains finite m eff < ∞ 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 wellapproximated by the lowest decoupled eigenstate |ψ 0 |0 (cf. Eq. (28)), which has O = 0 and provides the unique ground state due to the nonvanishing excitation energy ∆ g > 0 (see Eq. (93)). Note that, while ∆ g can be exponentially small as g is increased, it still remains nonzero in L → ∞ at any finite g. Thus, the ground state is not expected to exhibit the exact two-fold degeneracy and should remain delocalized. It is worthwhile to note that the same argument should also rule out the possibility of a phase transition for general gapped dispersions, which include some experimentally relevant situations, such as cavity array and (finite-size) open transmission lines.
Secondly, in certain gapless dispersions, the effective mass m eff exhibits the infrared divergence and grows polynomially as a function of system size L. One can also check that this leads to the polynomial divergence of ξ in the dressed emitter potential (23); the resulting effective potential V eff then has the unique minimum at Q = 0 (cf. Eq. (52)) and becomes infinitely tight in L → ∞. 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Ĥ An order of its ground state is characterized by the expectation value of Ξ U (see Eq. (92)). In the limit of a deep potential with large v, the first term in Eq. (94) is dominant, and the ground state should exhibit the two-fold degeneracy at Ξ U = ±d in L → ∞. This leads to the localized phase with O > 0. In contrast, in the opposite limit of a shallow potential, the second (harmonic) term in Eq. (94) dominates over the potential term, and the vacuum state gives the unique ground state, which has Ξ U = 0 and leads to the delocalized phase with O = 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 O 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 ∆.
To sum up, one expects a continuous quantum phase transition between the localized phase at strong v (resp. small ∆) and the delocalized phase at weak v (resp. large ∆); see the blue dashed vertical arrow in Fig. 8. Finally, the present consideration of the full QED Hamiltonian in nonperturbative regimes also reveals an intriguing new possibility beyond what is commonly expected before. Namely, as inferred from the nonmonotonic g dependence of the coupling coefficients ζ n in Eq. (4) and accordingly ξ 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 all the other energy scales. Physically, the origin of this favoring of the delocalized phase can be traced back to the diamagneticÂ 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 [119]; this is reminiscent of what has been discussed in the context of the no-go theorems of the superradiant transition [68][69][70][71][72][73].

B. Case study of the circuit QED Hamiltonian
One may expect that the present QED Hamiltonian with the double-well potential should feature the similar physics as known for the spin-boson model. For instance, it is often supposed that the two-level projection of the PZW Hamiltonian (40) should allow for the spin-boson description of the waveguide QED systems. However, as discussed before, such level-truncation procedure cannot in general be justified at strong couplings, and we must carefully reexamine the ground-state properties of nonperturbative waveguide QED systems separately from the simplified spin-boson description. Indeed, in the previous section, we point out that the fullfledged QED systems should exhibit a new feature, which is not present in the usual spin-boson models [120], such as the reentrant transition into the delocalized phase at sufficiently large light-matter coupling.
In this section, we concretely demonstrate these results in the case of resistively shunted Josephson junctions by using the functional renormalization group (FRG) analysis. For the sake of convenience, we here switch to the notation familiar with the circuit QED community. Specifically, we consider a microscopic circuit Hamiltonian (we set = 1) [121]: where ϕ is the Josephson junction phase andN = −i∂/∂ϕ is the charge operator, which play roles as the positionQ and momentumP operators in the previous notation, respectively. The form of the Hamiltonian (95) precisely corresponds to the Coulomb-gauge-type Hamiltonian in Eq. (4) obtained after the Bogoliubov transformation. The charging energy is denoted by E C and the potential term V (ϕ) is chosen to be the double-well potential in the same way as before: In practice, such potential is routinely realized in flux qubits by combining the inductive energy and flux-tuned Josephson energy. The coupling coefficient ζ m and the environmental frequency Ω m of mode m ∈ {1, 2, . . . , M − 1} are given by where W is the environmental cutoff frequency and α is the dimensionless parameter characterizing the coupling strength; the latter is related to the shunt resistance R via α = R Q /R with R Q = h/(4e 2 ) being the quantum of resistance. The characteristic coupling strength g and environmental frequency ω defined in Eqs. (18) and (17) basically correspond to √ αE C W and W in the present notation, respectively. Thus, for instance, the ESC regime corresponds to the region αE C W . We note that, when the environmental cutoff W satisfies the weak coupling condition αE C W , the lowenergy limit (m M ) of Eq. (95) reproduces the Caldeira-Leggett description of the Ohmic dissipation [18,121,122]. Since we are here interested in the ESC regime αE C W , we need to analyze the original circuit QED Hamiltonian (95) without taking such limit.
Before providing a quantitative analysis, we illustrate the general features by transforming to the AD frame: where we define the dimensionless coupling strength (which basically corresponds to g/ω in the previous notation) by and useĤ U =Û †Ĥ CÛ withÛ = exp(−iNΞ) andΞ = m i(b † m −b m )ζ m /Ω m . Importantly, the effective "mass" in the transformed frame shows the infrared divergence in η ≤ 1, while it remains finite in the ESC regime η > 1. Thus, following our arguments above, we expect that as the coupling strength is increased there occurs the delocalizationlocalization transition in η ≤ 1, while the ground state should ultimately exhibit the reentrant transition to the delocalized phase above η ∼ 1.
To make these predictions concrete, in Fig. 9 we show the ground-state phase diagram of the circuit QED system (95), which is obtained by the FRG analysis with the local potential approximation (see, e.g., Refs. [121,123,124] for technical details). In the limit of deep potential barrier 1/v → 0, ϕ 0 always remains to be positive during the RG flows and thus the localized phase is realized. As the potential barrier becomes shallow (1/v increases), the value of ϕ 0 is eventually renormalized to zero in the IR limit and the transition to the delocalized phase occurs. Notably, the behavior of the transition point drastically changes around η = 1 (vertical dashed line); in η < 1, the localized phase expands as the coupling strength η is increased, while the ordering is suppressed in the ESC regime η > 1. This is consistent with our arguments based on the AD frame above and also with the nonmonotonic dependence of ζ m on the coupling strength α in Eq. (97). We expect that these results might be tested by recent experiments realizing galvanic coupling of Josephson junctions to a high- impedance long transmission line [125,126].

VII. SUMMARY AND DISCUSSIONS
We analyzed equilibrium and dynamical properties of lightmatter 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. Below we summarize our key findings.
We first analyzed 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 [51,79] and plasmonic crystals [36] or (at least) expected to be accessible in the near future in view of rapid developments in achieving stronger light-matter coupling regimes [60,62]. To explore those nonperturbative regimes in setups of atoms coupled to photonic crystals, one can utilize Rydberg atoms [127][128][129] and/or the collective √ N enhancement of the light-matter coupling by assembling a large number of weakly coupled components [130] as discussed in Sec. V C. 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 argued that the ground-state physics can be understood from the perspective of the transverse-field Ising model (85) but with suitably renormalized parameters. Finally, we analyzed the case of gapless photonic dispersions and showed 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. There, we found that in certain cases the diamagnetic term leads to the suppression of the symmetry-broken phase, which is reminiscent of no-go theorems for superradiant transition in cavity QED (see e.g., Refs. [131,132]). One surprising consequence of this analysis is the appearance of a reenterant transition into the symmetryunbroken phase at sufficiently strong coupling, which was absent in the simplified descriptions such as the spin-boson model. These results are also confirmed by the FRG analysis of the circuit QED system.
It is interesting to analyze the ground-state properties of multi-dipole waveguide QED systems in further detail. In particular, the full understanding of a possible superradiant-type transition in the present multi-dipole systems remains as an intriguing open question. At a qualitative level, our analysis appears to indicate that the tendency to superradiance (i.e., localized phase) will be the strongest at intermediate coupling strengths. For instance, as far as the collective mode plays a dominant role and relative motion can be neglected, we expect that the analysis in Sec. VI can be extended to the multi-dipole cases via replacing the effective mass m eff by the collective one M eff (cf. Eq. (80)). If the bosonic dispersion is gapless and M eff exhibits the infrared divergence, then the ground state may exhibit the ordering akin to the superradiant phase, in a close analogy with the delocalization-localization transition for the single-dipole cases discussed in this paper. One of the conceptual advantages in our approach in this respect is to connect these seemingly unrelated phenomena.
Taking the limit of many emitters should provide alternative way to realize a quantum phase transition. In particular, one may use the two-level effective model (85) to determine the ground-state properties in such cases. This mapping to the transverse-field Ising model suggests 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.
Several further open questions also 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. They essentially broaden the absorption spectra shown in Fig. 6(g-i) and can affect the dynamics as excitations acquire finite lifetimes. While the waveguide coupling efficiency can be made close to the unit fidelity in, e.g., microwave superconducting devices [9], those loss effects are still ubiquitous in atoms coupled to photons in optical domains [133]. These issues can be addressed by combining the present nonperturbative QED formalism with the standard framework of Markovian open systems [134]. Second, instead of the exact diagonalization performed here, one can apply more efficient numerical methods, such as the matrix-productstates calculations [45,47,52,85] or a hybrid variational approach [135], 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. [136]) whose inclusion is often important depending on the cavity geometry and the coupling strength. We also note that the present formalism can be extended to higherdimensional systems (see, e.g., the Supplementary Materials of Ref. [65]). One intriguing direction is to explore a possible extension of these formalisms to the case in which matter degrees of freedom consist of indistinguishable quantum particles.
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. In particular, in the case of localized N identical emitters, our results obtained for a single-emitter system are expected to remain the same (except for the √ N enhancement) as far as the relative motion does not play a significant role (see Sec. V C). In this respect, the predicted vacuum-induced suppression of the potential barrier in Sec. III C may lie at the heart of the enhanced chemical reactivity observed in polaritonic chemistry [28][29][30]. More generally, our study reveals that the mass enhancement is one of the universal features of strongly interacting lightmatter systems. This naturally associates with the higher density of states, which could lead to enhancements of certain many-body properties, including superconductivity or ferromagnetism [137][138][139][140][141][142][143][144][145]. and rewrite the quadratic part ofĤ C in Eq. (1) as The last term in Eq. (A3) can readily be diagonalized by an orthogonal transformation, where an orthogonal matrix O satisfies Eq. (7). The quadratic photon Hamiltonian then becomes where we introduce the annihilation operator after the orthogonal transformation bŷ b n = Ω n 2 which gives Eq. (5). Using these squeezed photon operatorŝ b n , the total Hamiltonian in the Coulomb gauge is now given by Eq. (4) in the main text.

Appendix B: Numerical diagonalization
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, a e x=0 =â x=0 , a e x>0 = 1 interact with the emitter, we neglect contributions from the odd modes. The photon Hamiltonian in the real-space basis is then written aŝ The vector field is expressed in this basis aŝ which results in the electromagnetic amplitudes (cf. Eq. (3)) We use the frequencies ω p and the amplitudes f p to diagonalize the photon Hamiltonian including theÂ 2 term as explained in Appendix A. Then, we perform the unitary transformation (9) and arrive at the Hamiltonian in the AD framê To numerically diagonalize this Hamiltonian efficiently, we employ the few-photon ansatz (31) and express the matrix elements ofĤ U in terms of the basis |ψ α emitter ⊗ |n 0 n 1 · · · n (L−1)/2 photon (B7) with level truncations, α = 1, 2, . . . , α c , We recall that |ψ α are single-particle eigenstates of the renormalized emitter Hamiltonian (22), and |n 0 n 1 · · · n (L−1)/2 is a many-body bosonic Fock state in terms ofb n operators with n j = 0, 1, · · · . The corresponding Hilbert-space dimension is which grows polynomially with the system size L. Figure 10 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 α c = O(10) to achieve the accuracy with an error below ∼ 1%. When only the low-energy spectrum is of interest, one can use the Lanczos method to further reduce the computational cost. term. To do so, we introduce the position-dependent multiemitter coupling strengths by and rewrite the quadratic part as (aside constant) where we recall that the conjugate operatorsX k andP k are defined by Eqs.
where we define the squeezed photon operatorsb n in the same manner as in Eq. (A6). In terms of these new photon operators, the Coulomb-gauge Hamiltonian is expressed aŝ where we define ζ nj = m j Ω n × k g c kj S XX kn − g s kj S P X kn + i g c kj S XP kn − g s kj S P P kn .

(C8)
We now introduce the multi-emitter extension of the asymptotically decoupling unitary transformation bŷ where we recallΞ j = n i(ξ njb † n −ξ * njb n ) with the displacement parameters This transformation acts on the photon and emitter operators asÛ †b and transforms Eq. (C7) to the following form: which gives Eq. (67) in the main text. Here, we introduce the effective mass m eff,j for each emitter and the emitter-emitter coupling µ ij as The expressions (69) and (70) in Sec. V follow from the relations (C6) and (C8).