# Bound states in microwave QED: Crossover from waveguide to cavity regime

N. Pradeep Kumar  ARC Centre of Excellence for Engineered Quantum Systems,
School of Mathematics and Physics, The University of Queensland, Saint Lucia,
Queensland 4072, Australia  Andrés Rosario Hamann  ARC Centre of Excellence
for Engineered Quantum Systems, School of Mathematics and Physics, The
University of Queensland, Saint Lucia, Queensland 4072, Australia  Department
of Physics, ETH Zürich, CH-8093 Zürich, Switzerland  Rohit Navarathna  ARC
Centre of Excellence for Engineered Quantum Systems, School of Mathematics and
Physics, The University of Queensland, Saint Lucia, Queensland 4072, Australia
Maximilian Zanner  Center for Quantum Physics, and Institute for Experimental
Physics,University of Innsbruck, A-6020 Innsbruck, Austria  Institute for
Quantum Optics and Quantum Information of the Austrian Academy of Sciences,
A-6020 Innsbruck, Austria  Mikhail Pletyukhov  Institute for Theory of
Statistical Physics, RWTH Aachen University, 52056 Aachen, Germany  Arkady
Fedorov<EMAIL_ADDRESS>ARC Centre of Excellence for Engineered Quantum
Systems, School of Mathematics and Physics, The University of Queensland,
Saint Lucia, Queensland 4072, Australia

(August 27, 2024)

###### Abstract

Light-matter interaction at the single-quantum level is the heart of many
regimes of high fundamental importance to modern quantum technologies. Strong
interaction of a qubit with a single photon of an electromagnetic field mode
is described by the cavity/circuit electrodynamics (QED) regime which is one
of the most advanced platforms for quantum computing. The opposite regime of
the waveguide QED, where qubits interact with a continuum of modes in an
infinite one-dimensional space, is also at the focus of recent research
revealing novel quantum phenomena. Despite the demonstration of several key
features of waveguide QED, the transition from an experimentally realizable
finite-size system to the theoretically assumed infinite device size is
neither rigorously justified nor fully understood. In this paper, we formulate
a unifying theory which under a minimal set of standard approximations
accounts for physical boundaries of a system in all parameter domains.
Considering two qubits in a rectangular waveguide which naturally exhibits a
low frequency cutoff we are able to account for infinite number of modes and
obtain an accurate description of the waveguide transmission, a life-time of a
qubit-photon bound state and the exchange interaction between two qubit-photon
bounds states. For verification, we compare our theory to experimental data
obtained for two superconducting qubits in a rectangular waveguide
demonstrating how the infinite size limit of waveguide QED emerges in a
finite-size system. Our theory can be straightforwardly extended to other
waveguides such as the photonic crystal and coupled cavity arrays.

##  I Introduction

Engineering interaction between an atom or other quantum emitters with
quantized electromagnetic fields serve as the building block of modern quantum
technologies. In particular, cavity quantum electrodynamics (cQED) that
describes the strong coupling between an atom and a single mode of a high Q
cavity is now a well understood theory and underlies the development of
quantum computers based on superconducting qubits [1, 2]. The natural
extension of cQED leads to multimode-cQED, where several discreet modes of the
cavity are simultaneously hybridized with an atom [3, 4]. The multimode cavity
satisfies the condition that $\kappa\ll\delta\omega$, where $\kappa$ is the
linewidth of the modes and $\delta\omega$ is the free spectral range that
measures the mode spacing [5]. In the opposite case of $\kappa\gg\delta\omega$
results in a continuum of overlapping modes that can be utilized to study
waveguide-QED phenomenon that emulates the dynamics of an atom embedded in a
free space [6, 7]. In contrast to the standing waves of a cavity, waveguides
support a large bandwidth of propagating waves that enables long-range photon
mediated interaction between distant quantum emitters [8, 9]. These properties
find direct applications in quantum communications for realizing quantum
networks [10, 11], non-reciprocal photon transmission [12, 13] and routing
[14].

Waveguides can also be engineered to support non-trivial dispersion relation
that can result in photonic bandgaps or stopbands. When the qubit frequency
lies in the vicinity of the band edge, a pair of dressed states are predicted
to emerge. One of these dressed state merges into the continuum while the
other dressed state may be pushed into the bandgap wherein the photonic part
of the dressed state becomes an evanescent wave exponentially localized around
the qubit [15, 16, 17, 18]. Such a dressed state is referred to as the qubit-
photon bound state (QPBS) and has been observed in the experiments with
photonic crystals [19, 20] and metametarials [21, 22]. It has been proposed
that the interaction between an array of such bound states can be engineered
for simulating tunable spin-exchange interaction and quantum impurity models
[23, 24, 25, 26, 22].

A common approach to generate stopbands is to use photonic crystals and sub-
wavelength [19, 20] meta-materials [27, 21] built using an array of coupled
cavities, which are typically fabricated as 2D co-planar structure on a chip
that provides transverse confinement to the travelling photons. As an
alternative, this can also be achieved in 3D waveguides [28, 29] which has
naturally occurring low-frequency cutoff due to its geometry. Furthermore
coupling between the qubits and the waveguide can be enhanced by tailoring the
electric dipole of the superconducting qubit in the 3D Transmon architecture
[30].

Experimental results are usually explained by postulating infinite size of a
device and neglecting coupling to the input and output ports. However,
physical systems employed in experiments are finite dimensional and interacts
with measurement apparatus. For the passband the infinite nature is
experimentally justified by using arguments of impedance matching at the
coupling ports. This argument is especially questionable for a 3D waveguide
whose configuration is conceptually indistinguishable from a 3D cavity and
signature of the standing waves always remain. Similarly, photonic crystals
and coupled-cavity array implementations use a finite number of elementary
cells and the finite nature of the system is even more pronounced. Thus, a
question of validity of the infinite limit as opposed to the description which
utilise the single-mode (or multi-mode) Jaynes-Cummings model [5] remains
largely an open question.

In this work we developed a unifying theory of qubits interacting with
multiple number of electromagnetic modes coupled to the ports. We treat both
the ports and the qubits in a similar footing, by considering the dipole type
coupling with the waveguide. Furthermore, we apply the Markov approximation
(only for the ports), the rotating wave approximation and considering only
single-photon excitation subspace, we apply the Green’s function formalism
[31, 32] which allows one to obtain analytical results for all quantities of
interest. We then apply this formalism to a specific case of a rectangular 3D
waveguide with two qubits where we are able to sum over infinite number of
waveguide modes. We evaluate transmission through the waveguide (both empty
and with the qubits), lifetime of qubit-photon bound states in the stopband
and exchange interaction of two qubit-photon bound states. Our theory shows
how both the infinite limit and the opposite limit of the single-mode Jaynes-
Cummings interaction emerges for a spatially finite system in different
parameter domains. A particular outcome of our theory is the generalised
formula for the Purcell limit and the exchange interaction between qubit which
can include any (up to infinity) number of modes.

We compare our theory to experimental results for two superconducting qubits
in a copper rectangular waveguide and demonstrate how the finite-size effects
become negligible and the infinite waveguide limit results in an exponential
scaling for bound state interaction strength. In the final part of the paper
we show how our formalism can be easily extended to any arbitrary waveguide
with a known dispersion relation. In particular, we discuss the case of
photonic crystal where summation over infinite number of mode is not even
necessary.

##  II General theory of qubit in a finite size waveguide

Let us model a generic waveguide taking into account its finite size, the
length of the waveguide along the $z$ direction is denoted by $L$. The
Hamiltonian of the free waveguide can be written as,

$H_{0}=\sum_{k}\varepsilon_{k}a_{k}^{\dagger}a_{k}\quad k=\frac{l\pi}{L},\quad
l=1,2,\dots,$ (1)

where $k$ is the wavevector along the $z$ direction of the waveguide and
$\varepsilon_{k}$ is the dispersion relation associated with the waveguide,
and $a_{k}(a_{k}^{\dagger})$ is the annihilation (creation) operator of the
corresponding mode. Note that both the boundary conditions and the dispersion
relation can be kept arbitrary at this point.

As shown in Fig. (1), the waveguide is coupled to the (left and right) ports
at the points $z_{L}=-\frac{L}{2}+d_{e}$ and $z_{R}=+\frac{L}{2}-d_{e}$ as
well as to the two qubits at the points $z_{1}$ and $z_{2}$. We describe this
system by the full Hamiltonian

$\displaystyle
H=H_{0}+\sum_{s=L,R}(H_{s}+H_{s,int})+\sum_{j=1,2}(H_{j}+H_{j,int}).$ (2)

Here we keep the number of qubits to be two for brevity as this allows us to
calculate all the single-qubit effects and the most fundamental cooperative
effects such as the exchange interaction between two qubits. This
configuration is also reflected in our experiment with a 3D waveguide as
presented below.

Hereby the ports are modelled as structureless Markovian reservoirs with
linear dispersion,

$\displaystyle H_{s}=\int
d\omega_{s}\,\omega_{s}\,c^{\dagger}_{\omega_{s}}\,c_{\omega_{s}},$ (3)

in terms of the corresponding continuum field operators obeying the canonical
commutation relations
$[c_{\omega_{s}},c_{\omega^{\prime}_{s}}^{\dagger}]=\delta(\omega_{s}-\omega^{\prime}_{s})$.
The ports’ coupling to the waveguide is of the dipole type treated in the
rotating wave approximation (RWA)

$\displaystyle
H_{s,int}=\sum_{k}f_{k}a_{k}^{\dagger}\sqrt{L}\psi_{k}(z_{s})\sqrt{\frac{\Gamma_{s}}{\pi}}\int
d\omega_{s}c_{\omega_{s}}+h.c.,$ (4)

where $\psi_{k}(z)$ is the $k$-th eigenfunction of the waveguide Hamiltonian
(1). Here, we assume that all modes of the same port are equally coupled to a
given waveguide’s mode $k$ — this assumption supports the Markovian modelling
of the ports. In turn, the waveguide’s dipole moment may have a dispersion in
the waveguide’s mode index $k$. Typically this dependence is
$\propto\sqrt{\varepsilon_{k}}$, so we define
$f_{k}=\sqrt{\frac{\varepsilon_{k}}{\omega_{c}}}$, where $\omega_{c}$ is some
relevant frequency scale (e.g. the low-energy cutoff frequency, like it
appears in our following application).

The two qubits (labelled by $j=1,2$) with the transition frequencies
$\omega_{q,j}$ are described by the effective Hamiltonians

$\displaystyle
H_{j}=\left(\omega_{q,j}-i\frac{\gamma_{a,j}}{2}\right)\frac{1+\sigma_{z}^{(j)}}{2},$
(5)

which also take into account their nonradiative decay at the rates
$\gamma_{a,j}$. The qubits’ coupling to the waveguide is analogous to (4) —
the dipole coupling in the RWA:

$\displaystyle
H_{j,int}=\sum_{k}f_{k}a_{k}^{\dagger}\sqrt{L}\psi_{k}(z_{j})g_{j}\sigma_{-}^{(j)}+h.c.,$
(6)

where $\sigma_{z}^{(j)}$ and
$\sigma_{-}^{(j)}=\frac{1}{2}[\sigma_{x}^{(j)}-i\sigma_{y}^{(j)}]$ are the
standard Pauli matrices.

We note that the presented model is rather general and can be also applied to
a lattice realization of the waveguide. In this case the sum over $k$ in (1)
is finite, and the coordinate $z_{n}=na$ is discrete and given by multiples of
the lattice constant $a$.

Figure 1: Schematic of a typical waveguide QED set-up employed in experiments.
The waveguide is finite in size with boundaries located at $\pm L/2$. The
waveguide is driven from left ($c_{in}$) to right ($c_{out}$) via input-output
ports located at $z_{L,R}$ with th corresponding coupling strengths denoted by
$\Gamma_{L,R}$. Pair of two-level quantum emitters ($\omega_{q,1,2}$) are
coupled to the waveguide and placed symmetrically at $\pm d/2$ with respect to
the origin and decays into the waveguide at the rate $g_{1,2}$.

###  II.1 Calculating transmission

With the help of the Hamiltonian (2) we derive (see details in Appendix A) the
Heisenberg equations of motion for the (ports’ and waveguide’s) field
operators and the qubits’ operators. Assuming that the incident field from the
left port enters the waveguide at $t=0$, we define the input field

$\displaystyle c_{in}(t)=\int d\omega_{L}c_{\omega_{L}}(t+0^{+})$ (7)

in terms of the left port operators. Since the transmitted field is measured
in the right port, we define the output field

$\displaystyle c_{out}(t)=\int d\omega_{R}c_{\omega_{R}}(t-0^{+})$ (8)

in terms of the right port operators. The infinitesimal shifts of the time
arguments in the above expressions reflect the causality property — an
important feature of the input-output formalism [33]. The transmission
amplitude $S_{RL}(\omega)$ relates the Fourier images of (7) and (8),

$\displaystyle\tilde{c}_{out}(\omega)=S_{RL}(\omega)\tilde{c}_{in}(\omega).$
(9)

In the single-excitation approximation, which is justified for a weak pulse
intensity and which is based on the assumption that both qubits remain close
to their ground states in the time dynamics, that is
$\langle\sigma_{z}^{(j)}(t)\rangle\approx-1$, we establish the following
relation (see again in Appendix A for more details)

$\displaystyle
S_{RL}(\omega)=-2i\sqrt{\Gamma_{R}\Gamma_{L}}LG(z_{R},z_{L};\omega).$ (10)

Here

$\displaystyle
G(z,z^{\prime};\omega)=\sum_{k,k^{\prime}}f_{k}\psi_{k}(z)G_{kk^{\prime}}(\omega)f_{k^{\prime}}\psi_{k^{\prime}}(z^{\prime})$
(11)

is the waveguide Green’s function written in the coordinate representation.
Its momentum-space counterpart

$\displaystyle
G_{kk^{\prime}}(\omega)=\left(\frac{1}{\omega-\hat{\varepsilon}-\Sigma(\omega)}\right)_{kk^{\prime}},\quad\hat{\varepsilon}_{kk^{\prime}}=\varepsilon_{k}\delta_{kk^{\prime}},$
(12)

is dressed by the self-energy

$\displaystyle\Sigma_{kk^{\prime}}(\omega)=\Sigma_{kk^{\prime}}^{(p)}+\Sigma_{kk^{\prime}}^{(q)}(\omega),$
(13)

which takes into account the effects of the interaction with the ports

$\displaystyle\Sigma_{kk^{\prime}}^{(p)}=-iL\sum_{s=L,R}\Gamma_{s}f_{k}\psi_{k}(z_{s})f_{k^{\prime}}\psi_{k^{\prime}}(z_{s}),$
(14)

and with the qubits

$\displaystyle\Sigma_{kk^{\prime}}^{(q)}(\omega)=L\sum_{j=1,2}g_{j}^{2}\frac{f_{k}\psi_{k}(z_{j})f_{k^{\prime}}\psi_{k^{\prime}}(z_{j})}{\omega-\omega_{q,j}+i\frac{\gamma_{a,j}}{2}}.$
(15)

Note that (14) is frequency independent, which reflects the Markovian
modelling of the ports.

An alternative representation for $S_{RL}(\omega)$ can be worked out which
will be applied later in Section III.2 for the analysis of qubit-photon bound
states. In particular, in order to identify these states with resonant peaks
in the transmission data below the waveguide’s cutoff frequency, we
equivalently rewrite $S_{RL}(\omega)$ as

$\displaystyle S_{RL}(\omega)=$
$\displaystyle-2iL\sqrt{\Gamma_{R}\Gamma_{L}}G^{(p)}(z_{R},z_{L};\omega)$ (16)
$\displaystyle-2iL^{2}\sqrt{\Gamma_{R}\Gamma_{L}}\sum_{j,j^{\prime}}g_{j}g_{j^{\prime}}G^{(p)}(z_{R},z_{j};\omega)$
$\displaystyle\qquad\qquad\times\hat{G}_{jj^{\prime}}(\omega)G^{(p)}(z_{j^{\prime}},z_{L};\omega).$
(17)

Here $G^{(p)}(z,z^{\prime};\omega)$ is the waveguide Green’s function dressed
solely by the port-induced self-energy $\Sigma_{kk^{\prime}}^{(p)}$: One
obtains $G^{(p)}(z,z^{\prime};\omega)$ by replacing $G_{kk^{\prime}}\to
G_{kk^{\prime}}^{(p)}$ and $\Sigma_{kk^{\prime}}\to\Sigma_{kk^{\prime}}^{(p)}$
in (12) and (13), respectively. In addition, we introduce the dressed Green’s
function of the qubits,

$\displaystyle\hat{G}_{jj^{\prime}}(\omega)=\left(\frac{1}{[\hat{G}^{(0)}(\omega)]^{-1}-\hat{\Sigma}(\omega)}\right)_{jj^{\prime}},$
(18)

which is written in terms of the qubits’ bare Green’s function

$\displaystyle\hat{G}_{jj^{\prime}}^{(0)}(\omega)=\frac{\delta_{jj^{\prime}}}{\omega-\omega_{q,j}+i\frac{\gamma_{a,j}}{2}},$
(19)

and the qubits’ self-energy

$\displaystyle\hat{\Sigma}_{jj^{\prime}}(\omega)=g_{j}g_{j^{\prime}}LG^{(p)}(z_{j},z_{j^{\prime}};\omega),$
(20)

emerging due to their interaction with the waveguide.

Thus, the qubit-photon bound states are expected to manifest themselves as
poles of (18). We refer to Appendix B for the proof of the equivalence between
(16), (17) and (10).

###  II.2 Calculating the port-dressed Green’s function of the waveguide

Recapitulating the definition of $G_{kk^{\prime}}^{(p)}(\omega)$ in terms of
the Dyson equation

$\displaystyle
G_{kk^{\prime}}^{(p)}(\omega)=G_{kk^{\prime}}^{(0)}(\omega)+\sum_{k^{\prime\prime},k^{\prime\prime\prime}}G^{(0)}_{kk^{\prime\prime}}(\omega)\Sigma_{k^{\prime\prime}k^{\prime\prime\prime}}^{(p)}G_{k^{\prime\prime\prime}k^{\prime}}^{(p)}(\omega),$
(21)

where
$G_{kk^{\prime}}^{(0)}(\omega)=\frac{\delta_{kk^{\prime}}}{\omega-\varepsilon_{k}}$
is the bare Green’s function of the waveguide, we explicitly solve (21) in
Appendix C. The obtained solution allows us to represent (omitting for brevity
the frequency argument)

$\displaystyle G^{(p)}(z,z^{\prime})=G^{(0)}(z,z^{\prime})$
$\displaystyle-\frac{iL}{D}\Gamma_{R}G^{(0)}(z,z_{R})[1+i\Gamma_{L}LG^{(0)}(z_{L},z_{L})]G^{(0)}(z_{R},z^{\prime})$
$\displaystyle-\frac{iL}{D}\Gamma_{L}G^{(0)}(z,z_{L})[1+i\Gamma_{R}LG^{(0)}(z_{R},z_{R})]G^{(0)}(z_{L},z^{\prime})$
$\displaystyle-\frac{L^{2}}{D}\Gamma_{R}\Gamma_{L}G^{(0)}(z,z_{R})G^{(0)}(z_{R},z_{L})G^{(0)}(z_{L},z^{\prime})$
$\displaystyle-\frac{L^{2}}{D}\Gamma_{L}\Gamma_{R}G^{(0)}(z,z_{L})G^{(0)}(z_{L},z_{R})G^{(0)}(z_{R},z^{\prime}),$
(22)

where

$\displaystyle D$
$\displaystyle=[1+i\Gamma_{R}LG^{(0)}(z_{R},z_{R})][1+i\Gamma_{L}LG^{(0)}(z_{L},z_{L})]$
$\displaystyle+L^{2}\Gamma_{R}\Gamma_{L}G^{(0)}(z_{R},z_{L})G^{(0)}(z_{L},z_{R}),$
(23)

and we also defined the coordinate representation of the bare waveguide
Green’s function

$\displaystyle
G^{(0)}(z,z^{\prime})=\sum_{k}f_{k}^{2}\frac{\psi_{k}(z)\psi_{k}(z^{\prime})}{\omega-\varepsilon_{k}}.$
(24)

It is remarkable that for the special choice of spatial arguments $z=z_{R}$,
$z^{\prime}=z_{L}$ the expression (22) gets considerably simplified:

$\displaystyle G^{(p)}(z_{R},z_{L})=\frac{G^{(0)}(z_{R},z_{L})}{D}.$ (25)

For a weak coupling to the ports we can approximate (22) by its leading order
expansion in $\Gamma$’s

$\displaystyle G^{(p)}(z,z^{\prime})\approx G^{(0)}(z,z^{\prime})$
$\displaystyle-iL\Gamma_{R}G^{(0)}(z,z_{R})G^{(0)}(z_{R},z^{\prime})$
$\displaystyle-iL\Gamma_{L}G^{(0)}(z,z_{L})G^{(0)}(z_{L},z^{\prime}),$ (26)

which already includes the imaginary part necessary for a description of the
qubit-photon bound state resonance broadening. This becomes explicit after
inserting (26) into (20).

##  III Qubits in rectangular waveguide

Up to this point we have nowhere used the waveguide’s specific properties, in
particular its eigenspectrum $\varepsilon_{k}$ and its eigenmodes
$\psi_{k}(z)$. Instead we have expressed all quantities which are necessary
for the transmission calculation in terms of the bare Green’s function of the
empty waveguide (24) written in the coordinate representation. In this section
we demonstrate the application of our approach to model a rectangular
waveguide shown in Fig. (2). It will be a subject of our subsequent
experimental study in the regard of the qubit-photon bound states that arise
in the stopband.

Figure 2: (a) Schematic of a rectangular waveguide containing two 3D Transmon
qubits $Q_{1}$ and $Q_{2}$ separated by a distance $d$. When the qubits’
frequencies lie in the stopband, exponentially localized photonic bound state
can be induced, which are centered at the qubits’ positions. The bound states
can be probed via input-output ports separated by a distance $z_{R}-z_{L}$.
(b) The actual device used in the experiment consisting of copper rectangular
waveguide with co-axial connectors used for input/output ports. (c) Photon
dispersion curve for the dominant $\text{TE}_{10}$ mode. The cut-off frequency
$\omega_{c}/2\pi=6.5213$ GHz and the qubit frequency
$\omega_{q}\,(=\omega_{q,1}=\omega_{q,2})$ that lies in the stopband.

###  III.1 Calculating bare Green’s function of the waveguide

The rectangular waveguide supports two types of eigenmodes, the Transverse
Electric ($\text{TE}_{m,n}$) modes and Transverse Magnetic ($\text{TM}_{m,n}$)
modes, with the fundamental mode being $\text{TE}_{1,0}$. For the rest of this
section we assume both the ports and the qubits are predominantly coupled to
the fundamental mode although coupling to $\text{TM}_{m,n}$ modes can be in
principle taken into account as well. The photon dispersion associated with
the $\text{TE}_{1,0}$ mode is given by [28]

$\varepsilon_{k}=\sqrt{c^{2}k^{2}+\omega_{c}^{2}}.$ (27)

The eigenfunction of the rectangular metal waveguide are given by the
following normal modes,

$\psi_{k}(z)=\sqrt{\frac{2}{L}}\sin\,[k\,(z+L/2)],$ (28)

satisfying the open boundary conditions $\psi_{k}(-L/2)=\psi_{k}(+L/2)=0$,
which are equivalent to the vanishing tangential components of the electric
field ($E_{x}=E_{y}=0$) and the vanishing normal component of the magnetic
field ($B_{z}=0$) at the endpoints of the waveguide.

Let us evaluate the bare Green’s function of the empty waveguide as defined in
Eq. (24). For the continuum model of our present interest, the sum in (24)
includes infinitely many terms. A natural idea to truncate it at some large
longitudinal mode index $N_{r}$ has however a flaw: The summation convergence
with respect to increasing $N_{r}$ is very slow for the frequencies lying in
the stopband. This drawback does not allow us to accurately quantify the
transmission through the empty waveguide at $\omega<\omega_{c}$, as will be
shortly illustrated.

Instead of truncating the sum in (24), we first rewrite it as

$\displaystyle
G^{(0)}(z,z^{\prime};\omega)=\frac{\omega}{\omega_{c}}\sum_{k}\frac{\psi_{k}(z)\psi_{k}(z^{\prime})}{\omega-\varepsilon_{k}}-\frac{1}{\omega_{c}}\delta(z-z^{\prime}),$
(29)

where the delta function emerges due the basis completeness relation
$\sum_{k}\psi_{k}(z)\psi_{k}(z^{\prime})=\delta(z-z^{\prime})$.

The second contribution to (29) is only efficient for $z\approx z^{\prime}$,
and it should be regularized on the length scale of the contact size (which
has been so far treated as point-like). However, due to its frequency
independence the second term can only contribute to inessential energy level
renormalization effects. Therefore, we completely neglect it.

As we are interested in the transmission spectrum at frequencies
$\omega\approx\omega_{c}$, we can also approximate the first term contributing
to (30) and thus obtain

$\displaystyle
G^{(0)}(z,z^{\prime};\omega)=\sum_{k}\frac{\psi_{k}(z)\psi_{k}(z^{\prime})}{\omega-\varepsilon_{k}}.$
(30)

Note that the both approximations made above can be alternatively deduced from
the assumption that $f_{k}=1$, i.e. the waveguide’s dipole moment does not
depend on the mode index $k$.

To calculate the infinite sum in (30), we exploit the method frequently used
in the condensed matter physics which is known as the Matsubara frequency
summation [34]. By the analogy with the bosonic Matsubara frequencies
$\omega_{l}=\frac{2\pi l}{\beta}$ we define a formal analogue
$\beta=\frac{2L}{c}$ of the inverse temperature in the original method. Then
we represent

$\displaystyle
c\,G^{(0)}(z,z^{\prime};\omega)=\frac{1}{\beta}\sum_{i\omega_{l}}g(i\omega_{l}),$
(31)

where in the conventionally denoted summation $\sum_{i\omega_{l}}$ the range
of all integer $l$ from $-\infty$ to $+\infty$ is implied, and $g(w)$ is a
newly introduced function

$\displaystyle
g(w)=\frac{e^{-w\frac{|z-z^{\prime}|}{c}}-e^{-w\frac{|z+z^{\prime}+L|}{c}}}{\omega-\sqrt{-w^{2}+\omega_{c}^{2}}},$
(32)

which can be analytically continued in the complex plane of $w$. It has the
two branch cuts $(-\infty,-\omega_{c})$ and $(\omega_{c},+\infty)$ along the
real axis of $w$, as well as the two poles
$w_{\pm}=\pm\sqrt{\omega_{c}^{2}-\omega^{2}}\equiv\pm\bar{w}$.

The main prescription of the (bosonic) Matsubara frequency summation consists
in replacing the sum by the complex contour integral

$\displaystyle\frac{1}{\beta}\sum_{i\omega_{l}}g(i\omega_{l})=\frac{1}{2\pi
i}\oint dw\frac{g(w)}{1-e^{-\beta w}},$ (33)

where the integration contour embraces the above mentioned singularities of
the $g(w)$ function in the complex plane running in the clockwise direction.
An evaluation of this integral (see in Appendix D) results in the following
expression

$\displaystyle c\,G^{(0)}(z,z^{\prime};\omega)=$
$\displaystyle\frac{\omega}{\bar{w}}F(z,z^{\prime};\bar{w})$ (34)
$\displaystyle+$
$\displaystyle\frac{1}{\pi}\,\int_{\omega_{c}}^{\infty}d\bar{\omega}\frac{\sqrt{\bar{\omega}^{2}-\omega_{c}^{2}}}{\omega^{2}+\bar{\omega}^{2}-\omega_{c}^{2}}F(z,z^{\prime};\bar{\omega}),$
(35)

where the auxiliary function

$\displaystyle F(z,z^{\prime};w)$
$\displaystyle=\frac{e^{w\frac{|z-z^{\prime}|}{c}}-e^{w\frac{|z+z^{\prime}+L|}{c}}}{1-e^{\beta
w}}$
$\displaystyle-\frac{e^{-w\frac{|z-z^{\prime}|}{c}}-e^{-w\frac{|z+z^{\prime}+L|}{c}}}{1-e^{-\beta
w}}$ (36)

obeys the correct boundary conditions
$F(\pm\frac{L}{2},z^{\prime})=F(z,\pm\frac{L}{2})=0$.

The pole contribution (34) to $c\,G^{(0)}(z,z^{\prime};\omega)$ is dominant.
It features the exponentially decaying spatial behaviour in the stopband
$\omega<\omega_{c}$ ($\bar{w}$ is real), and the oscillating spatial behaviour
in the passband $\omega>\omega_{c}$ ($\bar{w}$ is imaginary).

The branch cut contribution (35) might be only important for $z\approx
z^{\prime}$ (otherwise the integrand is suppressed by quickly decaying
exponential terms). So we get an approximation for the equal-point Green’s
function

$\displaystyle c\,G^{(0)}(z,z;\omega)\approx$
$\displaystyle-\frac{\omega}{\bar{w}}-\frac{1}{\pi}\,\int_{\omega_{c}}^{\omega_{hf}}d\bar{\omega}\frac{\sqrt{\bar{\omega}^{2}-\omega_{c}^{2}}}{\omega^{2}+\bar{\omega}^{2}-\omega_{c}^{2}}$
(37) $\displaystyle=$
$\displaystyle-\frac{\omega}{\bar{w}}-\frac{1}{\pi}\ln\frac{2\omega_{hf}}{\omega_{c}}$
(38)
$\displaystyle+\frac{\omega}{\bar{w}}\frac{1}{\pi}\arcsin\frac{\bar{w}}{\omega_{c}}+O\left(\frac{\omega_{c}}{\omega_{hf}}\right).$
(39)

In the above integral (see in Appendix E for details of its evaluation) we
have introduced the high-frequency cutoff $\omega_{hf}\gg\omega_{c}$ in order
to regularize its logarithmic divergence (giving thus the leading branch cut
contribution which is constant in frequency). The calculated subleading
contribution is nearly constant ($\approx\frac{1}{\pi}$) near the threshold
$\omega\to\omega_{c}$ ($\bar{w}\to 0$). As it was discussed on the way from
(29) to (30), such constant contributions are negligible.

Overall we can omit the branch cut contribution (35) and reliably approximate
$c\,G^{(0)}(z,z^{\prime};\omega)$ by the pole contribution (34) alone. This
approximation gives consistent physical results in different limiting cases
considered below.

####  III.1.1 Long waveguide limit

In the long waveguide limit which is defined by the condition
$\bar{w}\frac{L}{c}\gg 1$ for $\omega<\omega_{c}$, the condition
$|\omega_{c}-\omega|\gg\frac{c^{2}}{L^{2}\omega_{c}}$ automatically holds and
the waveguide’s Green’s function can be approximated by

$\displaystyle
G^{(0)}(z,z^{\prime};\omega)\approx\frac{\omega}{c\bar{w}}F(z,z^{\prime};\bar{w})\approx-\frac{\omega}{c\bar{w}}e^{-\bar{w}\frac{|z-z^{\prime}|}{c}},$
(40)

that is by the translationally invariant expression for the infinite waveguide
length, which has no signatures of the endpoints.

####  III.1.2 Short waveguide limit

For a rather short waveguide’s length $L\sim\frac{c\pi}{\omega_{c}}$ (i.e. the
longitudinal and transverse sizes of the waveguide are comparable), such that
the condition
$|\omega-\varepsilon_{\frac{\pi}{L}}|\ll|\varepsilon_{\frac{2\pi}{L}}-\varepsilon_{\frac{\pi}{L}}|\sim\frac{c\pi}{L}$
is fulfilled, we expect that the expression (30) is well approximated by the
first term in the sum, i.e.

$\displaystyle G^{(0)}(z,z^{\prime};\omega)$
$\displaystyle\approx\frac{2}{L}\frac{\sin[\frac{\pi}{L}(z+\frac{L}{2})]\sin[\frac{\pi}{L}(z^{\prime}+\frac{L}{2})]}{\omega-\varepsilon_{\frac{\pi}{L}}}$
$\displaystyle=-\frac{1}{L\bar{\delta}}$
$\displaystyle\left[\cos\frac{\pi\,|z-z^{\prime}|}{L}-\cos\frac{\pi\,|z+z^{\prime}+L|}{L}\right],$
(41)

where $\bar{\delta}=\varepsilon_{\frac{\pi}{L}}-\omega$. Note that
$\bar{\delta}$ is defined with respect to the fundamental waveguide mode
$\varepsilon_{\frac{\pi}{L}}$, not with respect to $\omega_{c}$.

Let us show that the same result can be recovered from the term (34).
Approximating

$\displaystyle\bar{w}$
$\displaystyle=\sqrt{\omega_{c}^{2}-\omega^{2}}=\sqrt{\omega_{c}^{2}-(\varepsilon_{\frac{\pi}{L}}-\bar{\delta})^{2}}$
(42)
$\displaystyle\approx\sqrt{\omega_{c}^{2}-\varepsilon_{\frac{\pi}{L}}^{2}+2\varepsilon_{\frac{\pi}{L}}\bar{\delta}}\approx
i\frac{c\pi}{L}\left(1-\frac{L^{2}}{c^{2}\pi^{2}}\varepsilon_{\frac{\pi}{L}}\bar{\delta}\right),$
(43)

we insert this into (34) and obtain

$\displaystyle G^{(0)}(z,z^{\prime};\omega)\approx\frac{\omega L}{ic^{2}\pi}$
$\displaystyle\times$
$\displaystyle\left[\frac{e^{i\frac{\pi|z-z^{\prime}|}{L}}-e^{i\frac{\pi|z+z^{\prime}+L|}{L}}}{1-e^{-2\pi
i\frac{L^{2}}{c^{2}\pi^{2}}\varepsilon_{\frac{\pi}{L}}\bar{\delta}}}-\frac{e^{-i\frac{\pi|z-z^{\prime}|}{L}}-e^{-i\frac{\pi|z+z^{\prime}+L|}{L}}}{1-e^{2\pi
i\frac{L^{2}}{c^{2}\pi^{2}}\varepsilon_{\frac{\pi}{L}}\bar{\delta}}}\right]$
$\displaystyle\approx-\frac{\omega}{L\varepsilon_{\frac{\pi}{L}}\bar{\delta}}\left[\cos\frac{\pi|z-z^{\prime}|}{L}-\cos\frac{\pi|z+z^{\prime}+L|}{L}\right].$
(44)

To the leading order in $\bar{\delta}$, this expression is equivalent to (41).
The smallness condition on $\bar{\delta}$, justifying the expansions made
above, reads

$\displaystyle|\bar{\delta}|\ll\frac{c^{2}\pi^{2}}{L^{2}\varepsilon_{\frac{\pi}{L}}}\sim\frac{c\pi}{L}\sim|\varepsilon_{\frac{2\pi}{L}}-\varepsilon_{\frac{\pi}{L}}|,$
(45)

which accords with the initially made physical assumptions for this limiting
case.

###  III.2 Qubit-photon bound states below cutoff

####  III.2.1 One qubit — one bound state

For a single qubit $j=1$ in a 3D waveguide we obtain on the basis of Eq. (17)
an approximate transmission in the stopband by omitting the background empty
waveguide contribution:

$\displaystyle S_{RL}(\omega)$
$\displaystyle\approx-2i\sqrt{\Gamma_{R}\Gamma_{L}}g_{1}^{2}$
$\displaystyle\times\frac{L^{2}G^{(p)}(z_{R},z_{1};\omega)G^{(p)}(z_{1},z_{L};\omega)}{\omega-\omega_{q,1}+i\frac{\gamma_{a,1}}{2}-g_{1}^{2}LG^{(p)}(z_{1},z_{1};\omega)}.$
(46)

Next, we approximate $G^{(p)}$ by $G^{(0)}$ in the numerator, and $G^{(p)}$ by
(26) in the denominator. The latter reads then as

$\displaystyle\omega-\omega_{q,1}-\gamma\,c\,G^{(0)}(z_{1},z_{1};\omega)$ (47)
$\displaystyle+$ $\displaystyle i\left[\frac{\gamma_{a,1}}{2}+\frac{\gamma
L}{c}\sum_{s=L,R}\Gamma_{s}|c\,G^{(0)}(z_{1},z_{s};\omega)|^{2}\right],$ (48)

where $\gamma=\frac{g_{1}^{2}L}{c}$. Assuming a weak resonance broadening
quantified by (48), we obtain from (47) the qubit-photon bound state equation

$\displaystyle\omega_{b}-\omega_{q,1}-\gamma\,c\,G^{(0)}(z_{1},z_{1};\omega_{b})=0.$
(49)

Expanding the denominator near the bound state frequency $\omega_{b}$, and
setting $\omega=\omega_{b}$ in the numerator, we find the following Lorentzian
approximation for the bound state resonance

$\displaystyle S_{RL}(\omega)$
$\displaystyle\approx-2i\sqrt{\Gamma_{R}\Gamma_{L}}\frac{\gamma
L}{c}Z(\omega_{b})$
$\displaystyle\times\frac{c\,G^{(0)}(z_{R},z_{1};\omega_{b})\,c\,G^{(0)}(z_{1},z_{L};\omega_{b})}{\omega-\omega_{b}+i\frac{\Gamma_{b}}{2}},$
(50)

in terms of the linewidth $\Gamma_{b}$ which is a sum of both non-radiative
decay $\Gamma_{q}$ of the qubit and the radiative decay $\Gamma_{r}$,

$\displaystyle\Gamma_{b}\approx
2Z(\omega_{b})\left[\frac{\Gamma_{q,1}}{2}+\frac{\gamma
L}{c}\sum_{s=L,R}\Gamma_{s}|c\,G^{(0)}(z_{1},z_{s};\omega_{b})|^{2}\right]$
(51)

and the quasiparticle weight [34]

$\displaystyle
Z(\omega_{b})=\frac{1}{1-\gamma\frac{\partial}{\partial\omega}\,c\,G^{(0)}(z_{1},z_{1};\omega)}\bigg{|}_{\omega\to\omega_{b}}<1.$
(52)

For $Z\to 1$ the bound state has more weight in the qubit excited state, while
for $Z\ll 1$ it is mostly photon-like (see below for a justification of this
interpretation). Calculating this factor for our model, we obtain the result

$\displaystyle
Z(\omega_{b})=\frac{1}{1+\frac{\gamma\omega_{c}^{2}}{(\omega_{c}^{2}-\omega_{b}^{2})^{3/2}}},$
(53)

in which we see that the Markov approximation (consisting in the neglect of
the frequency dependence in the qubits’ self-energy) is justified for the
bound state frequency far below $\omega_{c}$ (provided that
$\gamma\ll\omega_{c}$), while for $\omega_{b}$ approaching $\omega_{c}$ the
Markov approximation breaks down.

Making the Lorentzian approximation for the transmission probability

$\displaystyle|S_{RL}(\omega)|^{2}\approx\frac{A(\Gamma_{b}/2)^{2}}{(\omega-\omega_{b})^{2}+(\Gamma_{b}/2)^{2}},$
(54)

we also define the amplitude at the resonance

$\displaystyle A=|S_{RL}(\omega_{b})|^{2}$
$\displaystyle\approx\Big{|}4Z(\omega_{b})\frac{\sqrt{\Gamma_{R}\Gamma_{L}}}{\Gamma_{b}}\frac{\gamma
L}{c}\Big{|}^{2}$ (55)
$\displaystyle\times\Big{|}c\,G^{(0)}(z_{1},z_{R};\omega_{b})\,c\,G^{(0)}(z_{1},z_{L};\omega_{b})\Big{|}^{2}.$

It is also worth mentioning that the representation (113) gives us an
immediate access to the photonic part $\psi_{phot}^{(1)}(z)$ of the bound
state wavefunction. Neglecting the port couplings, we identify

$\displaystyle|\psi_{phot}^{(1)}(z)|^{2}$
$\displaystyle=\lim_{\omega\to\omega_{b}}[(\omega-\omega_{b})G(z,z;\omega)]$
(56) $\displaystyle=Lg_{1}^{2}Z(\omega_{b})|G^{(0)}(z,z_{1};\omega_{b})|^{2}.$
(57)

Observing that

$\displaystyle\int_{-L}^{L}dz|\psi_{phot}^{(1)}(z)|^{2}$
$\displaystyle=-Z(\omega_{b})\,\gamma\,c\frac{\partial}{\partial\omega_{b}}G^{(0)}(z_{1},z_{1};\omega_{b})$
(58) $\displaystyle=1-Z(\omega_{b}),$ (59)

we justify the interpretation of $Z(\omega_{b})$ given after its definition in
Eq. (52). In the long waveguide limit we apply the approximation (40) to get

$\displaystyle|\psi_{phot}^{(1)}(z)|^{2}\propto
e^{-2\sqrt{\omega_{b}^{2}-\omega_{c}^{2}}\frac{|z-z_{1}|}{c}}=e^{-2\frac{|z-z_{1}|}{\xi(\omega_{b})}},$
(60)

that is $\psi_{phot}(z)$ is exponentially localized, with the localization
length

$\displaystyle\xi(\omega_{b})=\frac{c}{\sqrt{\omega_{c}^{2}-\omega_{b}^{2}}},$
(61)

near the qubit position $z_{1}$ — this spatial profile is actually sketched in
Fig. (2a).

####  III.2.2 Two qubits — two bound states

In the presence of the two qubits we have two bound-state resonances which can
be identified with the poles of the qubits’ Green’s function (18). To neglect
broadening of the resonances, we approximate
$\hat{\Sigma}_{jj^{\prime}}(\omega)$ in (20) by

$\displaystyle\hat{\Sigma}_{jj^{\prime}}(\omega)\approx\hat{\Sigma}_{jj^{\prime}}^{(0)}(\omega)=g_{j}g_{j^{\prime}}LG^{(0)}(z_{j},z_{j^{\prime}};\omega)$
(62)

as well as neglect non-radiative decay rates $\gamma_{a,j}$ in (19). Then the
bound states are found from the equation

$\displaystyle 0$
$\displaystyle=\det[\hat{G}(\omega)]^{-1}\approx\det\left([\hat{G}^{(0)}(\omega)]^{-1}-\hat{\Sigma}^{(0)}(\omega)\right)$
$\displaystyle\approx(\omega-\omega_{q,1}-\hat{\Sigma}_{11}^{(0)}(\omega))(\omega-\omega_{q,2}-\hat{\Sigma}_{22}^{(0)}(\omega))$
$\displaystyle-\hat{\Sigma}_{12}^{(0)}(\omega)\hat{\Sigma}_{21}^{(0)}(\omega).$
(63)

Due to the off-diagonal components $\hat{\Sigma}_{12}^{(0)}(\omega)$ and
$\hat{\Sigma}_{21}^{(0)}(\omega)$ of the qubits’ self-energy, which describe
an effective waveguide-mediated exchange interaction between the qubits, the
two bound states repel each other. The energy splitting between them
quantifies the strength $J$ of the exchange interaction, and it is the subject
of our next consideration.

For the symmetric setup with $z_{1}=-z_{2}=-\frac{d}{2}$ and $g_{1}=g_{2}$, we
have the identities
$\hat{\Sigma}_{11}^{(0)}(\omega)=\hat{\Sigma}_{22}^{(0)}(\omega)$ and
$\hat{\Sigma}_{12}^{(0)}(\omega)=\hat{\Sigma}_{21}^{(0)}(\omega)$. Tuning
$\omega_{q,2}$ to the value $\omega_{q,2}=\omega_{q,1}\equiv\omega_{q}$ in
order to have a minimally possible splitting between the two bound states, we
get two separate equations for each bound state

$\displaystyle\omega_{b,1}-\omega_{q}-\hat{\Sigma}_{11}^{(0)}(\omega_{b,1})-\hat{\Sigma}_{12}^{(0)}(\omega_{b,1})=0,$
(64)
$\displaystyle\omega_{b,2}-\omega_{q}-\hat{\Sigma}_{11}^{(0)}(\omega_{b,2})+\hat{\Sigma}_{12}^{(0)}(\omega_{b,2})=0.$
(65)

Defining the splitting $\Delta=\omega_{b,2}-\omega_{b,1}$ as well as the
middle point $\bar{\omega}_{b}=\frac{\omega_{b,2}+\omega_{b,1}}{2}$ we derive
the following self-consistent equation

$\displaystyle\Delta=\hat{\Sigma}_{-}^{(0)}\left(\bar{\omega}_{b}+\frac{\Delta}{2}\right)-\hat{\Sigma}_{+}^{(0)}\left(\bar{\omega}_{b}-\frac{\Delta}{2}\right),$
(66)

with
$\hat{\Sigma}_{\pm}^{(0)}(\omega)=\hat{\Sigma}_{11}^{(0)}(\omega)\pm\hat{\Sigma}_{12}^{(0)}(\omega)$,
which implicitly defines $\Delta$ as a function of $\bar{\omega}_{b}$.

An approximation to (66) by expanding its right-hand side up to the linear
order in $\Delta$ yields

$\displaystyle\Delta(\bar{\omega}_{b})$
$\displaystyle\approx-2\hat{\Sigma}_{12}^{(0)}(\bar{\omega}_{b})Z(\bar{\omega}_{b})\equiv
J.$ (67)

Alternatively, this expression can be represented (see Appendix F) as

$\displaystyle\Delta(\bar{\omega}_{b})$ $\displaystyle\approx
2|\langle\psi_{phot}^{(1)}|(\bar{\omega}_{b}-H_{0})|\psi_{phot}^{(2)}\rangle|,$
(68)

that is (two times) the overlap between the photonic contributions
$\psi_{phot}^{(1)}$ and $\psi_{phot}^{(2)}$ to the corresponding bound states
at their degeneracy point $\bar{\omega}_{b}$, subtracting the correction term
$\langle\psi_{phot}^{(1)}|H_{0}|\psi_{phot}^{(2)}\rangle$ to avoid the double-
counting of $H_{0}$ (since it is used each time for determining both
$\psi_{phot}^{(1)}$ and $\psi_{phot}^{(2)}$ independently of each other).

In the long waveguide limit $L\gg\xi(\bar{\omega}_{b})$ the coupling between
the two qubits is approximated by

$\displaystyle J$ $\displaystyle\approx
2\gamma\frac{\bar{\omega}_{b}\,\xi(\bar{\omega}_{b})}{c}e^{-\frac{d}{\xi(\bar{\omega}_{b})}}\frac{1}{1+\gamma\frac{\omega_{c}^{2}}{c^{3}}\xi^{3}(\bar{\omega}_{b})}.$
(69)

In the single-mode limit we obtain the coupling between the two qubits by
approximating (67) using (41) and (62):

$\displaystyle J\approx
2Z(\bar{\omega}_{b})\frac{\tilde{g}_{1}\tilde{g}_{2}}{\varepsilon_{\frac{\pi}{L}}-\bar{\omega}_{b}},$
(70)

where

$\displaystyle\tilde{g}_{j}=g_{j}\,\sqrt{L}\psi_{\frac{\pi}{L}}(z_{j}),$ (71)

and

$\displaystyle Z(\bar{\omega}_{b})$
$\displaystyle=\left[1+\frac{\tilde{g}_{j}^{2}}{(\varepsilon_{\frac{\pi}{L}}-\bar{\omega}_{b})^{2}}\right]^{-1}.$
(72)

In particular, $Z(\bar{\omega}_{b})\approx 1$ for
$|\varepsilon_{\frac{\pi}{L}}-\bar{\omega}_{b}|\gg|\tilde{g}_{j}|$.

Figure 3: Comparison of the Purcell decay rate of a qubit coupled to a cavity
with different number of modes $N_{r}$.

##  IV Generalized Purcell effect

The spontaneous emission rate of a qubit depends on both the nature of the
electromagnetic environment to which it is coupled as well as the detuning
between the frequency of the qubit and the modes of the environment. In our
formalism, the most general formula for calculating the Purcell induced decay
is given by,

$\Gamma_{r,N_{r}\rightarrow\infty}=2Z(\omega_{b})\frac{\gamma
L}{c}\sum_{s=L,R}\Gamma_{s}|c\,G^{(0)}(z_{1},z_{s};\omega_{b})|^{2}.$ (73)

The Green’s function $G^{(0)}$ carries the information about the environment
as seen by the qubit. For a rectangular waveguide that supports the infinite
number of modes, the Green’s function is given by Eqs. (34), (35).

Similarly, for a qubit coupled to a multimode environment with a finite number
of resonances $N_{r}$, such as a photonic crystal and a cavity array with
$N=N_{r}$ cavities, the Green’s function is represented as a discrete sum (see
an exemplary evaluation in the framework of the tight-binding model in
Appendix G). The decay rate can be then estimated by the following formula

$\Gamma_{r,N_{r}}=2Z(\omega_{b})\frac{\gamma
L}{c}\sum_{s=L,R}\Gamma_{s}\left|c\sum_{l=1}^{N_{r}}\frac{\psi_{k_{l}}(z_{1})\psi_{k_{l}}(z_{s})}{\omega_{b}-\varepsilon_{k_{l}}}\right|^{2}.$
(74)

In the case of a single mode cavity ($N_{r}=1$), where the qubit in only
coupled to the fundamental resonance $l=1$, we obtain the bound state decay
rate by approximating Eq. (73) with Eq. (41),

$\displaystyle\Gamma_{r,N_{r}=1}$ $\displaystyle\approx
Z(\omega_{b})\frac{\tilde{g}_{j}^{2}}{(\varepsilon_{\frac{\pi}{L}}-\omega_{b})^{2}}\,\kappa,$
(75)

where $\tilde{g}_{j}$ is defined in (71), and

$\displaystyle\kappa$
$\displaystyle=\sum_{s=L,R}2\Gamma_{s}\,|\sqrt{L}\psi_{\frac{\pi}{L}}(z_{s})|^{2}$
(76)

is the cavity decay rate. Taking into account (72), which holds in the single-
mode case, we observe that near $\varepsilon_{\frac{\pi}{L}}=\omega_{b}$ the
decay rate (75) is regularized,

$\displaystyle\Gamma_{r,N_{r}=1}$
$\displaystyle\approx\frac{\tilde{g}_{j}^{2}}{(\varepsilon_{\frac{\pi}{L}}-\omega_{b})^{2}+\tilde{g}_{j}^{2}}\,\kappa,$
(77)

and coincides with the cavity decay rate $\kappa$ at
$\varepsilon_{\frac{\pi}{L}}=\omega_{b}$. In turn, for
$|\varepsilon_{\frac{\pi}{L}}-\omega_{b}|\gg|\tilde{g}_{j}|$ it holds
$Z(\omega_{b})\approx 1$, and we obtain the well known formula [35] for
Purcell induced decay

$\displaystyle\Gamma_{r,N_{r}=1}$
$\displaystyle\approx\frac{\tilde{g}_{j}^{2}}{(\varepsilon_{\frac{\pi}{L}}-\omega_{b})^{2}}\,\kappa.$
(78)

In Fig. (3) we have compared the life time of a qubit coupled to a cavity with
different number of modes. Below the fundamental frequency, the propagating
modes are suppressed and the qubit predominantly couples to evanescent modes.
Therefore, in the limit of $N_{r}\rightarrow\infty$, the qubit has the largest
lifetime. In contrast, for frequencies above the fundamental mode the qubit
life is modified due to presence of higher harmonic modes which contributes to
the density of states available at the qubit frequency. Consequently the qubit
life time decreases with increase in the number of modes.

##  V Experiment

Figure 4: Transmission of the empty waveguide measured at room temperature
(black solid line). It is compared to the theoretical results computed using
two methods. The dot-dashed lines result from the numerical simulations of
$S_{RL}$ using the truncated sums in (30) at $N_{r}=10,\ 50\ \text{and}\ 100$
resonances. The dashed blue line is calculated using the exactly evaluated sum
in (30) (i.e. at $N_{r}\rightarrow\infty$). Exploiting its good agreement with
the experimental data, we estimate on the basis of its analytical expression
the values of $\Gamma_{L,R}/2\pi\approx 225$ MHz. The cut-off frequency is at
$\omega_{c}/2\pi=6.5213$ GHz.

###  V.1 Waveguide characterization

We first characterize the waveguide without qubits through a transmission
measurement at room temperature and fit the experimental data to theoretical
values given by Eq. (16). Note that the coupling between the waveguide and the
ports $\Gamma_{s=L,R}$ are the only free parameters in the fit. We also let
$\Gamma_{L}=\Gamma_{R}=\Gamma_{L,R}$, which is a reasonable assumption since
we use the identical co-axial connectors for both the ports. Furthermore, we
observe that the transmission amplitude reaches unity in the passband of the
waveguide indicating that there is no left-right asymmetry in the couplings.
In Fig. (4) we have compared the transmission data to two different
theoretical results. The first one is based on the approximation where the
waveguide’s Green’s function is represented by a sum of a finite number of
resonances $N_{r}$ truncating the sum in Eq. (30) — effectively this is
equivalent to replacing the continuum model by a lattice counterpart with
$N_{r}$ sites and the lattice constant $L/N_{r}$. In turn, the full sum
($N_{r}\rightarrow\infty$) provides the exact description of the continuum
model leading to the analytic expression Eq. (34) for the (bare) Green’s
function. We clearly see that the continuum limit best represents the data and
provides us with the estimate of $\Gamma_{L,R}/2\pi\approx 225$ MHz. We also
note that the small mismatch in fitting the passband can be caused by our
neglect of the weak frequency dependence of the coupling between the waveguide
modes and the ports which we made in deriving the formula for the
transmission. However, the minor effect of this couplings’ dispersion is
inessential in the stopband, where all subsequent measurements are performed.

Figure 5: (a) Measurement data and (b) is the simulation of waveguide
spectroscopy reveals the dressed QPBS appearing as a peak in the stopband. The
frequency of the QPBS is shifted lower by a self energy $\Sigma_{11}$ from the
qubit frequencies, which in turn appears as a dip. Note that the bound state
always exists in the stopband, for qubit frequencies close to the cutoff it is
photon-like ($Z\approx 0$) and for frequencies deep in the stopband it is
qubit-like ($Z\approx 1$). (c) and (d) are linewidth of the bound state
measured from spectroscopy data for qubits $Q_{1}$ and $Q_{2}$ placed
symmetrically from the center of the waveguide. (e) and (f) shows the
attenuation of peak amplitude as a function of bound state frequency which has
similar trend as that of linewidth. The only free parameter in the fitting is
$\Gamma_{L,R}/2\pi\approx 525$ MHz.

###  V.2 Probing QPBS in the stopband

In order to detect QPBS in our experiment, we performed a spectroscopy
measurement by sweeping the frequency $\omega_{q,1}$ of the first qubit
($Q_{1}$) through the cutoff into the stopband of the waveguide and observed
the scattering properties. Through out the experiment the frequency
$\omega_{q,2}$ of the second qubit ($Q_{2}$) is placed in the passband and
hence far detuned from $\omega_{q,1}$. As shown in Fig. (5a), when we tune the
bare qubit frequency $\omega_{q,1}$ close to $\omega_{c}$, we observe a peak
in transmission signal in the stopband as a result of the formation of a
single photon dressed bound state that is localized around the position
$z_{1}$ of the qubit and extends to the ends of the waveguide. The
localization length $\xi$ given by (61) depends on the bound state frequency
$\omega_{b}$, which differs from the qubit frequency by the self-energy term,
$\omega_{b}-\omega_{q,1}=\hat{\Sigma}_{11}(\omega_{b})\approx\hat{\Sigma}_{11}(\omega_{q,1})$.
Thus, $\xi$ is tuned in situ by changing $\omega_{q,1}$.

Although the density of states vanishes in the stopband, the bound state
acquires a finite linewidth due to the coupling to the input and output ports.
The strength of this coupling depends on $\xi$, therefore measuring linewidth
of the bound state serves as an indirect measure of the localization length.
We extracted linewidth as a function of the bound state frequency $\omega_{b}$
for both $Q_{1}$ and $Q_{2}$ from independent spectroscopy measurements. In
Figs. (5a,b), we compared experimental data with our theoretical estimate
given in Eq. (51) and find good agreement for both the qubits. Since we
measure the output of the waveguide through the right port, the extent of the
bound state induced by $Q_{1}$ which is close to the port is large compared to
the bound state induced by $Q_{2}$. This results in smaller linewidth for the
bound state localized around $Q_{2}$ which can be observed by comparing Fig.
(5a) and (5b). We also observe that when the frequency of the bound states are
$\omega_{b}/2\pi\approx 6.1$ GHz, the linewidth has the lowest value of
$\Gamma_{b}/2\pi\approx 1$ MHz for both the bound states, which is
approximately close to the bare linewidth ($\Gamma_{q}$) of the qubits. This
indicates that in the span of $\approx 400$ MHz from $\omega_{c}$ the bound
states are completely localized in our system.

Similar to the linewidth, the transmission amplitude of the bound state
resonance also depends on $\omega_{b}$. In Figs. (5e,f) we show the
attenuation of the bound state signal as function of the frequency and find a
good agreement between the theoretical prediction (see Eq. (55)) and the data
for both $Q_{1}$ and $Q_{2}$. The qubit-waveguide coupling strength $g_{1,2}$
needed for the fitting is obtained from the avoided-crossing data which will
be discussed in section (V.4). Therefore the fit obtained for both linewidth
and amplitude has no free parameters which shows the validity of the theory
for our experimental results.

###  V.3 QPBS as an effective two-level system

Inducing a localized photonic mode in the bandgap can also be achieved by
introducing dislocations in the unit cell of a photonic crystal. This
principle underlies the working of defect mode lasers where the localized mode
acts as a high finesse cavity [36, 37]. In contrast, QPBS acts as a non-linear
defect as it inherits the anharmonic nature of the bare qubit seeding it. In
order to probe the internal structure of the bound state, we detuned the qubit
by $450$ MHz below the cutoff frequency such that the bound state frequency is
at $\omega_{b}/2\pi=6.077$ GHz and the corresponding linewidth is measured to
be $\Gamma_{b}/2\pi=0.143$ MHz. When the qubit frequency lies deep in the
stopband one can make the Markovian approximation to eliminate the waveguide
modes and treat the bound state as an effective two-level system with finite
anharmonicity (see Appendix H).

When probed with a microwave tone, the bound state responds by emitting photon
out of both the input $c_{in}$ and the output $c_{out}$ ports. Since the input
signal is strongly attenuated in the stopband, measurement of the radiation
coming out of the output port $\langle c_{out}\rangle$ will be a direct
measure of the bound state polarization $\langle\sigma_{-}^{(1)}\rangle$ as
shown in Eq. (164). This should be contrasted to the measurement in the
passband where the radiation will be a sum of the incoming radiation and the
qubit emission described in Eq. (162) and (163). The large spacial separation
between the ports exponentially suppresses the coupling through the waveguide
modes and minimizes direct cross-talk between the ports enabling direct
measurement of the bound state polarization even in the presence of the probe
field.

To elucidate the qubit nature of the QBPS we apply an excitation pulse of
fixed length $\tau_{p}=1.3\ \mu s$ to input port while measuring the emission
from the output port. In Fig. 6(a) one can see the emission as function of
time and driving frequency $\omega_{d}$ exhibiting the characteristic chevron
pattern of Rabi oscillations (at frequency $\Omega_{r}$) for the pulse
duration. More specifically, Fig. 6(b) shows a horizontal line cut of which
shows decaying Rabi oscillations for the duration of the pulse $\tau_{p}$ and
decays exponential to zero after that. We also notice that, when the drive
frequency becomes resonant with the bound state frequency $\omega_{b}$, the
Rabi oscillations vanishes and the emission has a constant value for the time
period $\tau_{p}$ of the drive. The amplitude of this emission depends the
drive strength, dissipation rate and the detuning
$\delta_{d}=\omega_{d}-\omega_{b}$ as shown in Eq. (167). The steady state
thus obtained is not unique to bound states but a simple feature of a driven-
dissipate two level system which cannot be reveled in a typical circuit-QED
setup due to the absence of direct coupling to $\langle\sigma_{-}\rangle$
component of the qubit. The numerical simulation of the Lindblad master
equation for the qubit reduced density matrix
$\frac{d\rho_{q}}{dt}=-i[H_{\text{eff}},\rho_{q}]+\Gamma_{b}\left(\sigma_{-}^{(1)}\rho_{q}\sigma_{+}^{(1)}-\frac{1}{2}\\{\sigma_{+}^{(1)}\sigma_{-}^{(1)},\rho_{q}\\}\right)$,
with
$H_{\text{eff}}=-\frac{\delta_{d}}{2}\sigma_{z}^{(1)}+\frac{\Omega_{r}}{2}\sigma_{x}^{(1)}$,
which is equivalent to (165) and (166) upon identifying
$\langle\hat{\sigma}^{(1)}_{\pm,z}(t)\rangle=\text{Tr}[\sigma^{(1)}_{\pm,z}\rho_{q}(t)]$,
reproduces all features that we observe in the experimental data as shown in
Fig. (6).

Figure 6: (a) Chevron pattern of the Rabi oscillations obtained by sweeping
drive frequency in the vicinity of the bound state frequency
$\omega_{b}/(2\pi)=6.077$ GHz. The length of the drive pulse is $\tau_{p}=1.3\
\mu s$ during which we observe oscillations. The signal decays at the rate
$\propto 2\Gamma_{b}$ after the pulse is turned off. (b) Numerical simulation
of the Chevron pattern with $\Gamma_{b}/\Omega_{r}\sim 0.28$ reproduces
qualitatively the features in the experimental data. (c) Horizontal line-cut
(orange dashed line) of (a). The signal recorded here is proportional to
direct emission form the bound state. Figure 7: (a) Avoided crossing measured
by fixing qubit $\omega_{q_{1}}$ and sweeping qubit $\omega_{q_{2}}$, when
both the bound states are on resonance the coupling is $J=60$ MHz. (b) Same as
(a) but the bound state is further detuned from the cutoff frequency, here
$J=15$ MHz. (c) and (d) are the simulation of the spectroscopy experiment that
reproduces all the features, the dip in transmission happens when the bound
states comes into resonance with the qubit. (e) and (f) shows the avoided
crossing measured at different bare bound state frequencies $\bar{\omega}_{b}$
for two different distance of separation between the qubits. For $d=44.95$ mm,
the formula obtained from the infinite waveguide limit gives an accurate
description of the avoided crossings. We use this the data to estimate the
qubit-waveguide coupling $g_{1}/2\pi=g_{2}/2\pi$ to be $62.53$ MHz. When the
inter-qubit distance is reduced to $d=22.95$ mm we notice that the approximate
formula starts to breaks down for avoided crossing measured close to the
cutoff.

###  V.4 Non-radiative dipole-dipole interaction between two QPBS

So far in our experiment, we have considered the formation of a single QPBS
centered at the qubit co-ordinate when $\omega_{q}<\omega_{c}$. However, many
interesting quantum effects in waveguide QED arise due to distinctive photon
mediated interaction between multiple quantum emitters. In this section, we
study the interaction between bound states that can be established by allowing
the individual photonic wavefunctions to overlap. Such type of interaction can
be described by transverse field exchange Hamiltonian of the form
$H_{I}=\frac{1}{2}J[\sigma_{+}^{(1)}\sigma_{-}^{(2)}+h.c.]$. In order to
measure the coupling $J$ between the qubits we tuned their transition
frequencies into the stopband such that two bound states are induced at the
qubits’ positions. The amount of overlap between the two bound states depends
on the detuning between the bound state frequencies and the cut-off frequency.
When both the bound states are on resonance a single virtual excitation is
exchanged between them leading to the formation of symmetric ($\ket{+}$) and
anti-symmetric ($\ket{-}$) pair of states:
$\ket{\pm}=(\ket{ge}\pm\ket{eg})/\sqrt{2}$, where $\ket{g}$ and $\ket{e}$ are
the ground and first exited levels of the bound states. We measured the
coupling $J$ with respect to the ”bare” bound state frequency
$\tilde{\omega}_{b}=(\omega_{b,1}+\omega_{b,2})/2$ through spectroscopy
measurements shown in Fig. (7). When detuning
$\delta/2\pi=(\omega_{c}-\tilde{\omega}_{b})/2\pi$ is $50$ MHz, the measured
coupling strength $J/2\pi=66$ MHz as shown in Fig. (7a). In comparison, the
size of the avoided crossing reduces to $J/2\pi=15$ MHz when the qubit is
tuned deep into the stopband where the detuning $\delta/2\pi=15$ MHz as shown
in Fig. (7b). We also utilized the general two qubit transmission model
described in Eq. (17) and reproduced all the features that we observe in the
experiment as shown in Fig. (7a & b).

The coupling between two QPBS can be shown (see Eq. (69)) to decay
exponentially with $J\propto e^{-|z_{1}-z_{2}|/\xi}$ in the infinite-waveguide
limit (applicable for $L\gg\xi$), where $\xi$ is the in situ tunable
localization length, see Eq. (61). This is in contrast to interactions
mediated by propagating modes in passband which is oscillatory. In order to
show the tunable nature of this interaction, we have measured the size of the
avoided crossings when both the qubits are on-resonance as a function
$\tilde{\omega}_{b}$. Note that the interaction strength not only depend on
detuning but also on the distance $d$ of separation between the qubits. As
shown in Fig. (7), we measured the coupling for two different distances and
observed that $J$ scales exponentially as predicted by Eq. (69). Furthermore,
by comparing Fig. (7(e & f) we precisely find when the infinite waveguide
limit is a good approximation. In our system, when the inter-qubit distance is
large $d=44.95$ mm we find that infinite limit gives the accurate description
of the avoided crossings. On the other hand, it breaks down when the inter-
qubit distance is small $d=22.95$ mm. The only free parameters in the fitting
are the qubit-waveguide couplings $g_{1,2}$, we assume that both the qubits
have approximately same coupling strength $g_{1}\approx g_{2}=g$ and estimated
it to be $g/2\pi\approx 62.53$ MHz which we have also used in fitting Fig.
(5). It is important to note that although the exponential nature of the
interaction may indicate that this interaction would vanish in the long
waveguide limit, $\xi$ can be made equal to the length of the total system by
appropriately choosing $\delta$ and $d$, hence the interaction is effectively
long-range.

##  VI Summary and conclusions

In summary we developed a theory which allows us to calculate all the
parameters of interest for a system of two qubits coupled to a waveguide of
finite length. More specifically, knowing the waveguide dispersion relation
$\varepsilon_{k}$ and its wavefunctions $\psi_{k}(z)$ for a given waveguide
length $L$ it is possible to account for boundaries through the calculation of
the waveguide Green’s function given by Eq. (24). Our formulae then cover all
parameter regimes as well as properties of the qubits and waveguide above and
below cutoff.

For the systems like 1D arrays of coupled cavities or 1D photonic crystals the
Green’s function can be found by direct summation analytically (for a few
cavities/cells) or numerically (for a large number of cavities/cells). For the
3D waveguide we take into account an infinite number of modes analytically by
Eqs. (34)-(36).

Using the qubit-photon states below cutoff as a benchmark we show that the
infinite limit arises when $\sqrt{\omega_{c}^{2}-\omega^{2}}L/c\gg 1$. Given
that the localization length of the photonic part of a wavefunction of a bound
state is $\xi=c/\sqrt{\omega_{c}^{2}-\omega^{2}}$ we have a simple physical
interpretation: the infinite limit holds as far as the photonic wave-functions
of the qubit-photon bound state does not reach the boundaries of the
waveguide. In this limit the exponential spacial decay of the waveguide
Green’s function arises (40) which is then casted into exponential
localization of the photonic part of the qubit-photon bound state and the
experimentally observed exponential dependence of various parameters on the
qubits detuning.

It is also interesting to note that in the infinite limit the ports are
completely decouple from the problem and neither validity of the infinite
limit nor the physics below cutoff depends on the coupling strength to the
ports. Thus, the infinite limit will apply even for the case of the small
coupling where the transmission in the passband will not be flat and will show
sharp resonance transmission peaks.

Eqs. (34-36) can also reproduce the Jaynes-Cummings model behaviour for the
opposite short waveguide limit. However, it should be noted that that once the
detuning of the qubit from the fundamental mode increases the short cavity
limit will break down and more modes will be taken into account. It is
especially visible for our generalised Purcell formula (73) describing the
life-time of the qubit-bound state deviate from the conventional Purcell
formula (75) as the qubit detuning grows. We compare our theoretical results
to the experimental data where we coupled two superconducting transmon qubits
to a rectangular copper waveguide. Here we focused on the regime where the
qubit frequencies are tuned below the cutoff frequency in the stopband of the
3D waveguide. In this regime, both the linewidth and the amplitude of the
transmission peaks corresponding to the qubit-photon bound states below cutoff
as well as the exchange interaction strength between the two qubit-photon
bound states show exponential dependence as predicted by the infinite limit
formulae.

Overall, our results provide an important step for waveguide QED regime with
the missing rigorous justification of the infinite-size limit for finite-size
physical implementations. Furthermore, out 3D waveguide QED system serve as an
alternative platform for studying interacting qubit-photon bound states. In
the future work a high fidelity dispersive readout of the bound state can be
performed by inserting a planar readout resonator inside the waveguide[38],
thereby effectively realizing circuit-QED in our set-up. We have also shown
the tunable nature of interaction between two bound states through long range
non-radiative dipole-dipole interaction. The exponential nature of these
interaction can be potentially used for performing quantum simulations of
quantum spin-models in a 3D-architecture [39].

###### Acknowledgements.

We thank E. Shahmoon for fruitful discussions. The authors were supported by
the Australian Research Council Centre of Excellence for Engineered Quantum
Systems (EQUS, CE170100009).

##  Appendix A Heisenberg equations of motion and derivation of the expression
(10) for the transmission

With the help of the Hamiltonian (2) we derive the Heisenberg equations of
motion

$\displaystyle\frac{d}{dt}a_{k}(t)$ $\displaystyle=-i\varepsilon_{k}a_{k}(t)$
$\displaystyle-i\sqrt{L}\sum_{s=L,R}f_{k}\psi_{k}(z_{s})\sqrt{\frac{\Gamma_{s}}{\pi}}\int
d\omega_{s}c_{\omega_{s}}(t)$
$\displaystyle-i\sqrt{L}\sum_{j=1,2}f_{k}\psi_{k}(z_{j})g_{j}\sigma_{-}^{(j)}(t),$
(79) $\displaystyle\frac{d}{dt}c_{\omega_{s}}(t)$
$\displaystyle=-i\omega_{s}c_{\omega_{s}}(t)$
$\displaystyle-i\sqrt{\frac{\Gamma_{s}}{\pi}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{s})a_{k}(t),$
(80) $\displaystyle\frac{d}{dt}\sigma_{-}^{(j)}(t)$
$\displaystyle=-i\left(\omega_{q,j}-i\frac{\gamma_{a,j}}{2}\right)\sigma_{-}^{(j)}(t)$
$\displaystyle+i\sum_{k}g_{j}\sqrt{L}f_{k}\psi_{k}(z_{j})\sigma_{z}^{(j)}(t)a_{k}(t).$
(81)

Integrating the port field operators

$\displaystyle c_{\omega_{s}}(t)$
$\displaystyle=c_{\omega_{s}}(0)e^{-i\omega_{s}t}$ (82)
$\displaystyle-i\sqrt{\frac{\Gamma_{s}}{\pi}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{s})\int_{0}^{t}dt^{\prime}e^{-i\omega_{s}(t-t^{\prime})}a_{k}(t^{\prime}),$

we next evaluate

$\displaystyle\int d\omega_{R}c_{\omega_{R}}(t)$
$\displaystyle=-i\sqrt{\frac{\Gamma_{R}}{\pi}}\int
d\omega_{R}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})$
$\displaystyle\times\int_{0}^{t}dt^{\prime}e^{-i\omega_{R}(t-t^{\prime})}a_{k}(t^{\prime})$
(83)
$\displaystyle=-2i\sqrt{\pi\Gamma_{R}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})$
$\displaystyle\times\int_{0}^{t}dt^{\prime}\delta(t-t^{\prime})a_{k}(t^{\prime})$
(84)
$\displaystyle=-i\sqrt{\pi\Gamma_{R}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})a_{k}(t).$
(85)

Note that we used $\int_{0}^{t}dt^{\prime}\delta(t-t^{\prime})=\frac{1}{2}$,
and set $c_{\omega_{R}}(0)=0$, which means no input pulse in the right port.
In turn,

$\displaystyle\int d\omega_{L}c_{\omega_{L}}(t)$ $\displaystyle=\int
d\omega_{L}c_{\omega_{L}}(0)e^{-i\omega_{L}t}$
$\displaystyle-i\sqrt{\pi\Gamma_{L}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{L})a_{k}(t).$
(86)

For the input field we analogously obtain

$\displaystyle c_{in}(t)$ $\displaystyle=\int
d\omega_{L}c_{\omega_{L}}(t+0^{+})$ (87) $\displaystyle=\int
d\omega_{L}c_{\omega_{L}}(0)e^{-i\omega_{L}(t+0^{+})}$
$\displaystyle-i2\pi\sqrt{\frac{\Gamma_{L}}{\pi}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{L})$
$\displaystyle\times\int_{0}^{t}dt^{\prime}\delta(t+0^{+}-t^{\prime})a_{k}(t^{\prime}).$
(88)

The last term vanishes, since $t^{\prime}=t+0^{+}$ lies outside the
integration range, and thus obtain

$\displaystyle c_{in}(t)=\int d\omega_{L}c_{\omega_{L}}(0)e^{-i\omega_{L}t}.$
(89)

Next, we define the output field

$\displaystyle c_{out}(t)$ $\displaystyle=\int
d\omega_{R}c_{\omega_{R}}(t-0^{+})$ (90)
$\displaystyle=-i2\pi\sqrt{\frac{\Gamma_{R}}{\pi}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})$
$\displaystyle\times\int_{0}^{t}dt^{\prime}\delta(t-0^{+}-t^{\prime})a_{k}(t^{\prime})$
(91)
$\displaystyle=-i2\sqrt{\pi\Gamma_{R}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})a_{k}(t).$
(92)

Note that in contrast to (85) we get now the twice larger contribution, since
$t^{\prime}=t-0^{+}<t$ entirely lies inside the integration range.

Performing the Fourier transform

$\displaystyle\tilde{a}_{k}(\omega)=\int_{0}^{\infty}dta_{k}(t)e^{i(\omega+i0^{+})t},$
(93)

we obtain

$\displaystyle\tilde{c}_{out}(\omega)=-2i\sqrt{\pi\Gamma_{R}}\sum_{k}\sqrt{L}f_{k}\psi_{k}(z_{R})\tilde{a}_{k}(\omega).$
(94)

Inserting the solutions (85) and (86) into (79), we obtain

$\displaystyle\frac{d}{dt}a_{k}(t)$
$\displaystyle=-i\varepsilon_{k}a_{k}(t)-i\sum_{k^{\prime}}\Sigma_{kk^{\prime}}^{(p)}a_{k^{\prime}}(t)$
(95)
$\displaystyle-i\sum_{j=1,2}\sqrt{L}f_{k}\psi_{k}(z_{j})g_{j}\sigma_{-}^{(j)}(t)$
(96)
$\displaystyle-i\sqrt{L}f_{k}\psi_{k}(z_{L})\sqrt{\frac{\Gamma_{L}}{\pi}}c_{in}(t),$
(97)

where the port-induced self-energy $\Sigma_{kk^{\prime}}^{(p)}$ is defined in
(14). In the Fourier representation this equation reads

$\displaystyle(\omega-\varepsilon_{k})\tilde{a}_{k}(\omega)-\sum_{k^{\prime}}\Sigma_{kk^{\prime}}^{(p)}\tilde{a}_{k^{\prime}}(\omega)$
$\displaystyle-\sqrt{L}\sum_{j=1,2}f_{k}\psi_{k}(z_{j})g_{j}\tilde{\sigma}_{-}^{(j)}(\omega)$
$\displaystyle=\sqrt{\frac{\Gamma_{L}}{\pi}}\sqrt{L}f_{k}\psi_{k}(z_{L})\tilde{c}_{in}(\omega).$
(98)

In the equation (81) we approximate $\sigma_{z}^{(j)}(t)\approx-1$, which is
justified for a weak power of the incident pulse. Thus we get the equations

$\displaystyle\frac{d}{dt}\sigma_{-}^{(j)}(t)$
$\displaystyle=-i\left(\omega_{q,j}-i\frac{\gamma_{a,j}}{2}\right)\sigma_{-}^{(j)}(t)$
(99) $\displaystyle-i\sum_{k}g_{j}\sqrt{L}f_{k}\psi_{k}(z_{j})a_{k}(t).$ (100)

Rewriting it in the Fourier representation, we eliminate
$\tilde{\sigma}_{-}^{(j)}(\omega)$:

$\displaystyle\tilde{\sigma}_{-}^{(j)}(\omega)$
$\displaystyle=\frac{1}{\omega-\omega_{q,j}+i\frac{\gamma_{a,j}}{2}}\sum_{k}g_{j}\sqrt{L}f_{k}\psi_{k}(z_{j})\tilde{a}_{k}(\omega).$
(101)

Inserting these relations into (98) yields

$\displaystyle\sum_{k^{\prime}}[(\omega-\varepsilon_{k})\delta_{kk^{\prime}}-\Sigma_{kk^{\prime}}(\omega)]\tilde{a}_{k^{\prime}}(\omega)$
$\displaystyle=\sqrt{\frac{\Gamma_{L}}{\pi}}\sqrt{L}f_{k}\psi_{k}(z_{L})\tilde{c}_{in}(\omega),$
(102)

where $\Sigma_{kk^{\prime}}(\omega)$ is the full self-energy of the waveguide
defined in (13).

Introducing the fully dressed waveguide Green’s function (12), we solve (102)
for $\tilde{a}_{k}(\omega)$:

$\displaystyle\tilde{a}_{k}(\omega)=\sqrt{\frac{\Gamma_{L}}{\pi}}\sum_{k^{\prime}}G_{kk^{\prime}}(\omega)\sqrt{L}f_{k^{\prime}}\psi_{k^{\prime}}(z_{L})\tilde{c}_{in}(\omega).$
(103)

Inserting the result into (94), we obtain the relation (9) between the input
and the output fields in terms of the transmission amplitude

$\displaystyle S_{RL}(\omega)$
$\displaystyle=-2iL\sqrt{\Gamma_{R}\Gamma_{L}}\sum_{k,k^{\prime}}f_{k}\psi_{k}(z_{R})G_{kk^{\prime}}(\omega)f_{k^{\prime}}\psi_{k^{\prime}}(z_{L}),$

which accords with (10).

##  Appendix B Transmission representation (16), (17)

To achieve the representation (16), (17) for the transmission amplitude
$S_{RL}(\omega)$, we attempt to solve the Dyson equation for the fully dressed
Green’s function $G_{kk^{\prime}}$ of the waveguide

$\displaystyle G_{kk^{\prime}}$
$\displaystyle=G_{kk^{\prime}}^{(p)}+\sum_{k^{\prime\prime},k^{\prime\prime\prime}}G^{(p)}_{kk^{\prime\prime}}\Sigma_{k^{\prime\prime}k^{\prime\prime\prime}}^{(q)}G_{k^{\prime\prime\prime}k^{\prime}}$
(104) $\displaystyle\equiv
G_{kk^{\prime}}^{(p)}+\sum_{j,j^{\prime}}\bar{Q}^{(p)}_{kj}\hat{G}_{jj^{\prime}}^{(0)}Q_{j^{\prime}k^{\prime}}.$
(105)

Hereby the port-dressed Green’s function $G_{kk^{\prime}}^{(p)}$ of the
waveguide is defined

$\displaystyle
G^{(p)}_{kk^{\prime}}(\omega)=\left(\frac{1}{\omega-\hat{\varepsilon}-\Sigma^{(p)}}\right)_{kk^{\prime}}$
(106)

in terms of the port-induced self-energy (14). The qubit-induced self-energy
$\Sigma_{kk^{\prime}}^{(q)}$ is given in (15). In addition, we define
$\tilde{\psi}_{k}(z)=f_{k}\psi_{k}(z)$,

$\displaystyle\bar{Q}^{(p)}_{kj}$
$\displaystyle=\sum_{k^{\prime\prime}}G^{(p)}_{kk^{\prime\prime}}\sqrt{L}\tilde{\psi}_{k^{\prime\prime}}(z_{j})g_{j},$
(107) $\displaystyle Q_{j^{\prime}k^{\prime}}^{(p)}$
$\displaystyle=g_{j^{\prime}}\sum_{k^{\prime\prime}}\sqrt{L}\tilde{\psi}_{k^{\prime\prime}}(z_{j^{\prime}})G_{k^{\prime\prime}k^{\prime}}^{(p)},$
(108)

and the analogous quantities without the superscript $(p)$. Finally, the bare
qubits’ Green’s function $\hat{G}^{(0)}_{jj^{\prime}}$ is given in (19).

Appropriately convoluting (104), we arrive at the equation

$\displaystyle
Q_{j^{\prime}k^{\prime}}=Q_{j^{\prime}k^{\prime}}^{(p)}+\sum_{j,j^{\prime\prime}}\hat{\Sigma}_{j^{\prime}j}\hat{G}_{jj^{\prime\prime}}^{(0)}Q_{j^{\prime\prime}k^{\prime}},$
(109)

where $\hat{\Sigma}_{jj^{\prime}}$ is the full qubits’ self-energy given in
(20).

Inverting the matrix in the qubits’ space in (109), we obtain

$\displaystyle
Q_{jk^{\prime}}=\sum_{j^{\prime}}\left(\frac{1}{1-\hat{\Sigma}\hat{G}^{(0)}}\right)_{jj^{\prime}}Q_{j^{\prime}k^{\prime}}^{(p)}.$
(110)

Inserting this expression into (105), we establish

$\displaystyle G_{kk^{\prime}}$
$\displaystyle=G_{kk^{\prime}}^{(p)}+\sum_{j,j^{\prime}}\bar{Q}^{(p)}_{kj}\left(\hat{G}^{(0)}\frac{1}{1-\hat{\Sigma}\hat{G}^{(0)}}\right)_{jj^{\prime}}Q_{j^{\prime}k^{\prime}}^{(p)}$
(111)
$\displaystyle=G_{kk^{\prime}}^{(p)}+\sum_{j,j^{\prime}}\bar{Q}^{(p)}_{kj}\hat{G}_{jj^{\prime}}Q_{j^{\prime}k^{\prime}}^{(p)},$
(112)

where $\hat{G}_{jj^{\prime}}$ is the fully dressed qubits’ Green’s function
given in (18). Hence, we find the coordinate-space solution

$\displaystyle G(z,z^{\prime};\omega)$
$\displaystyle=G^{(p)}(z,z^{\prime};\omega)$
$\displaystyle+L\sum_{j,j^{\prime}}g_{j}g_{j^{\prime}}G^{(p)}(z,z_{j};\omega)$
$\displaystyle\qquad\times\hat{G}_{jj^{\prime}}(\omega)G^{(p)}(z_{j^{\prime}},z^{\prime};\omega),$
(113)

which immediately leads to the representation (16), (17) of the transmission
amplitude (10).

##  Appendix C Solution of the Dyson equation (21)

We notice that the Eq. (21) is analogous to Eq. (104): The latter is obtain
from the former by the replacements

$\displaystyle G_{kk^{\prime}}\to G_{kk^{\prime}}^{(p)},$ (114) $\displaystyle
G_{kk^{\prime}}^{(p)}\to G^{(0)}_{kk^{\prime}},$ (115)
$\displaystyle\hat{G}_{jj^{\prime}}^{(0)}\to-i\pi\delta_{jj^{\prime}},\quad
j,j^{\prime}\to L,R,$ (116) $\displaystyle
g_{j}\to\sqrt{\frac{\Gamma_{j}}{\pi}},$ (117)
$\displaystyle\hat{\Sigma}_{jj^{\prime}}\to\frac{\sqrt{\Gamma_{j}\Gamma_{j^{\prime}}}}{\pi}LG^{(0)}(z_{j},z_{j^{\prime}}).$
(118)

Therefore we can exploit the solution (113), making in it the same
replacements. Most of them are obvious, and we comment only on the matrix
$\hat{G}_{jj^{\prime}}$: In the $(L,R)$ basis it reads

$\displaystyle\hat{G}=$
$\displaystyle\left(\begin{array}[]{cc}-\frac{1}{i\pi}-\frac{\Gamma_{L}}{\pi}LG^{(0)}(z_{L},z_{L})&-\frac{\sqrt{\Gamma_{L}\Gamma_{R}}}{\pi}LG^{(0)}(z_{L},z_{R})\\\
-\frac{\sqrt{\Gamma_{R}\Gamma_{L}}}{\pi}LG^{(0)}(z_{R},z_{L})&-\frac{1}{i\pi}-\frac{\Gamma_{R}}{\pi}LG^{(0)}(z_{R},z_{R})\end{array}\right)^{-1}$
(121) $\displaystyle=-\frac{i\pi}{D}$ (122)
$\displaystyle\\!\\!\times\left(\begin{array}[]{cc}1+iL\Gamma_{R}G^{(0)}(z_{R},z_{R})&-iL\sqrt{\Gamma_{L}\Gamma_{R}}G^{(0)}(z_{L},z_{R})\\\
-iL\sqrt{\Gamma_{R}\Gamma_{L}}G^{(0)}(z_{R},z_{L})&1+iL\Gamma_{L}G^{(0)}(z_{L},z_{L})\end{array}\right)$
(125)

together with $D$ defined in (23). Expanding the sums over $j$ and
$j^{\prime}$ in (113) we obtain (22).

##  Appendix D Further details on evaluation of the Matsubara sum (33)

To switch from the sum in the left-hand side of (33) to the integral in its
right-hand side, which is performed along the contours shown in Fig. (8), it
is also necessary to make sure that at $\text{Re}\,w\neq 0$ the integrand goes
to zero,

$\displaystyle\Big{|}\frac{g(w)}{1-e^{-\beta w}}\Big{|}\to 0,$ (126)

faster than $\frac{1}{|w|}$ as $|w|\to\infty$.

Excluding the special point $z=z^{\prime}$ (which we consider separately
later), we first inspect the case $\text{Re}\,w>0$ and observe the
exponentially decaying terms $e^{-|\text{Re}\,w|\frac{|z-z^{\prime}|}{c}}$,
$e^{-|\text{Re}\,w|\frac{|z+z^{\prime}+L|}{c}}$.

Figure 8: Integration contour for the integral in the right-hand side of (33)
in the complex plane of $w$. It embraces the poles $w=\pm\bar{w}$ of $g(w)$
(the red crosses) as well as the branch cuts of $g(w)$ (the branching points
$w=\pm\omega_{c}$ are indicated by the fat black points).

For $\text{Re}\,w<0$ we obtain the exponentially decaying terms
$e^{-|\text{Re}\,w|(\frac{2L}{c}-\frac{|z-z^{\prime}|}{c})}$,
$e^{-|\text{Re}\,w|(\frac{2L}{c}-\frac{|z+z^{\prime}+L|}{c})}$, provided that
$z$ and $z^{\prime}$ do not appear on the waveguide’s endpoints. The latter
condition is however not restrictive, since we know that the bare Green’s
function vanishes if one of its arguments approaches $+\frac{L}{2}$ or
$-\frac{L}{2}$.

Having established the sufficient conditions to deform the integration contour
to the shape shown in Fig. (8), we easily find the pole contribution (34) by
evaluating the corresponding residua values. In turn, to evaluate the branch
cut integrals we observe that

$\displaystyle
g(\bar{\omega}+i0^{+})-g(\bar{\omega}-i0^{+})=(e^{-\bar{\omega}\frac{|x-x^{\prime}|}{c}}-e^{-\bar{\omega}\frac{|x+x^{\prime}+L|}{c}})$
$\displaystyle\times\left[\frac{1}{\omega+i\,\text{sgn}(\bar{\omega})\sqrt{\bar{\omega}^{2}-\omega_{c}^{2}}}-\frac{1}{\omega-i\,\text{sgn}(\bar{\omega})\sqrt{\bar{\omega}^{2}-\omega_{c}^{2}}}\right].$
(127)

This observation leads to the expression (35).

In the special case $z=z^{\prime}$, the sum in the left-hand side of (33)
diverges logarithmically, and we regularize it by subtracting from it the
analogous sum with $\omega=0$. Defining the function

$\displaystyle\tilde{g}(w)$
$\displaystyle=\frac{1}{\omega-\sqrt{-w^{2}+\omega_{c}^{2}}}+\frac{1}{\sqrt{-w^{2}+\omega_{c}^{2}}}$
(128)
$\displaystyle=\frac{\omega}{\sqrt{-w^{2}+\omega_{c}^{2}}}\frac{1}{\omega-\sqrt{-w^{2}+\omega_{c}^{2}}},$
(129)

we see that it produces the same pole contribution (34) as $g(w)$ does at
$z=z^{\prime}$ (when neglecting the regular contribution $\propto
e^{-w\frac{|z+z^{\prime}+L|}{c}}$). Eventually this observation means that the
proposed regularization is equivalent to introducing the high-frequency cutoff
$\omega_{hf}$ in (37), which also captures the logarithmically diverging
frequency-independent contribution.

##  Appendix E Evaluation of the integral in (37)

Let us perform the integral in (37) by making the variable change
$\bar{\omega}=\omega_{c}\cosh\lambda$:

$\displaystyle-\frac{1}{\pi}\,\int_{\omega_{c}}^{\omega_{hf}}d\bar{\omega}\frac{\sqrt{\bar{\omega}^{2}-\omega_{c}^{2}}}{\omega^{2}+\bar{\omega}^{2}-\omega_{c}^{2}}$
$\displaystyle=-\frac{1}{\pi}\,\int_{0}^{\text{arccosh}\frac{\omega_{hf}}{\omega_{c}}}d\lambda\frac{\sinh^{2}\lambda}{(\frac{\omega}{\omega_{c}})^{2}+\sinh^{2}\lambda}$
(130)
$\displaystyle\approx-\frac{1}{\pi}\ln\frac{2\omega_{hf}}{\omega_{c}}+\frac{1}{\pi}\frac{\omega^{2}}{\omega_{c}^{2}}\int_{0}^{\infty}\frac{2d\lambda}{2(\frac{\omega}{\omega_{c}})^{2}-1+\cosh
2\lambda},$ (131)

where the last approximation is valid for $\omega_{hf}\gg\omega_{c}$.

The next change of variable $y=e^{2\lambda}$ leads to

$\displaystyle\int_{0}^{\infty}\frac{2d\lambda}{2(\frac{\omega}{\omega_{c}})^{2}-1+\cosh
2\lambda}$
$\displaystyle=\int_{-\infty}^{\infty}\frac{d\lambda}{2(\frac{\omega}{\omega_{c}})^{2}-1+\cosh
2\lambda}$ $\displaystyle=\int_{0}^{\infty}\frac{dy}{(y+\cos
2\alpha)^{2}+\sin^{2}2\alpha}$ $\displaystyle=\frac{2\alpha}{\sin
2\alpha}=\frac{\omega_{c}^{2}}{\omega\bar{w}}\arcsin\frac{\bar{w}}{\omega_{c}},$
(132)

where we have parameterized $\omega=\omega_{c}\cos\alpha$ and
$\bar{w}=\sqrt{\omega_{c}^{2}-\omega^{2}}=\omega_{c}\sin\alpha$. Combining
(132) with (131), we obtain the result shown in (39).

##  Appendix F Explanation of the formula (68)

Traditionally, the bound state in the presence of a single qubit is found from
the eigenvalue problem
$H|\psi_{BS}^{(1)}\rangle=\omega_{b}|\psi_{BS}^{(1)}\rangle$, where in the
Hamiltonian (2) the coupling to the ports and to the second qubit is
neglected. Representing

$\displaystyle|\psi_{BS}^{(1)}\rangle=\alpha^{(1)}|eg\rangle|\text{vac}\rangle+|gg\rangle\sum_{k}\beta_{k}^{(1)}a_{k}^{\dagger}|\text{vac}\rangle,$
(133)

where $|\text{vac}\rangle$ is the photonic vacuum state, we obtain the
following equations for the coefficients $\alpha^{(1)}$ and $\beta_{k}^{(1)}$:

$\displaystyle(\omega_{b}-\omega_{q,1})\alpha^{(1)}$
$\displaystyle=g_{1}\sum_{k}\sqrt{L}\psi_{k}(z_{1})\beta_{k}^{(1)},$ (134)
$\displaystyle(\omega_{b}-\varepsilon_{k})\beta_{k}^{(1)}$
$\displaystyle=g_{1}\sqrt{L}\psi_{k}(z_{1})\alpha^{(1)}.$ (135)

It follows

$\displaystyle|\psi_{BS}^{(1)}\rangle=\alpha^{(1)}\left[|eg\rangle|\text{vac}\rangle+|gg\rangle\sum_{k}\frac{g_{1}\sqrt{L}\psi_{k}(z_{1})}{\omega_{b}-\varepsilon_{k}}a_{k}^{\dagger}|\text{vac}\rangle\right].$
(136)

In particular, from this expression we recover the photonic part (57) of the
bound state wavefunction

$\displaystyle\psi_{phot}^{(1)}(z)$
$\displaystyle=\alpha^{(1)}\sum_{k}\frac{g_{1}\sqrt{L}\psi_{k}(z_{1})\psi_{k}(z)}{\omega_{b}-\varepsilon_{k}}$
(137) $\displaystyle=\sqrt{Z(\omega_{b})L}g_{1}G^{(0)}(z,z_{1};\omega_{b}),$
(138)

with $\alpha^{(1)}=\sqrt{Z(\omega_{b})}$.

Analogously we find a bound state emerging due to the coupling of the
waveguide only with the second qubit:

$\displaystyle|\psi_{BS}^{(2)}\rangle=\alpha^{(2)}\left[|ge\rangle|\text{vac}\rangle+|gg\rangle\sum_{k}\frac{g_{2}\sqrt{L}\psi_{k}(z_{2})}{\omega_{b}-\varepsilon_{k}}a_{k}^{\dagger}|\text{vac}\rangle\right].$
(139)

In the symmetric setup with $g_{1}=g_{2}$, $\omega_{q,1}=\omega_{q,2}$, and
$z_{1}=-z_{2}$, the energy of this bound state has the same value $\omega_{b}$
as for the state (136).

When the both qubits are coupled to the waveguide, we can find the energy
splitting of the two bound states from the overlap of the photonic
contributions to the states (136) and (139). In addition, we must subtract the
matrix element of $H_{0}$, given in (1), since this Hamiltonian has been
already used twice, that is in the eigenvalue problems for each bound state.
Thereby we get

$\displaystyle\langle\psi_{BS}^{(1)}|(\omega_{b}-H_{0})|\psi_{BS}^{(2)}\rangle=\langle\psi_{phot}^{(1)}|(\omega_{b}-H_{0})|\psi_{phot}^{(2)}\rangle$
(140)
$\displaystyle=Z(\omega_{b})g_{1}g_{2}L\sum_{k}\frac{\psi_{k}(z_{1})\psi_{k}(z_{2})}{\omega_{b}-\varepsilon_{k}}\equiv
Z(\omega_{b})\hat{\Sigma}_{12}^{(0)}(\omega_{b}).$

##  Appendix G Lattice realization of the waveguide

Let us model a cavity array or a photonic crystal by a chain of $N$ sites
(labelled by $n=1,\ldots,N$). Treating it in the tight-binding approximation
we introduce the nearest-neighbor hopping amplitude $-t$ and the uniform on-
site energy $\omega_{c}+2t$. The energy spectrum (see Fig. (9)) and the
eigenfunctions of this model are

$\displaystyle\varepsilon_{k}$ $\displaystyle=\omega_{c}+2t-2t\cos ka,$ (141)
$\displaystyle\psi_{k}(z_{n})$ $\displaystyle=\sqrt{\frac{2}{L}}\sin(kz_{n}),$
(142)

where $k=\frac{l\pi}{L}$ is labelled by integer $l$, $1\leq l\leq N$. Hereby
we introduced the chain’s length $L=(N+1)a$ in terms of the lattice constant
$a$, as well as the $n$th site coordinate $z_{n}=na$. We note the
normalization $a\cdot\sum_{n=1}^{N}\psi_{k}^{2}(z_{n})=1$.

All observables discussed in the paper are expressed via the core object — the
waveguide’s bare Green’s function (29). Its lattice analogue reads

$\displaystyle G^{(0)}(z_{n},z_{n^{\prime}};\omega)\equiv
G_{nn^{\prime}}^{(0)}(\omega)$ (143) $\displaystyle=$
$\displaystyle\frac{\omega}{\omega_{c}}\sum_{l=1}^{N}\frac{\psi_{\frac{l\pi}{L}}(na)\psi_{\frac{l\pi}{L}}(n^{\prime}a)}{\omega-\varepsilon_{\frac{l\pi}{L}}}-\frac{\delta_{nn^{\prime}}}{\omega_{c}a}$
(144) $\displaystyle\approx$
$\displaystyle\sum_{l=1}^{N}\frac{\psi_{\frac{l\pi}{L}}(na)\psi_{\frac{l\pi}{L}}(n^{\prime}a)}{\omega-\varepsilon_{\frac{l\pi}{L}}}.$
(145)

This finite sum can be easily evaluated numerically.

For large $N$, we can replace the sum by the integral

$\displaystyle G_{nn^{\prime}}^{(0)}(\omega)$
$\displaystyle\approx\frac{2}{\pi}\int_{0}^{\pi/a}dk\frac{\sin(kna)\sin(kn^{\prime}a)}{\omega-\omega_{c}-2t+2t\cos
ka}$ (146)
$\displaystyle=\frac{1}{a}\int_{-\pi}^{\pi}\frac{d\bar{k}}{2\pi}\frac{e^{i\bar{k}|n-n^{\prime}|}-e^{i\bar{k}|n+n^{\prime}|}}{\omega-\omega_{c}-2t+2t\cos\bar{k}}$
(147) $\displaystyle\equiv I_{|n-n^{\prime}|}-I_{|n+n^{\prime}|}.$ (148)

Figure 9: Dispersion relation of a waveguide modelled as cavity array. The the
detuning is given by $\delta$ and the size of the passband is $2t$.

Evaluating the integral $I_{\bar{n}}$ at integer $\bar{n}\geq 0$ and
$\omega_{c}-\omega\equiv\delta>0$ (in practice, the condition should be
$\delta\gg ta^{2}\frac{\pi^{2}}{L^{2}}$ in order to justify the integral
approximation (146)), we obtain

$\displaystyle I_{\bar{n}}$
$\displaystyle\equiv\frac{1}{a}\int_{-\pi}^{\pi}\frac{d\bar{k}}{2\pi}\frac{e^{i\bar{k}\,\bar{n}}}{\omega-\omega_{c}-2t+2t\cos\bar{k}}$
(149)
$\displaystyle=-\frac{1}{a\sqrt{(4t+\delta)\delta}}\left(\frac{2t}{2t+\delta+\sqrt{(4t+\delta)\delta}}\right)^{\bar{n}}.$
(150)

For the broad bandwidth $4t\gg\delta$ we further approximate

$\displaystyle I_{\bar{n}}$
$\displaystyle\approx-\frac{1}{a\sqrt{4t\delta}}e^{-\bar{n}\sqrt{\frac{\delta}{t}}},$
(151)

and then

$\displaystyle G_{nn^{\prime}}^{(0)}(\omega)\approx
I_{|n-n^{\prime}|}\approx-\frac{1}{a\sqrt{4t\delta}}e^{-|n-n^{\prime}|\sqrt{\frac{\delta}{t}}}.$
(152)

In general, the low-energy approximation
$\varepsilon_{k}\approx\omega_{c}+\frac{k^{2}}{2m_{\text{eff}}}$ in terms of
the effective mass $m_{\text{eff}}=\frac{1}{2ta^{2}}$ directly leads to (152).
This formula is obtained from the integral

$\displaystyle G_{nn^{\prime}}^{(0)}(\omega)$
$\displaystyle\approx-\int_{-\infty}^{\infty}\frac{dk}{2\pi}\frac{e^{ika|n-n^{\prime}|}}{\delta+\frac{k^{2}}{2m_{\text{eff}}}}$
(153)

after closing the integration contour in the upper half-plane of complex $k$
and evaluating the residue value at the pole
$k_{p}=i\sqrt{2m_{\text{eff}}\delta}$.

To draw an analogy of the present tight-binding model with the 3D waveguide
model (27), we note the value $m_{\text{eff}}=\frac{\omega_{c}}{c^{2}}$ for
the latter case. Thereby we also recover (40) at $\omega_{c}\gg\delta>0$,
which coincides with (152) (up to the relabelling of the parameters). This
observation implies that the bound state properties near the cutoff frequency
$\omega_{c}$ in both the 3D waveguide model (27) and the cavity array model
(141) are physically equivalent.

##  Appendix H Bound state as a two-level system

For a single qubit detuned deep in the stopband, we derive an effective
description in terms of the two-level system coupled to the ports by
eliminating the waveguide modes.

From (98) it follows

$\displaystyle\tilde{a}_{k}(\omega)=\sqrt{L}\sum_{j=1,2}\sum_{k^{\prime}}G^{(p)}_{kk^{\prime}}(\omega)f_{k^{\prime}}\psi_{k^{\prime}}(z_{j})g_{j}\tilde{\sigma}_{-}^{(j)}(\omega)$
$\displaystyle+\sqrt{\frac{\Gamma_{L}}{\pi}}\sqrt{L}\sum_{k^{\prime}}G^{(p)}_{kk^{\prime}}(\omega)f_{k^{\prime}}\psi_{k^{\prime}}(z_{L})\tilde{c}_{in}(\omega).$
(154)

Deep in the stopband we make the Markov approximation
$G^{(p)}_{kk^{\prime}}(\omega)\approx G^{(p)}_{kk^{\prime}}(\omega_{q,1})$.
After this we transform $\tilde{a}_{k}(\omega)$ back to the time domain, and
insert $a_{k}(t)$ into (100), relaxing the weak-intensity approximation
$\sigma_{z}^{(j)}\approx-1$.

For a single qubit $j=1$ this gives

$\displaystyle\frac{d}{dt}\sigma_{-}^{(1)}(t)$
$\displaystyle\approx-i\left(\omega_{q,1}-i\frac{\Gamma_{q,1}}{2}-\hat{\Sigma}_{11}(\omega_{q,1})\right)\sigma_{-}^{(1)}(t)$
(155)
$\displaystyle+ig_{1}\sigma_{z}^{(1)}(t)\sqrt{\frac{\Gamma_{L}}{\pi}}LG^{(p)}(z_{1},z_{L};\omega_{q,1})c_{in}(t).$
(156)

Averaging over an initial coherent state and denoting

$\displaystyle
2g_{1}\sqrt{\frac{\Gamma_{L}}{\pi}}LG^{(p)}(z_{1},z_{L};\omega_{q,1})\langle
c_{in}(t)\rangle=\Omega_{r}e^{-i\omega_{d}t}$ (157)

in terms of the Rabi frequency $\Omega_{r}$ and the driving frequency
$\omega_{d}$, we obtain the equation

$\displaystyle\frac{d}{dt}\langle\sigma_{-}^{(1)}(t)\rangle$
$\displaystyle\approx-i\left(\omega_{b}-i\frac{\Gamma_{b}}{2}\right)\langle\sigma_{-}^{(1)}(t)\rangle$
(158)
$\displaystyle+\frac{i}{2}\Omega_{r}e^{-i\omega_{d}t}\langle\sigma_{z}^{(1)}(t)\rangle.$
(159)

In addition, we derive

$\displaystyle\frac{d}{dt}\langle\sigma_{z}^{(1)}(t)\rangle$
$\displaystyle\approx-\Gamma_{b}[1+\langle\sigma_{z}^{(1)}(t)\rangle]$ (160)
$\displaystyle-i\Omega_{r}\langle
e^{-i\omega_{d}t}\sigma_{+}^{(1)}(t)-e^{i\omega_{d}t}\sigma_{-}^{(1)}(t)\rangle.$
(161)

Inserting the approximate $a_{k}(t)$ into (92) we obtain the effective input-
output relation

$\displaystyle c_{out}(t)$
$\displaystyle=-2i\sqrt{\Gamma_{R}\Gamma_{L}}LG^{(p)}(z_{R},z_{L};\omega_{b})c_{in}(t)$
(162)
$\displaystyle-2ig_{1}\sqrt{\pi\Gamma_{R}}LG^{(p)}(z_{R},z_{1};\omega_{b})\sigma_{-}^{(1)}(t).$
(163)

In the stopband the value of $G^{(p)}(z_{R},z_{L};\omega_{b})$ is negligible,
therefore $\langle\sigma_{-}^{(1)}(t)\rangle$ is a direct measure of $\langle
c_{out}(t)\rangle$. In the co-rotating frame, where
$\langle\hat{\sigma}_{-}^{(1)}(t)\rangle=\langle\sigma_{-}^{(1)}(t)\rangle
e^{i\omega_{d}t}$,
$\langle\hat{\sigma}_{z}^{(1)}(t)\rangle=\langle\sigma_{-}^{(1)}(t)\rangle$,
and $\langle\hat{c}_{out}(t)\rangle=\langle c_{out}(t)\rangle
e^{i\omega_{d}t}$, we relate

$\displaystyle\langle\hat{c}_{out}(t)\rangle\approx-2ig_{1}\sqrt{\pi\Gamma_{R}}LG^{(p)}(z_{R},z_{1};\omega_{b})\langle\hat{\sigma}_{-}^{(1)}(t)\rangle.$
(164)

As well we get the following equations for the qubit observables

$\displaystyle\frac{d}{dt}\langle\hat{\sigma}_{-}^{(1)}(t)\rangle$
$\displaystyle\approx
i\left(\delta_{d}+i\frac{\Gamma_{b}}{2}\right)\langle\hat{\sigma}_{-}^{(1)}(t)\rangle+\frac{i}{2}\Omega_{r}\langle\hat{\sigma}_{z}^{(1)}(t)\rangle,$
(165) $\displaystyle\frac{d}{dt}\langle\hat{\sigma}_{z}^{(1)}(t)\rangle$
$\displaystyle\approx-\Gamma_{b}[1+\langle\hat{\sigma}_{z}^{(1)}(t)\rangle]-i\Omega_{r}\langle\hat{\sigma}_{+}^{(1)}(t)-\hat{\sigma}_{-}^{(1)}(t)\rangle,$
(166)

where $\delta_{d}=\omega_{d}-\omega_{b}$. On their basis we find the steady
state value

$\displaystyle|\langle\hat{\sigma}_{-}^{(1)}\rangle_{ss}|=|\langle\sigma_{-}^{(1)}\rangle_{ss}|=\frac{\Omega_{r}\sqrt{4\delta_{d}^{2}+\Gamma_{b}^{2}}}{2\Omega_{r}^{2}+4\delta_{d}^{2}+\Gamma_{b}^{2}}.$
(167)

## References

  * Wallraff  _et al._ [2004] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nature 431, 162 (2004).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Sundaresan  _et al._ [2015] N. M. Sundaresan, Y. Liu, D. Sadri, L. J. SzHocs, D. L. Underwood, M. Malekakhlagh, H. E. Türeci, and A. A. Houck, Physical Review X 5, 021035 (2015).
  * Chakram  _et al._ [2021] S. Chakram, A. E. Oriani, R. K. Naik, A. V. Dixit, K. He, A. Agrawal, H. Kwon, and D. I. Schuster, Physical Review Letters 127, 107701 (2021).
  * Martínez  _et al._ [2019] J. P. Martínez, S. Léger, N. Gheeraert, R. Dassonneville, L. Planat, F. Foroughi, Y. Krupko, O. Buisson, C. Naud, W. Hasch-Guichard, S. Florens, I. Snyman, and N. Roch, npj Quantum Information 5, 1 (2019).
  * 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, Science 327, 840 (2010).
  * Abdumalikov  _et al._ [2011] A. A. Abdumalikov, O. V. Astafiev, Y. A. Pashkin, Y. Nakamura, and J. S. Tsai, Physical Review Letters 107, 043604 (2011).
  * Van Loo  _et al._ [2013] A. F. Van Loo, A. Fedorov, K. Lalumiere, B. C. Sanders, A. Blais, and A. Wallraff, Science 342, 1494 (2013).
  * Lalumière  _et al._ [2013] K. Lalumière, B. C. Sanders, A. F. van Loo, A. Fedorov, A. Wallraff, and A. Blais, Physical Review A 88, 043806 (2013).
  * Kimble [2008] H. J. Kimble, Nature 453, 1023 (2008).
  * Sipahigil  _et al._ [2016] A. Sipahigil, R. E. Evans, D. D. Sukachev, M. J. Burek, J. Borregaard, M. K. Bhaskar, C. T. Nguyen, J. L. Pacheco, H. A. Atikian, C. Meuwly, R. M. Camacho, F. Jelezko, E. Bielejec, H. Park, M. Lončar, and M. D. Lukin, Science 354, 847 (2016).
  * Pichler  _et al._ [2015] H. Pichler, T. Ramos, A. J. Daley, and P. Zoller, Physical Review A 91, 042116 (2015).
  * Hamann  _et al._ [2018] A. R. Hamann, C. Müller, M. Jerger, M. Zanner, J. Combes, M. Pletyukhov, M. Weides, T. M. Stace, and A. Fedorov, Physical Review Letters 121, 123601 (2018).
  * Hoi  _et al._ [2011] I.-C. Hoi, C. M. Wilson, G. Johansson, T. Palomaki, B. Peropadre, and P. Delsing, Physical Review Letters 107, 073601 (2011).
  * Bykov [1975] V. P. Bykov, Soviet Journal of Quantum Electronics 4, 861 (1975).
  * John and Wang [1990] S. John and J. Wang, Physical Review Letters 64, 2418 (1990).
  * Kofman  _et al._ [1994] A. Kofman, G. Kurizki, and B. Sherman, Journal of Modern Optics 41, 353 (1994).
  * Calajó  _et al._ [2016] G. Calajó, F. Ciccarello, D. Chang, and P. Rabl, Physical Review A 93, 033833 (2016).
  * Liu and Houck [2016] Y. Liu and A. A. Houck, Nature Physics 13, 48 (2016).
  * Sundaresan  _et al._ [2019] N. M. Sundaresan, R. Lundgren, G. Zhu, A. V. Gorshkov, and A. A. Houck, Physical Review X 9, 011021 (2019).
  * Scigliuzzo  _et al._ [2021] M. Scigliuzzo, G. Calajò, F. Ciccarello, D. P. Lozano, A. Bengtsson, P. Scarlino, A. Wallraff, D. Chang, P. Delsing, and S. Gasparinetti, arXiv:2107.06852 (2021).
  * Zhang  _et al._ [2022] X. Zhang, E. Kim, D. K. Mark, S. Choi, and O. Painter, arXiv: 2206.12803 (2022).
  * John and Quang [1996] S. John and T. Quang, Physical Review Letters 76, 1320 (1996).
  * Douglas  _et al._ [2015] J. S. Douglas, H. Habibian, C.-L. Hung, A. V. Gorshkov, H. J. Kimble, and D. E. Chang, Nature Photonics 9, 326 (2015).
  * Shi  _et al._ [2016] T. Shi, Y.-H. Wu, A. González-Tudela, and J. Cirac, Physical Review X 6, 021027 (2016).
  * Shi  _et al._ [2018] T. Shi, Y.-H. Wu, A. González-Tudela, and J. I. Cirac, New Journal of Physics 20, 105005 (2018).
  * Mirhosseini  _et al._ [2018] M. Mirhosseini, E. Kim, V. S. Ferreira, M. Kalaee, A. Sipahigil, A. J. Keller, and O. Painter, Nature Communications 9, 1 (2018).
  * Pozar [2011] D. M. Pozar,  _Microwave engineering_ (John wiley & sons, 2011).
  * Shahmoon and Kurizki [2013] E. Shahmoon and G. Kurizki, Physical Review A 87, 033831 (2013).
  * Paik  _et al._ [2011] H. Paik, D. I. Schuster, L. S. Bishop, G. Kirchmair, G. Catelani, A. P. Sears, B. R. Johnson, M. J. Reagor, L. Frunzio, L. I. Glazman, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Physical Review Letters 107, 240501 (2011).
  * Asenjo-Garcia  _et al._ [2017] A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, Physical Review A 95, 033818 (2017).
  * Schneider  _et al._ [2016] M. P. Schneider, T. Sproll, C. Stawiarski, P. Schmitteckert, and K. Busch, Physical Review A 93, 10.1103/physreva.93.013828 (2016).
  * Gardiner  _et al._ [2004] C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise_ (Springer Science & Business Media, 2004).
  * Mahan [2000] D. G. Mahan,  _Many-Particle Physics_ (Springer New York, NY, 2000).
  * Houck  _et al._ [2008] A. A. Houck, J. A. Schreier, B. R. Johnson, J. M. Chow, J. Koch, J. M. Gambetta, D. I. Schuster, L. Frunzio, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Physical Review Letters 101, 080502 (2008).
  * Yablonovitch [1987] E. Yablonovitch, Physical Review Letters 58, 2059 (1987).
  * Painter  _et al._ [1999] O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O'Brien, P. D. Dapkus, and I. Kim, Science 284, 1819 (1999).
  * Zoepfl  _et al._ [2017] D. Zoepfl, P. R. Muppalla, C. Schneider, S. Kasemann, S. Partel, and G. Kirchmair, AIP Advances 7, 085118 (2017).
  * Dalmonte  _et al._ [2015] M. Dalmonte, S. Mirzaei, P. Muppalla, D. Marcos, P. Zoller, and G. Kirchmair, Physical Review B 92, 174507 (2015).