# Few- and many-body photon bound states in quantum nonlinear media: a spin-
model approach

G. Calajó1 and D. E. Chang1,2 1ICFO-Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  2ICREA-Institució Catalana de Recerca i Estudis Avançats,
08015 Barcelona, Spain

###### Abstract

The emergence of multi-photon bound states in quantum nonlinear media is an
intriguing phenomenon that has attracted significant theoretical and
experimental interest. Such bound states can emerge in quite different
settings, such as arrays of quantum emitters coupled to photonic waveguides
and Rydberg nonlinear media. The theoretical tools to characterize the bound
states are equally diverse, ranging from Bethe ansatz to effective field
theories. Here, we propose and utilize a spin-model formulation of quantum
atom-light interactions as a general, unified framework to study bound states,
and apply this formalism to several distinct systems. We show how to obtain
the bound-state dispersion relation within the two-excitation subspace, and
also numerically investigate how the properties of higher multi-photon bound
states manifest themselves in the propagation of large photon number pulses
and spatio-temporal correlations at the output. Interestingly, for all the
specific systems studied (chiral and bi-directional waveguides, and Rydberg
media), we find that the large-photon number limit always coincides with the
well-known semi-classical soliton phenomenon of self-induced transparency or
immediate generalizations thereof.

##  I Introduction

One of the most distinct predictions within classical nonlinear optics is the
emergence of solitons, whose shape does not change during propagation Agrawal
. From early on, there were also efforts to understand how these solitons
might emerge from a fully quantum theory, perhaps most prominently in
continuous Kerr nonlinear media lai1 ; lai2 . While the weak nonlinearities of
conventional media historically made quantum properties largely inaccessible,
in recent years there have been diverse systems, ranging from ensembles of
Rydberg atoms firstrev ; callum ; first ; LiangBS  to waveguide QED systems
coupled to quantum emitters roy_rev_mod ; review_lodahl ; Darrick_rev_mod ;
shen_fan_prl ; shen_fan ; Harold_strcorr ; mahmo_calajo ; Harold_bs ;
Prasadnat , where quantum few-photon bound states can truly be observed due to
strong interactions at the level of individual photons changnl . The
theoretical techniques by which bound states can be predicted and studied in
these different systems is equally diverse, varying from Bethe ansatz
shen_fan_prl ; shen_fan ; Bethe ; sahandprl  (when the interactions between
photons are local in character) to effective field theories bienas ; Efimov ;
magrebi  (particularly successful for Rydberg nonlinear media).

In recent years, an alternative formalism to capture quantum atom-light
interactions has gained attention. This formalism is based on the insight that
the light-matter polaritons that form in near-resonant propagation are in fact
almost entirely atomic in character in typical settings. Thus, the photonic
degrees of freedom can be effectively integrated out, to yield a quantum spin
model describing photon-mediated interactions between the atoms anapra ;
anaprx ; ana_masson ; rev_wqedarray . In principle, if the dynamics of the
spin model can be solved, the correlations of the outgoing quantum fields can
be readily obtained using an input-output equation Caneva ; MPSJames ; Caneva
; MPSJames ; glauber ; Blais_in . Within this framework, particularly in the
context of atomic arrays and in waveguide QED, the spin model has led to the
prediction of many interesting effects, including strongly subradiant states
Albrecht ; Loic ; ritch ; dimermolmer1 ; zhang ; buchler ; dimer_shermet ;
Pod_scat ; Pod_3 ; Lesanosky ; Yudson , photon mediated localized states
Pod_photon_loc , topological states Pod_topo , chaotic states caos  and
subradiant dimers dimermolmer2 ; dimer_poddu . Techniques to effectively map
continuous, macroscopic nonlinear media such as Rydberg ensembles to the spin
model have also been proposed and investigated james_ryd .

Here, we show that the spin model formalism can be used to achieve a unified
description of photon bound states within atomic nonlinear media in general.
In particular, we exploit this formalism to describe both the quantum
nonlinear propagation of an optical field in an interacting Rydberg atom
ensemble, and in an array of two-level atoms coupled to a waveguide, treating
both paradigmatic cases of bidirectional and chiral Stannigel2012 ;
Pichler2015 ; ChiralRev  (unidirectional) emission. Within the spin model, the
procedure to distinguish bound states is analogous to that of identifying
bound magnon excitations in condensed matter spin models Bethe ; wortis ;
Schneider ; Mattis ; bloch , a well-known problem in the context of quantum
magnetism that has been recently reconsidered in open systems settings petroBS
; cooper ; ramos . An important physical distinction, however, is the long-
range nature of the photon-mediated spin interactions in our model, and also
the driven open nature of light-matter interacting systems. In particular,
using the input-output equations, we are able to show how these bound states
within the atomic medium are excited by non-interacting optical pulses
incident from the outside, and how these bound states manifest themselves in
the spatio-temporal correlations of the optical field once it leaves the
medium. These bound states are also closely related to the dimers previously
identified within waveguide QED dimermolmer2 ; dimer_poddu ; rev_wqedarray .
However, while the emphasis of those previous works was on the highly
subradiant nature of discrete dimer states, here, we focus on the continuous
dispersion relations of such states, which describe the distinct propagation
of these states from individual, non-interacting photons, and which are
responsible for the interesting spatio-temporal correlations that arise in the
output.

We find that the spin model is versatile, both able to capture bound state
properties in regimes where Bethe ansatz is not valid and without employing
approximations typical of effective theories, and readily suitable to treat
the many-photon regime numerically via matrix product state (MPS)
representations Verstraete ; Schollw . In this many-photon limit, we generally
find a transition from true quantum many-body behavior, as characterized by
non-trivial spatio-temporal correlations, to a semi-classical soliton wave
reminiscent of self-induced transparency in atomic media, irrespective of the
underlying details of the model.

This paper is structured as follows. In Sec. II we introduce a generic model
for an ensemble of emitters coupled to a quasi-one-dimensional (1D) quantum
optical field. From here, we formulate the spin model, and present general
criteria to identify the photon bound states. In Sec. III we consider as a
first case an array of two-level atoms (TLA) coupled to either a bidirectional
or chiral waveguide. We analyze the single- and two-excitation spectrum,
followed by multi-photon propagation. In particular, in the many-photon limit,
we recover the semi-classical phenomenon of self-induced transparency (SIT).
In Sec. IV, we apply the spin model to an ensemble of Rydberg atoms. After
analyzing the two-excitation bound state dispersion relation, we make use of
an effective two-level atom description to recover, in the many-body regime, a
solitonic behavior similar to SIT, i.e. _Rydberg SIT_. We conclude and provide
an outlook of possible future directions in Sec. V.

##  II Model

Figure 1: Generic quantum nonlinear media. A free-space ensemble of atoms with
an unspecified level structure is coupled to photons in a focused propagating
beam via a particular two-level transition. The system is characterized by a
single-atom emission rate $\Gamma$ into the beam, and $\Gamma_{0}$ into other
modes (free space). Alternatively, we consider an atom array coupled to a one-
dimensional optical waveguide atoms. An individual atom emits into the left-
propagating modes, right-propagating modes, and free space at rates
$\Gamma_{L},\Gamma_{R}$, and $\Gamma_{0}$. respectively. The free-space
ensemble can be mapped to an effective waveguide model, as we discuss in the
main text. In both scenarios, multi-photon bound states (red shaded area) can
exist due to the nonlinear interactions induced by the atomic medium.

We consider a setup as illustrated in Fig. 1 where an ensemble of $N$ atoms is
coupled to a quasi-1D propagating field. We assume that the field couples to a
transition of frequency $\omega_{eg}$ involving ground and excited states
$|g\rangle$ and $|e\rangle$, respectively, while the atoms can have other
levels (e.g., driven by classical control fields) that need not be specified
now. The spontaneous emission rate of a single, excited atom into the 1D modes
occurs at a rate $\Gamma$, while it can also emit at a rate $\Gamma_{0}$ into
other modes besides the 1D continuum. This generic setup describes well a
number of systems consisting of atoms coupled to a physical waveguide, where
one might have $\Gamma_{0}\gg\Gamma$ (such as for atoms coupled to optical
nanofibers) Prasadnat ; ReitzPRL2013 ; Mitsch2014 ; Pablo2017 ,
$\Gamma_{0}\sim\Gamma$ (for atoms coupled to photonic crystal waveguides)
Arcari2014 ; Sollner2015 ; Goban2014 ; Hood2016 , or $\Gamma_{0}\ll\Gamma$
(superconducting qubits coupled to unstructured or photonic crystal
transmission lines) Sundaresan ; painter1 ; painter2 ; ustinov ; ustinov_topo
; marcoBS . As pointed out in Ref. MPSJames , free space ensembles interacting
with quasi-1D optical beams can also be mapped to this waveguide model, by
taking $\Gamma_{0}\gg\Gamma$ (to account for the weak free-space light
coupling) and finding suitable maps between the microscopic parameters of this
model and the macroscopic parameters of the physical system such as optical
depth (discussed in detail in Sec. IV). For our situations of interest, we can
consider the waveguide-coupled atoms to be in an array of lattice constant
$d$, as shown in Fig. 1. For efficient coupling $\Gamma_{0}\ll\Gamma$, this
avoids Anderson or many-body localization of light associated with disorder
Anderson ; ReviewManybodyloc ; NikosManybodyloc , while for weak coupling
$\Gamma_{0}\gtrsim\Gamma$, multiple scattering of light is highly suppressed
and thus the physics is largely insensitive to choice of positions.

###  II.1 Full Hamiltonian

Within the 1D model the full waveguide QED Hamiltonian of the system is given
by the following contributions

$H=H_{ph}+H_{a}+H_{in}.$ (1)

The fist term $H_{ph}$ is the photonic Hamiltonian describing right- and left-
propagating modes travelling in the channel with group velocity $c$. It
explicitly reads ($\hbar=1$)

$H_{ph}=-i\int dx\left[\hat{E}_{R}^{\dagger}(x)\frac{\partial}{\partial
x}\hat{E}_{R}(x)-\hat{E}_{L}^{\dagger}(x)\frac{\partial}{\partial
x}\hat{E}_{L}(x)\right]$ (2)

where $\hat{E}_{R(L)}(x)$ is the bosonic field operator annihilating a right-
going (left-going) photon at position $x$ and fulfilling the commutation rule
$[\hat{E}_{R(L)}(x),\hat{E}^{\dagger}_{R(L)}(x^{\prime})]=c\delta(x-x^{\prime})$.
Note that with this definition of the field operator
$\langle\hat{E}^{\dagger}_{R(L)}(x)\hat{E}_{R(L)}(x)\rangle$ gives us directly
the intensity of the field. $H_{a}$ is the atomic Hamiltonian which at this
stage is not fully specified. In particular, we will assume that the atom has
ground and excited states $|g,e\rangle$ whose transition of frequency
$\omega_{eg}$ couples to the quantum propagating field. Beyond that, the atom
could contain additional levels, dissipation from other photonic channels
(encoded in the rate $\Gamma_{0}$), auxiliary control fields acting on
different transitions, or other terms independent from the interaction of
atoms with the quantum propagating field. The atom-light interaction allows
for the creation or annihilation of excitations on the $n$-th atom through the
spin operators $\hat{\sigma}^{n}_{ge}=|g_{n}\rangle\langle e_{n}|$ and
$\hat{\sigma}^{n}_{eg}=|e_{n}\rangle\langle g_{n}|$. The corresponding
interaction Hamiltonian reads

$H_{in}=\sum_{n}\left[\hat{\sigma}_{ge}^{n}\left(\sqrt{\Gamma_{R}}\hat{E}_{R}^{\dagger}(x_{n})+\sqrt{\Gamma_{L}}\hat{E}_{L}^{\dagger}(x_{n})\right)+\rm
H.c.\right]$ (3)

where $\Gamma_{R}$ and $\Gamma_{L}$, with $\Gamma=\Gamma_{R}+\Gamma_{L}$, are
the single-atom decay rates associated respectively to the emission of right-
and left-propagating photons. Here we explicitly distinguish the emission into
the two directions to cover both the paradigmatic cases of bi-directional
emission and chiral (uni-directional) emission that will be discussed in the
rest of the paper.

###  II.2 Spin model

Standard treatments to tackle Hamiltonian (1) in the multi-excitation sector
(multiple photons and/or excited atoms) usually follow two different
strategies depending on the regime considered. First, in the limit of strong
dissipation, $\Gamma_{0}\gg\Gamma$, individual atoms have negligible
interaction with light and interesting phenomena instead arise by collective
coupling. Then, the ensemble of atoms can be treated as a continuous bosonic
field. Interesting nonlinearities, such as arising from Rydberg interactions,
can be added and be treated by effective field theories (so far, limited to a
few excitations) bienas ; Efimov ; magrebi , for example. Conversely, in the
waveguide QED regime with negligible loss, $\Gamma_{0}\sim 0$, the eigenstates
can be exactly computed using scattering theory formalism (e.g. S-matrix,
Bethe ansatz, etc.) shen_fan_prl ; shen_fan ; Bethe ; sahandprl  for a few
excitations, either for small atom number in general or for large atom number
and chiral waveguides.

To go beyond specific limitations on excitation number or system details, we
consider an alternative approach, where the photons are integrated out to
arrive at an effective spin model for the atoms. This approximation takes
advantage of the fact that in most physical systems of interest, the atom-
photon dynamics occur on a time scale longer than the photon propagation time
through the system, i.e. $\Gamma\ll c/Nd$, so that the light-mediated
interactions can be considered instantaneous. Equivalently, the highly
dispersive nature of the atoms due to the small linewidth causes the light-
matter polaritons to be almost entirely atomic in character. The full system
evolution is then given by the master equation (ME) for the reduced atomic
density operator Caneva ; MPSJames ; glauber

$\dot{\rho}=-i\left[(H_{\rm eff}+H_{\rm drive})\rho-\rho(H_{\rm eff}+H_{\rm
drive})^{\dagger}\right]+\mathcal{J}[\rho].$ (4)

Here the non-Hermitian collective evolution of the system is given by the
effective Hamiltonian

$H_{\rm
eff}=H_{a}-i\sum_{nm}A_{nm}\hat{\sigma}^{n}_{eg}\hat{\sigma}^{m}_{ge},$ (5)

where the matrix

$A_{nm}=\left[\Gamma_{L}e^{ik_{0}|x_{n}-x_{m}|}\theta_{nm}+\Gamma_{R}e^{ik_{0}|x_{n}-x_{m}|}\theta_{mn}\right]$
(6)

encodes the photon mediated atom-atom interactions and
$\theta_{nm}=\theta(x_{n}-x_{m})$, with $\theta_{nn}=1/2$, is the Heaviside
function. Hamiltonian $H_{\rm eff}$ is invariant for translation of the
resonant wave-vector $k_{0}=\omega_{eg}/c$, i.e. $k_{0}d=k_{0}d+2\pi$, so a
full description is obtained by considering $k_{0}d\in[-\pi,\pi]$. The
population recycling contribution to the evolution in Eq.(4) is given by

$\mathcal{J}[\rho]=\Gamma_{0}\sum_{n}\hat{\sigma}^{n}_{ge}\rho\hat{\sigma}_{eg}^{n}+\sum_{nm}\left[\left(A_{nm}+A^{*}_{nm}\right)\hat{\sigma}^{n}_{ge}\rho\hat{\sigma}_{eg}^{m}+\rm
H.c.\right].$ (7)

Finally in Eq.(4) we also explicitly add a driving term

$H_{\rm
drive}=\sum_{n}\sqrt{\Gamma_{R}}\left[E_{in}(t,x_{n})\hat{\sigma}_{eg}^{n}+\rm
H.c.\right],$ (8)

which couples the emitters to a right-propagating coherent state input field
$E_{in}(t,x)=\mathcal{E}_{in}(t)e^{ik_{0}x-i\omega_{\rm in}t}$ with
$\omega_{\rm in}$ being the central frequency of the driving field.

While such master equations describing photon-mediated dipole-dipole
interactions have long existed, in recent years it has been realized that one
can also use their solutions to re-construct the quantum field that has been
previously integrated out, thus constituting a complete model of atom-light
interactions. This field takes the form of an input-output relation Caneva ;
MPSJames ; glauber ; Blais_in

$\hat{E}_{R}(t)=\mathcal{E}_{\rm
in}(t)+i\sum_{n}\sqrt{\Gamma_{R}}e^{-ik_{0}x_{n}}\hat{\sigma}_{ge}^{n}(t),$
(9)
$\hat{E}_{L}(t)=i\sum_{n}\sqrt{\Gamma_{L}}e^{ik_{0}x_{n}}\hat{\sigma}_{ge}^{n}(t).$
(10)

In this work, we will on one hand consider the microscopic properties of the
photon bound states themselves, which can be derived from the eigenstates of
(5) alone (practically in the few-excitation limit). Separately, we will solve
the full master equation dynamics of Eq. (4) in the presence of the input
field, to see how the bound states manifest themselves in the outgoing field
properties of Eq. (9). This dynamics can be computed even in the many-
excitation limit, with an MPS-based quantum trajectories algorithm discussed
in Appendix A.

###  II.3 Photon bound states

Photon bound states are states with spatially correlated positions, which
propagate through the bulk of a nonlinear medium experiencing low distortion
mahmo_calajo . Despite the name, these bound states are actually almost
entirely atomic (spin-like) in nature when they propagate through the medium
chang_mirror ; chang_tao . The problem of identifying these states then
becomes similar to that of magnon bound excitations in spin chains Bethe ;
wortis . For example, within the two-excitation subspace we can follow
standard procedures to diagonalize the Hamiltonian (5) in the relative and
center-of-mass coordinate frame, which is convenient because the latter is
characterized simply by a plane wave of wavevector $K$. On the other hand, the
problem in the relative coordinate reduces to finding single-particle bound
states in an effective (and possibly non-trivial) impurity model, whose
energies depend parametrically on $K$. Due to the nonlinear interaction, this
bound state will have a different energy $E_{\rm BS}$ than the sum of
individual photon energies, which allows for their spectral identification:

$E_{\rm BS}(K)\neq J(q)+J(K-q).$ (11)

Here, $J(k)$ denotes the single-excitation dispersion relation and $q$ the
relative momentum.

A crucial difference compared to canonical condensed matter spin models lies
in the dissipative, long-range nature of the interactions encoded in the non-
Hermitian Hamiltonian (5). A number of studies dimermolmer1 ; dimermolmer2 ;
dimer_shermet ; dimer_poddu  have already pointed out that in a finite system,
this Hamiltonian can give rise to two-excitation dimers as eigenstates that
are generally lossy, physically due to the non-zero support of the wave
functions with the system boundaries and their subsequent radiation into the
empty waveguide. While these studies largely focused on how long-lived or
subradiant these states could be, here, we more generally will study the
dispersion relation of this continuum of bound states, and show that this
dispersion relation and the spatial properties of the bound states indeed
manifest themselves in the quantum nonlinear optics problem, in terms of
correlations of the output field given a multi-photon input field.

##  III Array of two-level atoms

As a first example of a quantum nonlinear medium, we consider an array of two-
level atoms (TLA) where the atomic Hamiltonian is simply given by
$H_{a}=(\omega_{eg}-i\Gamma_{0}/2)\sum_{n}\hat{\sigma}_{eg}^{n}\hat{\sigma}_{ge}^{n}$.
In particular, we focus on the ideal waveguide QED scenario of low
dissipation, $\Gamma_{0}\rightarrow 0$. In this regime there is a non-
negligible probability that two propagating photons can interact
simultaneously with the same atom inducing a nonlinear optical response.
Motivated by the previous discussion, we will first calculate the single-
excitation dispersion relation, in order to identify energy gaps in which two-
excitation bound states might exist, and proceed to calculate their spectrum.
We then investigate the many-excitation limit numerically. We will consider
both cases of a bidirectional waveguide, $\Gamma_{L}=\Gamma_{R}=\Gamma/2$, and
chiral waveguide coupling $\Gamma_{L}=0$, $\Gamma_{R}=\Gamma$.

###  III.1 Single excitation sector

####  III.1.1 Bidirectional ordered array

Figure 2: Single excitation dispersion relation $J_{k}$ (in units of the
emission rate $\Gamma$) versus dimensionless wavevector $kd$ for the (a)
bidirectional waveguide case with $k_{0}d=0.2$, and (b) the chiral case. The
vertical dashed lines indicate the values of the wavevector where the
dispersion relation diverges.

We work in a rotating frame in order to discard the atomic frequency
$\omega_{eg}$. Then, for a bidirectional ordered array with no additional
dissipation ($\Gamma_{0}=0$), Eq. (5) reads

$H_{\rm
eff}=-\frac{i}{2}\Gamma\sum_{mn}e^{ik_{0}|x_{m}-x_{n}|}\hat{\sigma}^{n}_{eg}\hat{\sigma}^{m}_{ge}.$
(12)

For an infinite lattice of atoms ($N\rightarrow\infty$), the single-excitation
sector is diagonalized by Bloch waves of wavevector $k$, $H_{\rm
eff}|\psi_{k}\rangle=J_{k}|\psi_{k}\rangle$. The dispersion relation $J_{k}$
has previously been calculated to be Albrecht

$J_{k}=\frac{\Gamma}{4}\left(\rm cot\left[\frac{(k_{0}+k)d}{2}\right]+\rm
cot\left[\frac{(k_{0}-k)d}{2}\right]\right),$ (13)

as plotted in Fig. 2(a). The dispersion relation (13) exhibits two distinct
branches, along with a band gap near the atomic resonance frequency where no
exciations are allowed. The dispersion relation reveals the polaritonic nature
of the excitations, being close to the atomic resonance frequency ($J_{k}\sim
0$) for wavevectors significantly different than the resonant wavevector of
light ($|k|\neq k_{0}$), while strongly hybridizing with light around
$|k|\approx k_{0}$. Due to the Markov approximation, the slope of the
polariton bands approaches the dashed vertical lines $|k|=k_{0}$ (see Fig. 2)
rather than a line of slope $c$. We also note that although the effective
Hamiltonian is in principle non-Hermitian, the dispersion relation is purely
real. This reflects the fact that in an infinite waveguide with no additional
dissipation, atoms and light must form completely lossless polaritons.

Figure 3: (a) Spectrum $E_{K}$ of two-excitation eigenstates of a bi-
directional waveguide, as obtained from Eq. (18), as a function of the center
of mass momentum $K$ and for different values of the lattice constant
$k_{0}d$. The blue dots indicate continuum states, whose boundaries are
denoted by the black curves. The solid red curves indicate the bound state
dispersion relation, while the red dots indicate the special values $Kd=0$ and
$Kd=\pm\pi$ where the bound-state properties can be computed analytically. The
red dashed lines indicate the values of $K=\pm 2k_{0}$ where the bound state
dispersion branch emerges from the continuum. The horizontal dashed black
lines in the first two panels correspond to twice the frequency of the input
field used in Fig. 5. (b) Population distribution of the most localized state
as a function of the dimensionless relative coordinate $r=x_{r}/d$, for
different values of total momentum $Kd$ and lattice constant $k_{0}d$. These
numerical calculations were performed for $N=150$.

####  III.1.2 Chiral array

In the chiral case the position dependent phases can be incorporated into the
spin operators, i.e.
$\hat{\sigma}_{ge}^{n}e^{ik_{0}x_{n}}\rightarrow\hat{\sigma}_{ge}^{n}$
Pichler2015 ; ChiralRev , and Eq. (5) reads:

$H_{\rm
eff}=-\frac{i}{2}\Gamma\sum_{n}\hat{\sigma}^{n}_{eg}\hat{\sigma}^{m}_{ge}-i\Gamma\sum_{n>m}\hat{\sigma}^{n}_{eg}\hat{\sigma}^{m}_{ge}.$
(14)

The second term (involving the sum $n>m$) enforces that a given atom does not
couple to atoms to the right, and thus encodes the unidirectional propagation
of excitations along the array. An infinite lattice is again diagonalized by
Bloch waves, with dispersion relation

$J_{k}=-\frac{\Gamma}{2}\rm cot\left[\frac{kd}{2}\right].$ (15)

The dispersion relation is plotted in Fig. 2(b) and is characterized by a
positive group velocity, $v^{(1)}_{g}(k)=\Gamma
d/(4\sin^{2}{\left(\frac{kd}{2}\right)})$, for every value of $k$.
Interestingly, Eq. (15) cannot be derived by directly diagonalizing the
Hamiltonian (14) for a finite system, as $H_{\rm eff}$ is a triangular matrix
with trivial eigenvalues. Instead, one should realize that the dispersion
relation in principle should be real, and thus diagonalize the Hermitian part
of the effective Hamiltonian:

$H_{\rm
h}=-\frac{i\Gamma}{2}\left(\sum_{n>m}\hat{\sigma}^{n}_{eg}\hat{\sigma}^{m}_{ge}-\sum_{n>m}\hat{\sigma}^{m}_{eg}\hat{\sigma}^{n}_{ge}\right).$
(16)

This observation will be useful to numerically find the dispersion relation of
multi-excitation bound states.

###  III.2 Two excitation sector and bound states

Similar to the single-excitation sector, we can look for solutions of $H_{\rm
eff}|\psi^{(2)}\rangle=E|\psi^{(2)}\rangle$, where the two-excitation
eigenstate generically takes the form
$|\psi^{(2)}\rangle=\sum_{m<n}c^{(2)}_{mn}|e_{m},e_{n}\rangle$. It is
convenient to re-parametrize the coordinates in terms of center-of-mass
($x_{\rm cm}=(x_{m}+x_{n})/2$) and relative coordinates
($x_{r}=|x_{n}-x_{m}|$). We further utilize Bloch’s theorem to write the
center-of-mass wave function in terms of a wavevector $K$,

$|\psi^{(2)}\rangle=\sum_{x_{\rm cm}}e^{iKx_{\rm cm}}f(x_{r})|x_{\rm
cm}-\frac{x_{r}}{2},x_{\rm cm}+\frac{x_{r}}{2}\rangle,$ (17)

where $f(x_{r})$ is a generic function of the relative coordinate. This form
allows one to reduce the two-excitation problem to a single-excitation one
involving just the relative coordinate, and derive its spectrum as a function
of $K$.

####  III.2.1 Bidirectional ordered array

For a bidirectional array, the Hamiltonian (12) for the relative coordinate is

$H^{K}=-i\frac{\Gamma}{2}\sum_{r,r^{\prime}>0}\sum_{\epsilon,\epsilon^{\prime}=\pm
1}e^{i(k_{0}+\epsilon\frac{K}{2})|x_{r}+\epsilon^{\prime}x_{r^{\prime}}|}\hat{\sigma}^{r}_{eg}\hat{\sigma}^{r^{\prime}}_{ge},$
(18)

which depends parametrically on the center of mass momentum $K$. As pointed
out in Ref. dimermolmer2 , the two-excitation bound states can be exactly
computed, for the special values $Kd=0$ and $Kd=\pm\pi$, by using a Fourier
transform ansatz,
$|q\rangle=\sum_{r>0}e^{iqx_{r}}\hat{\sigma}^{r}_{eg}|0\rangle$. These three
dimers have energies $\omega_{0}=2\Gamma\rm cot(k_{0}d)$ and
$\omega_{\pi}=\omega_{-\pi}=2\Gamma\rm cot(2k_{0}d)$, and are characterized by
a complex relative momentum, $q_{0}d=\log{\cos{(k_{0}d)}}$ and
$q_{\pm\pi}d=-\log{\cos{(2k_{0}d)}}$. These complex momenta lead to an
exponential localization of the wave function along the relative coordinate,
which explicitly reads $p_{0}(r)\propto e^{-2x_{r}q_{0}}$ and
$p_{\pm\pi}(r)\propto\cos^{2}(\frac{Kx_{r}}{2})e^{-x_{r}q_{\pm\pi}}$. Note
that while for the dimer at $Kd=0$ the two excitations have a relative
distance of $d$ given by the TLA saturation, for $Kd=\pm\pi$ they are
separated by $2d$ due to the oscillating term in the population.

These bound states occur as one excitation creates an effective propagating
defect in the array, which can localize and drag the other excitation
dimermolmer2 . Importantly, this defect-induced localization process can occur
only in presence of an energetic gap in the relative coordinate Hamiltonian
(18). To illustrate this, we numerically diagonalize (18) as a function of the
center of mass momentum $K$, to obtain the full two-excitation spectrum. We
implement this by truncating the single-particle problem of Eq. (18) to a
large, but finite set of sites, with $1\leq r,r^{\prime}\leq N-1$ and
diagonalizing the resulting $(N-1)\times(N-1)$ matrix. The center of mass
momenta are sampled at discrete points $K=2\pi m_{K}/N$, with
$m_{K}=-N/2,...N/2-1$ for even $N$.

In Fig. 3(a) we plot the real part of the eigenvalues, ${\rm Re}(E_{K})$, for
different values of the atomic distance $k_{0}d$. As anticipated in Sec. II,
we observe an energy continuum (blue eigenvalues) that correspond to unbound
states with energies given by

$\omega(K)=J(q)+J(K-q).$ (19)

Using our analytical expression (13) for the single-particle dispersion
relation, we plot the boundaries of the allowed range of values of $\omega(K)$
in black lines. This continuum exhibits a gap of forbidden energies, which
occur in the intervals $Kd\in[2k_{0}d,\pi]$ and $Kd\in[-\pi,-2k_{0}d]$ (red
dashed lines), and at the singular point $Kd=0$. Note that the eigenvalues
surrounding $Kd=0$ are characterized by a vanishing density of states but do
not exhibit a gap. Within the gaps, we find some discrete two-excitation bound
state energies that satisfy condition (11). In particular, besides the three
exact analytic bound states previously derived at $Kd=0$ and $Kd=\pm\pi$ (red
dots), we observe the formation of well-defined dispersion branches (red
lines).

The population distribution of these states is shown in Fig. 3(b) as a
function of the relative coordinate, for different representative values of
the total momentum $Kd$ and of the atomic separation $k_{0}d$. A full
localization is indeed observed when the parameter choices coincide with a
gap. However, the neighborhood around $Kd=0$ does not exhibit a true bound
state solution, as evidenced by the large-$r$ tail in the population of the
most localized state plotted in Fig. 3(b) for $Kd=0.1$.

Such a bound state dispersion relation, as far as we know, has not been
previously derived when the photonic degree of freedom is explicitly kept. The
reduced complexity is one of the strengths of the spin model approach, which
we will soon see also extends to numerically exploring the many-body limit.

####  III.2.2 Chiral array

We now consider the chiral waveguide. This scenario is solvable by Bethe
ansatz incorporating both photons and atoms mahmo_calajo , so we just briefly
review how the same results emerge from the spin model. Using the same
methodology as in the bidirectional case, we rewrite Hamiltonian (16) in the
relative coordinate frame:

$H^{K}=-i\frac{\Gamma}{2}\sum_{r,r^{\prime}>0}\sum_{\epsilon=\pm
1}\left[e^{-i\frac{K}{2}|x_{r}+\epsilon
x_{r^{\prime}}|}-e^{i\frac{K}{2}|x_{r}+\epsilon
x_{r^{\prime}}|}\right]\hat{\sigma}^{+}_{r}\hat{\sigma}^{-}_{r^{\prime}},$
(20)

and numerically compute the full spectrum, which is plotted in Fig. 3(a). We
clearly observe a discrete dispersion branch (red curve) for two-excitation
bound states, within the gap of the unbound continuum states (blue points)
whose boundaries are given by $\omega(K)=J(k)+J(K-q)$ (black curve).

Figure 4: (a) Spectrum $E_{K}$ of two-excitation eigenstates, as obtained from
Eq. (18), as a function of the center of mass momentum $K$. The red curves
indicate the bound state dispersion relation, while the blue dots indicate
continuum states whose boundaries are given by the black curves. (b)
Population distribution of the bound state as a function of the dimensionless
relative coordinate $r=x_{r}/d$, for different values of total momentum $Kd$.
These numerical calculations were performed for $N=150$.

In contrast to the bidirectional case, here the bound state never merges with
the continuum and exists for any value of $K$. One can use a Fourier transform
ansatz for the relative coordinate wave function in order to calculate the
following exact dispersion relation

$E_{\rm BS}(K)=-2\Gamma\rm cot\left[\frac{Kd}{2}\right].$ (21)

The wave function itself is given by $p(r)\propto e^{-x_{r}\kappa}$ with
$\kappa d=-i\log{\cos{\left(\frac{Kd}{2}\right)}}$, which is plotted in Fig.
3(b) for different values of $Kd$. The exact dispersion also allows to derive
the bound state group velocity,

$v^{(2)}_{g}(K)=\frac{\partial E_{\rm BS}(K)}{\partial K}=\frac{\Gamma
d}{\sin^{2}{\left(\frac{Kd}{2}\right)}}=4v^{(1)}_{g}(k).$ (22)

In other words, the bound state travels faster than a single excitation. This
result coincides with the Bethe ansatz calculation of Ref. mahmo_calajo ,
which is nor surprising, since the polariton propagation speed is almost
entirely dictated by the atomic dispersion, rather than the speed of light
itself.

###  III.3 Self-induced transparency

We will shortly discuss bound states in the many-excitation limit, but given
the close relation of those results to self-induced transparency (SIT), we
first briefly this effect. SIT is a semi-classical phenomenon, involving the
emergence of a soliton when the atoms are treated as two-level systems and the
field classically. In order to write down the SIT equations, we move to a
continuum description, mapping the spin operators to a continuous density,
i.e. $\hat{\sigma}^{n}\rightarrow\hat{\sigma}(x)/\rho$, where $\rho$ is the
linear density of the medium and the spin operators fulfill the commutation
relation
$[\hat{\sigma}_{ge}(x),\hat{\sigma}_{eg}(x^{\prime})]=-\hat{\sigma}_{z}(x)\delta(x-x^{\prime})$
with $\hat{\sigma}_{z}(x)$ being the Z Pauli matrix. Using the Heisenberg
equations with respect to the Hamiltonian (1) for a right propagating field
(we omit the $R$ label), and taking the expectation values
$E(x,t)=\langle\hat{E}(x,t)\rangle$ and
$\sigma(x,t)=\langle\hat{\sigma}(x,t)\rangle$, we get the mean field equations

$\begin{split}&\left[\frac{\partial}{\partial t}+c\frac{\partial}{\partial
x}\right]E(x,t)=-ic\sqrt{\Gamma}\sigma_{-}(x,t)\\\ &\frac{\partial}{\partial
t}\sigma_{-}(x,t)=i\sqrt{\Gamma}\sigma_{z}(x,t)E(x,t)\\\
&\frac{\partial}{\partial t}\sigma_{z}(x,t)=4\sqrt{\Gamma}{\rm
Im}[E(x,t)\sigma^{*}_{-}(x,t)],\end{split}$ (23)

where the external dissipation has been set to zero, $\Gamma_{0}=0$. This set
of equations (23) admits a solitonic solution for the photonic field:

$E(x,t)=\frac{n_{\rm ph}\sqrt{\Gamma}}{2}\operatorname{sech}\left[\frac{n_{\rm
ph}\Gamma}{2}\left(\frac{x}{v_{g}}-t\right)\right],$ (24)

with $v_{g}=(n^{2}_{\rm ph}\Gamma c)/(n^{2}_{\rm ph}\Gamma+4c\rho)$ being the
group velocity inside the medium and $n_{\rm ph}$ being the average photon
number in teh pulse. Translated to a 1D waveguide setting for a chain of $N$
atoms with linear density $\rho=1/d$, this leads to an overall pulse delay of
mahmo_calajo

$\tau_{n_{\rm ph}}=4N/(n^{2}_{\rm ph}\Gamma).$ (25)

This solitonic solution has the property that its integrated Rabi frequency
gives a 2$\pi$ pulse (area law) McCall1 ; McCall2 ; Bullough ; mahmo_calajo

$2\sqrt{\Gamma}\int dtE(x,t)=2\pi,$ (26)

letting each individual atom undergoing a full Rabi oscillation from the
ground to excited state and back. This condition implies that photons are not
taken away from the original pulse, if the Rabi oscillation occurs on a time
scale faster than the spontaneous emission rate $\Gamma_{0}$, allowing the
soliton to propagate in a transparent fashion through the medium.

###  III.4 Many-body pulse propagation

While analytically deriving the properties of few-excitation bound states
(beyond two) seems generally challenging within the spin model, this framework
still offers a route towards numerical investigations, even in the many-body
limit. In this section we utilize such numerics to explore the dynamics of
these many-body bound states.

The finite spatial extent of the bound states generally allows for their
excitation given a coherent state input pulse. For a chiral array, it has been
proven that a coherent pulse resonant with the atomic transition can
efficiently excite $n_{\rm ph}$-excitation bound states that experience a
photon number-dependent time delay, $\tau_{n_{\rm ph}}=4N/(n_{\rm
ph}^{2}\Gamma)$. In sufficiently long systems, the combination of photon
number uncertainty of the coherent state and the number-dependent delay
results in a train of correlated photons ordered by photon number at the
output mahmo_calajo . (Note that for $n_{\rm ph}=2$, this result coincides
with the bound state group velocity derived in Eq. (22)). For shorter systems
and/or sufficiently large average photon number, the number components cannot
separate. In that case, it can be shown that a coherent state distribution of
number states in fact coincides with the classical self-induced transparency
(SIT) soliton McCall1 ; McCall2 ; Bullough ; mahmo_calajo , thus establishing
a crossover from quantum to classical nonlinear optics. Indeed, one can note
that the predicted time delay coincides with the one given in Eq. (25)
obtained by the solution of the SIT mean-field equations.

We now explore the analogue of these effects in a bi-directional array. One
can already see that a major difference compared to the chiral case is the
emergence of a single-excitation band gap around the atomic resonance, as
shown in Fig. 2(a), which causes strong reflectance of weak resonant pulses.
On the other hand, the experimental observation of SIT does not rely on
chirality, and so one might expect that SIT will arise for strong enough input
pulses in the bi-directional case. In order to investigate the change in
behavior from weak to strong pulses, we utilize a previously developed matrix
product state (MPS) algorithm for this problem (see App. A and Ref. MPSJames ;
james_ryd ; mahmo_calajo  for further details). In this representation, the
maximum bond dimension needed to obtain convergence in the simulations,
$D_{\rm max}$, is directly connected to the entanglement entropy of the system
and it acts as an important figure of merit to understand if the dynamics is
highly correlated, $D_{\rm max}\gg 1$, or if a mean field approximation for
the atoms suffices, $D_{\rm max}\sim 1$.

####  III.4.1 Exciting many-body bound states in a bidirectional array

Figure 5: Few to many-body pulse propagation in a bi-directional array of
$N=30$ atoms. In (a)-(b) we plot the time-dependent output intensity
$I_{R}(t)$ and correlation functions $G^{(n)}(t)$, given a coherent Gaussian
input pulse with width $\sigma\Gamma=3$ and average photon number
$n_{ph}=2.0$, for two different values of the atomic distance $k_{0}d$. The
input frequencies are respectively $\omega_{\rm in}=0.33$ in (a) and
$\omega_{\rm in}=-0.4\Gamma$ in (b) and are indicated by the black horizontal
dashed lines of Fig. 3(a). In (c)-(d) we plot the output intensity $I_{R}(t)$
for a solitonic input pulse (see main text) with different strength
$\tilde{n}_{ph}$ as indicated in the panels. Here we fix $k_{0}d=0.35\pi$ and
$\omega_{\rm in}=-0.4\Gamma$. In (d) the black dashed line corresponds to the
mean field solitonic ansatz. In all plots the peak of the input pulse is
indicated by the grey vertical line and is set to $\Gamma t_{0}=10$ for
(a)-(c) and $\Gamma t_{0}=0.3$ for (d). The vertical dashed line in (b) and
(c) indicate instead the expected delay of the two-photon bound states. The
simulation has been performed with an MPS based quantum trajectories algorithm
involving $N_{t}=5000$ trajectories and maximum bond dimension $D_{\rm
max}=40$.

We solve the full emitter dynamics, as governed by the master equation (4),
within the MPS formalism. We specifically start with the atoms in the ground
state $|g\rangle^{\otimes N}$ at $t=0$, and driven by a Gaussian coherent
input pulse, peaked at time $t_{0}$, with amplitude $\mathcal{E}_{\rm
in}(t)=\sqrt{n_{\rm
ph}}e^{(t-t_{0})^{2}/(2\sigma^{2})}/(\sqrt{\sigma}\pi^{1/4})$, where $n_{\rm
ph}$ is the average number of photons in the pulse. The output field operator
is given by Eq. (9), from which we can compute the transmitted intensity
$I_{R}=\langle\hat{E}_{R}^{\dagger}(t)\hat{E}_{R}(t)\rangle$ and the equal-
time higher-order correlation functions
$G^{(m)}(t)=\langle\left[\hat{E_{R}}^{\dagger}(t)\right]^{m}\left[\hat{E_{R}}(t)\right]^{m}\rangle$
($m\geq 2$) in terms of the input field amplitude and emitter correlations.

Representative output intensities and correlation functions for such a setup
are plotted in Fig. 5(a)-(b), for an array of $N=30$ atoms with lattice
constants $k_{0}d=0.2\pi$ and $k_{0}d=0.35\pi$. Considering a large array is
crucial because the propagating bound states experience a time delay at the
output, which is proportional to the number of emitters and the inverse of
their group velocity evaluated at the frequency of the input pulse
$\omega_{\textrm{in}}$, i.e. $\tau=Nd/v_{g}(\omega_{\rm in})$. Note that the
bound state energies occur at frequencies away from the atomic resonance, so
we consider input fields with a slightly detuned central frequency, i.e.
$\omega_{\textrm{in}}\neq\omega_{eg}$. For the lattice constant
$k_{0}d=0.2\pi$, we set $\omega_{\textrm{in}}\approx 0.33\Gamma$, which is
situated within the single-excitation band gap of Fig. 2(a), which suppresses
transmission of the single-photon component of the pulse. For the case of
$k_{0}d=0.35\pi$, the chosen input frequency $\omega_{\textrm{in}}=-0.4\Gamma$
instead coincides with the flat, uppermost region of the lower dispersion
branch, leading to a large delay of the transmitted single-excitation
component. For both cases we set an average photon number for the pulse to be
$n_{\textrm{ph}}=2.0$, and the pulse width and peak time to be
$\sigma\Gamma=3$ and $\Gamma t_{0}=10$, respectively.

In the case of inter-atomic distance $k_{0}d=0.2\pi$, the two-excitation bound
state dispersion (Fig. 3(a)) is extremely flat and significantly deviates from
a simple linear slope over the bandwidth of the incoming pulse. Thus, the
excited bound states exhibits both a large delay and significant distortion at
the output, which suppresses the appearance of a clear peak associated to the
bound states in the output intensity of Fig. 5(a). One can also observe a
faster exiting peak in the higher-order correlation $G^{(3)}(t)$, which could
correspond to the excitationof higher photon number bound states. A similar
peak immediately following the input is also observed in Fig. 5(b) for the
lattice constant $k_{0}d=0.35\pi$. Importantly though, here we also resolve
another intensity peak with a significant delay $\Gamma(\tau-t_{0})\sim 30$,
and an associated peak in $G^{(2)}(t)$. This peak can be clearly associated to
the excitation of a two-photon bound state resonant with the input field, as
indicated by the dashed black line in the second panel of Fig. 3(a). The
numerically calculated group velocity $v_{g}(\omega_{\textrm{in}})$ at this
frequency can be used to predict the expected delay (vertical dashed line in
Fig. 5(b)), and coincides closely with the peaks in the intensity and
$G^{(2)}$.

In summary, as with the chiral case, it is possible to observe peaks in the
output field and its correlations, which are associated with the excitation of
few-photon bound states. However, unlike the chiral case, this observation
requires more fine tuning in the bi-directional array, due to the finite
bandwidth of the bound state band and its possibly strong group velocity
dispersion.

####  III.4.2 Towards the SIT limit

In contrast, we now investigate progressively higher photon number inputs, and
the robust transition toward the mean field solution of SIT. Anticipating such
a transition, we now drive the system with an input field exhibiting an SIT-
like solitonic shape McCall1 ; McCall2 ; Bullough ; mahmo_calajo ,
$\mathcal{E}_{\rm in}(t)=(\tilde{n}_{\rm
ph}\sqrt{\Gamma_{R}}/2)\operatorname{sech}(\tilde{n}_{\rm
ph}\Gamma_{R}(t-t_{0})/2)$, where $\tilde{n}_{\rm ph}=\sqrt{n_{\rm
ph}^{2}+(\omega_{\rm in}/\Gamma)^{2}}$ is the generalized Rabi amplitude, that
matches the SIT area law given in Eq. (26). This creates a large overlap
between the field and the many-photon bound state, assuming the latter
approaches the SIT solution.

In Fig. 5(c)-(d) we plot the resulting output intensity for photon numbers
ranging from small ($\tilde{n}_{\rm ph}=2$) to large. In Fig. 5(c), for small
$\tilde{n}_{\rm ph}$, we clearly observe the appearance of the two bound state
peaks previously discussed. Not surprisingly, the second, more delayed peak
associated with the two-photon bound state gets rapidly suppressed as the
photon number is increased. For even larger photon numbers as plotted in Fig.
5(d), we see that the output intensity gradually approaches the mean field
solitonic solution

$\bar{I}_{R}(t)=\frac{\tilde{n}^{2}_{\rm
ph}\Gamma_{R}}{4}\operatorname{sech}^{2}\left[\frac{\tilde{n}_{\rm
ph}\Gamma_{R}(t-t_{0}-\tau_{\tilde{n}_{\rm ph}})}{2}\right],$ (27)

which is delayed respect to the input by an amount $\tau_{\tilde{n}_{\rm
ph}}=4N/(\tilde{n}^{2}_{\rm ph}\Gamma_{R})$ mahmo_calajo . The convergence
toward this mean field solution also indicates that the medium becomes more
transparent to the large photon number pulse, and that the impedance matching
of the pulse to the atomic array becomes irrelevant despite the bidirectional
coupling to the waveguide. Finally, we comment that, similar to the chiral
case studied in Ref. mahmo_calajo , the transition toward classical behavior
is evidenced by the decreasing maximum bond dimension required for convergence
in the MPS simulations, from $D_{\rm max}\sim 40$, for $\tilde{n}_{\rm ph}\sim
1$, to $D_{\rm max}\sim 1$, for $\tilde{n}_{\rm ph}\gg 1$.

####  III.4.3 Experimental considerations

The bound states studied here have the advantage of not relying on a chiral
waveguide setup, which was the focus of Ref. mahmo_calajo . Although schemes
to realize chiral coupling with superconducting qubits coupled to microwave
transmission lines has been proposed Guimond ; Gheeraert , it is more routine
to realize bi-directional coupling in circuit QED setups either to microwave
transmission lines or coupled resonator arrays Sundaresan ; painter1 ;
painter2 ; ustinov ; ustinov_topo ; marcoBS . Such systems are ideal to
investigate the physics predicted above, as the artificial qubits can
experience extremely low loss into channels other than waveguide emission,
$\Gamma_{0}/\Gamma<10^{-2}$, and can be located at precise positions along the
waveguide.

##  IV Photon bound states in Rydberg media

We now consider a different quantum nonlinear medium consisting of a free
space ensemble of Rydberg atoms. The level scheme of the atoms is pictorially
shown in Fig. 6(a) where, besides the $|g\rangle$ and $|e\rangle$ states,
there is an additional long-lived Rydberg level $|s\rangle$, which couples to
$|e\rangle$ via a uniform, classical control field of Rabi frequency $\Omega$
and detuning $\delta$. The propagation of photons interacting with atoms on
the $|g\rangle$-$|e\rangle$ transition can become highly nonlinear, due to
strong van der Waals interactions between Rydberg excitations. In particular,
the van der Waals interaction is responsible for shifting the resonant energy
of a Rydberg level by an amount $V(x_{r})=C_{6}/x_{r}^{6}$ given the presence
of another Rydberg excitation nearby, with $x_{r}$ being the relative distance
between the two atoms. This can significantly modify the propagation of two
photons (or more precisely, Rydberg polaritons) that are closer than a Rydberg
blockade radius $r_{b}$, to be defined later. We first introduce the effective
model used to describe the system, and then we investigate the many-body
photon dynamics, showing that in this case a generalized version of SIT can
emerge.

###  IV.1 Model and spectrum

The three-dimensional Rydberg ensemble can be approximately treated as one
dimensional provided that the blockade radius is larger than the beam waist of
the photons exciting the $|g\rangle$-$|e\rangle$ transition, such that the
photon dynamics occurs within a single transverse mode. As anticipated in Sec.
II, this situation can then be mapped to a 1D array of atoms coupled to a
chiral waveguide bienas ; MPSJames , with a large additional and independent
excited-state dissipation rate $\Gamma_{0}\gg\Gamma$ to capture the emission
of photons into $4\pi$ modes other than the mode of interest. The chirality
suppresses reflection from the array (reproducing the negligible reflection of
a free-space ensemble). We will later see how the artificial parameters of our
microscopic model (such as $N,\Gamma,\Gamma_{0}$) can be matched to physical,
macroscopic quantities in an actual ensemble, and in particular, how a
numerically tractable atom number $N$ and moderate $\Gamma_{0}$ for our system
can be used to deduce the physics in an ensemble that exhibits much larger $N$
but simultaneously larger $\Gamma_{0}$.

Within the waveguide mapping, the corresponding spin model (see Sec. II.2) in
a rotating frame is given by

$\begin{split}H_{a}=&-\frac{i\Gamma_{0}}{2}\sum_{n}\hat{\sigma}_{eg}^{n}\hat{\sigma}_{ge}^{n}-\delta\sum_{n}\hat{\sigma}^{n}_{sg}\hat{\sigma}^{n}_{gs}+\Omega\sum_{n}\left(\hat{\sigma}_{es}^{n}+\hat{\sigma}_{se}^{n}\right)\\\
&+\sum_{n<m}V(|x_{n}-x_{m}|)\hat{\sigma}^{n}_{sg}\hat{\sigma}^{n}_{gs}\hat{\sigma}^{m}_{sg}\hat{\sigma}^{m}_{gs},\end{split}$
(28)

where $\hat{\sigma}^{n}_{gs}=|g\rangle\langle s|$
($\hat{\sigma}^{n}_{sg}=(\hat{\sigma}^{n}_{gs})^{\dagger}$) and
$\hat{\sigma}^{n}_{es}=|e\rangle\langle s|$
($\hat{\sigma}^{n}_{se}=(\hat{\sigma}^{n}_{es})^{\dagger}$) are the operators
associated with the $|g\rangle$-$|s\rangle$ and $|e\rangle$-$|s\rangle$
transitions of atom $n$, respectively, and $\delta=\omega_{c}-\omega_{se}$ is
the detuning of the control driving field $\Omega$ with respect to the $e$-$s$
transition. Here, we assume that dissipation on the Rydberg level can be
neglected owing to its long lifetime.

Figure 6: (a) Sketch of the three levels of a Rydberg atom involved in the
scheme. (b) Dispersion relation associated to (28) Hamiltonian for
$\Omega=2\Gamma$ and $\delta=3\Gamma$. The color scale is associated to the
fraction of population in the Rydberg state $|s\rangle$. The horizontal black
dashed line indicates the EIT ($E^{\rm EIT}$) resonance while the red dashed
line the Stark shifted one $E_{\rm st}$. (c) Two-excitation spectrum $E_{K}$
(relative to the stark shifted resonance $E_{\rm st}$) of the effective spin
model (29), as a function of the center of mass momentum $K$ for
$C_{6}/(d)^{6}=600\bar{\Gamma}$. In grey we have superimposed the equivalent
spectrum (continuum and bound states) obtained by solving the full model with
$\delta=-10\Gamma$ and $\Omega=\Gamma$ (see App. B). The black lines indicate
the boundaries of the continuum states, obtained from the single excitation
dispersion relation.

It is instructive to first investigate the single excitation spectrum, in the
case of a far off-resonant control field $|\delta|\gg\Omega$. A representative
dispersion relation, shown in Fig. 6(b), is characterized by the occurrence of
a band that exhibits electromagnetically induced transparency (EIT), which is
centered around the EIT resonance condition $E^{\rm EIT}=-\delta$ lukinEIT
(black dashed line). Around this resonance, the polariton consists mostly of a
Rydberg excitation, as indicated by the color scale of Fig. 6(b) associated to
the amount of population on the $|s\rangle$ state. The excited state
population $|e\rangle$ and its corresponding emission is suppressed via
interference in excitation pathways, allowing this “dark state” polariton to
propagate without loss lukinEIT .

Two or more dark state polaritons in a close vicinity of each other then
strongly interact via the Rydberg interaction $V$ in (28). This significantly
affects the transparency condition and results in a strong photon-photon
interaction, which can be either dissipative Peyr , for
$|\delta|\ll\Gamma_{0}$, or dispersive in nature first , for
$|\delta|\gg\Gamma_{0}$, depending on the detuning of the control field.
Quantum nonlinear optics has been extensively explored around the EIT
transition. Notably, in the dispersive regime, this includes the observation
of signatures of bound states, via photon bunching correlations in the
outgoing field given a cw input first ; LiangBS . Theoretically, these states,
which can be of rich and varying nature, have been characterized via effective
continuum theories bienas ; Efimov ; magrebi .

Here, in order to distinguish from previous work and to also explore
connections with SIT, we will consider the nonlinear effects that arise around
the narrow, effective two-level (Stark shifted) resonance, $E_{\rm
st}\sim-\delta-\frac{\Omega^{2}}{\delta}$, indicated by the red dashed line of
Fig. 6(b). Around this resonance, the excited state $|e\rangle$ can be
effectively eliminated, and $|s\rangle$ can be considered the new, effective
“excited state” of a two-level atom, with renormalized properties. This is
evidenced by the similarity between the single-excitation dispersion relation
in the vicinity of $E_{\rm st}$ (Fig. 6(b)), and that of Fig. 2(b) for a chain
of TLA coupled to a chiral waveguide. This intuition can be made more concrete
by formally eliminating the far-detuned excited state, obtaining the effective
Hamiltonian

$\begin{split}H_{\rm eff}=&\left(E_{\rm
st}-i\frac{\bar{\Gamma}_{0}+\bar{\Gamma}}{2}\right)\sum_{n}\hat{\sigma}^{n}_{sg}\hat{\sigma}^{n}_{gs}-i\bar{\Gamma}\sum_{n>m}\hat{\sigma}^{n}_{sg}\hat{\sigma}^{m}_{gs}\\\
&+\sum_{n<m}V(|x_{n}-x_{m}|)\hat{\sigma}^{n}_{sg}\hat{\sigma}^{n}_{gs}\hat{\sigma}^{m}_{sg}\hat{\sigma}^{m}_{gs}.\end{split}$
(29)

This is indeed equivalent to the one given in Eq. (14) with renormalized
emission rates $\bar{\Gamma}=\Gamma\Omega^{2}/\delta^{2}$ and
$\bar{\Gamma}_{0}=\Gamma_{0}\Omega^{2}/\delta^{2}$. The main differences with
respect to Eq. (14) are the additional spin-spin Rydberg interaction term, and
the (large) additional independent emission $\bar{\Gamma}_{0}$.

For more than one excitation, the Rydberg interaction allows for multiple
bound state modes. This is illustrated in Fig. 6(c), where we plot the real
part of the energy $E_{K}$ vs $K$ around the Stark shifted resonance $E_{\rm
st}$. A more in-depth discussion of these Rydberg bound states is provided in
App. B and App. C. However, the large additional dissipation $\Gamma_{0}$
acting on the bound states would make their observation challenging in free-
space ensembles (although perhaps Rydberg-like interactions could be emulated
within waveguide QED Efioptica ; PRXcaneva  where the ratio of
$\Gamma/\Gamma_{0}$ can be extremely high). Instead, in the following we will
investigate the many-body limit, where it is known that SIT is robust in
normal TLA ensembles with large dissipation, and examine how SIT changes in
the presence of the long-range Rydberg interactions.

###  IV.2 Many-body dynamics

Figure 7: (a)-(b) Output intensity (a) and atomic population of an atom in the
middle of the ensemble (b) as a function of time, for pulses of different
average photon number. The form of the input field is given in Eq. (30), with
the soliton parameter $\alpha$ determined by a variational approach
illustrated in the main text. The input pulse peak is centered at
$\bar{\Gamma}t_{0}=0.05$. In both plots we have added an external spontaneous
emission rate of $\bar{\Gamma}_{0}=5\bar{\Gamma}$ and we set
$C_{6}/d^{6}=10^{5}\bar{\Gamma}$. The continuous lines represents the results
obtained numerically while the black dashed lines represent instead the
delayed solitonic ansatz. In the inset of (a) we plot the overlap between the
numerical obtained output and the solitonic ansatz as a function of number of
photons in the pulse. (c) Sketch of the effective array of superatoms arising
in the Rydberg medium. (d)-(e) Estimated number of blockade atoms (d) and
group velocity of the pulse (e) as a function of photon number for different
Rydberg interaction strengths. The points are inferred from numerics, while
the continuous lines in (d) and (e) are obtained from our effective model
described in the main text. In (e) the black line represents the standard SIT
prediction for two-level atoms without Rydberg interactions. In all plots we
take an array of $N=30$ atoms and we performed an MPS simulation of the full
master equation by setting $D_{\rm max}=70$. In all plots regarding the
soliton characterization, i.e. the inset of (a), and subfigures (d) and (e),
we set $\bar{\Gamma}_{0}=0$ consistent with the standard SIT analysis.

Finding an analytic solution of the mean field equations (23) generalized to
Rydberg interactions is more involved than in standard SIT due to the
interactions. An approximate solution for these equations has been derived in
Ref. Rydberg-SIT  using a local field approximation that neglects two-body
correlations. However, such a mean-field approach fundamentally breaks down in
the case where $V$ is sufficiently large that a blockade radius emerges, a
correlation effect in which two $|s\rangle$ excitations cannot be excited
within a distance $r_{b}$ of one another due to the prohibitive energy cost.
Here, we carry out an MPS simulation with bond dimension $D\gg 1$ to capture
the full correlated dynamics, and we show that key aspects of the solitons
that appear can be well-described by an SIT theory with renormalized
parameters.

Numerically, we drive the system with the input field

$\mathcal{E}_{\rm in}(t)=\frac{n_{\rm
ph}\alpha\sqrt{\bar{\Gamma}}}{2}\operatorname{sech}\left(\frac{\bar{\Gamma}}{2}\alpha^{2}n_{\rm
ph}(t-t_{0})\right),$ (30)

which represents a generalization of the SIT solution (recovered for
$\alpha=1$) and our ansatz of what a soliton in the presence of Rydberg
interactions might look like. Note that when acting on a single atom, this
pulse violates the usual area law for the integrated Rabi frequency:

$2\sqrt{\bar{\Gamma}}\int dt\mathcal{E}_{\rm in}(t)=\frac{2\pi}{\alpha}.$ (31)

The “best fit” parameter $\alpha$ will in general be a function of the number
of photons and the Rydberg interaction, i.e. $\alpha:=\alpha(n_{\rm
ph},C_{6})$. We determine the optimal value $\alpha_{\textrm{opt}}$ with a
variational approach that maximizes the overlap $O=\int dt\,{\rm
min}\\{I_{R}(t),I_{\rm sol}(t)\\}/\sqrt{\int dtI_{R}(t)\cdot\int dtI_{\rm
sol}(t)}$ between the numeric and the expected output field intensity,
$I_{R}(t)$ and $I_{\rm sol}(t)=|\mathcal{E}_{\rm in}(t-\tau)|^{2}$
respectively, where the pulse delay $\tau=N/v_{g}$ and corresponding group
velocity are numerically determined based on the position of the output
intensity peak.

The inset of Fig. 7(a) shows, for a chain of $N=30$ atoms, a progressively
improving agreement between the output intensity and the expected soliton
solution, when the number of photons is increased, with the optimal overlap
reaching $O(\alpha_{\rm opt})\sim 99\%$ for $n_{\rm ph}=200$. For the
determined optimal values $\alpha_{\textrm{opt}}$, we plot in Fig. 7(a) the
corresponding output intensity at different pulse strengths. We also compare
the numerically obtained output intensity with the expected damped solitonic
solution $I_{\rm sol}(t)e^{-\bar{\Gamma}_{0}\tau}$ (black dashed line). This
approximately accounts for the external dissipation
$\bar{\Gamma}_{0}=5\bar{\Gamma}$ introduced in the numerical simulations, and
the two sets of curves are observed to agree well with one another.

The difference compared to usual SIT is made evident in Fig. 7(b), where we
plot the time-dependent atomic excited population, $p_{s}(t)$, of an atom in
the middle of the chain. We observe that, even in regimes of strong pulses,
the excited population does not go to unity and back to zero as expected of a
$2\pi$ pulse, confirming the violation of the area law.

We now present an effective theory for the observed phenomena. We begin by
noting that, if $V$ is strong enough and covers a sufficiently large range,
some portion of the atomic medium can be blockaded. In particular, given a
Rydberg excitation at the origin and the peak Rabi frequency
$\Omega_{r}=n_{\rm ph}\alpha_{\textrm{opt}}\bar{\Gamma}$ associated with the
pulse (30), the condition $V(r_{b})=\Omega_{r}$ sets the distance at which the
pulse can overcome the dispersive Rydberg energy shift and efficiently drive a
second Rydberg excitation. The resulting “blockade radius” is then
$r_{b}=\left(C_{6}/\Omega_{r}\right)^{1/6}$.

The blockade region consisting of $N_{b}=2r_{b}$ atoms acts like an effective,
single two-level “super-atom” Hofferberth ; Hofferberth2 , within which only a
single collective Rydberg excitation can be generated. Furthermore, this
excitation experiences a collectively enhanced emission rate $\Gamma_{\rm
s}=N_{b}\bar{\Gamma}$ into the probe mode, while the dissipation
$\bar{\Gamma}_{0}$ into other modes remains fixed. This super-atom picture
suggests that a generalized version of SIT should emerge, where the parameter
in Eq. (30) is set to $\alpha_{\textrm{opt}}=\sqrt{N_{b}}$ to account for the
enhanced coupling, and in turn, Eq. (31) for this choice of $\alpha$ states
that each super-atom experiences a full $2\pi$-pulse. The whole Rydberg medium
can then be described as an array of $N/N_{b}$ super-atoms, pictorially shown
in Fig. 7(c), and the standard SIT time delay in Eq. (25) is modified to
become $\tau=\frac{N}{N_{b}}\frac{4}{\bar{\Gamma}n^{2}_{\rm ph}N_{b}}$. This
equation can be recast in terms of a more experimentally relevant quantity,
the optical depth per blockade radius
$D_{b}=2N_{b}\bar{\Gamma}/\bar{\Gamma}_{0}$ MPSJames ; james_ryd , such that

$\bar{\Gamma}_{0}\tau=\frac{N}{N_{b}}\frac{8}{D_{b}n^{2}_{\rm ph}}.$ (32)

Note that $D_{b}$, $N/N_{b}$, the photon number in the pulse, and the free-
space spontaneous emission rate $\bar{\Gamma}_{0}$ of the effective two-level
transition all have well-defined meaning in an experimental setup of a Rydberg
ensemble, which thus provides a connection between our microscopic spin model
and a physical system.

In order to gauge the accuracy of this effective description, in Fig. 7(d) we
plot in solid lines the expected number of atoms within a blockade radius
$N_{b}$ versus photon number $n_{\rm ph}$ assuming our effective theory is
correct, by solving the equation $r_{b}=N_{b}/2=\left(C_{6}/n_{\rm
ph}\bar{\Gamma}\sqrt{N_{b}}\right)^{1/6}$. This is plotted for several
different values of $C_{6}$, as indicated by different colors. On the other
hand, from the numerically determined value for $\alpha_{\textrm{opt}}$, we
can use $N_{b}=\alpha_{\textrm{opt}}^{2}$ to infer $N_{b}$ (points). It is
seen that the two approaches agree well with each other. Likewise, in Fig.
7(e), we plot with points the group velocity, as given from the numerically
determined delay by $v_{g}=N/\tau$. We then plot the same quantity
$v_{g}=N/\tau$ in solid lines, where $\tau$ is given by the effective
description (32), and where the expected $N_{b}$ is taken from the solid lines
of Fig. 7(d). Again, good agreement is seen. Finally, for comparison, we plot
the the expected group velocity for normal SIT (without Rydberg interactions)
in black, which is seen to be slower than in the blockaded case.

We believe that the discrepancies between the numerical results and the
effective description in Fig. 7(d)-(e) can be mainly attributed to three
sources of errors. The first is the intrinsic discreteness of our numerical
model. For example, the moderate number of atoms $N_{b}\lesssim 10$ within a
blockade region in our simulations suggests that moving to a continuum
description should not be entirely accurate. A second error comes from the
numerical imperfections in the evaluation of the optimal value
$\alpha_{\textrm{opt}}$. Specifically, due to the large simulation complexity,
the optimal value was obtained using a limited number of sampling points
($\sim 20$) for $\alpha$. An additional (non-numerical) discrepancy is that,
similar to the case of SIT with normal TLA, the approach from quantum many-
photon bound states to semi-classical SIT is a gradual one as a function of
increasing photon number. For the lower range of photon numbers used in our
simulations, it is possible that additional quantum many-body features are
present, which would be interesting to pinpoint and explore further in future
work.

We now consider the experimental feasibility of observing this many-body
dynamics in a Rydberg ensemble. State-of-the-art experiments LiangBS ;
Hofferberth ; Hofferberth2  allow for a large optical depth per blockade
radius, $D_{b}\gtrsim 1$, and multiple blockade regions in the ensemble,
$N/N_{b}\sim 10$. For moderate photon numbers $n_{\rm ph}\sim 10^{2}$ in the
pulse, it is then possible to acquire the time delay before full absorption
takes place, $\bar{\Gamma}_{0}\tau\ll 1$. Separately, from Eqs. (30) and (32),
the ratio of the delay to the temporal width of the input pulse is given by
$\tau/t_{\rm in}\sim(2/n_{\rm ph})(N/N_{b})$. For parameters like above, this
allows for the delay to be a reasonable, detectable fraction of the input
pulse width. Finally, as pointed out in Ref. Rydberg-SIT , this Rydberg SIT
dynamics could also be probed without making use of an intermediate
$|e\rangle$ state, by directly driving a ground-Rydberg state transition.

##  V Conclusions

In summary, we have presented a unified method to treat photon bound states
within an atomic nonlinear medium, based upon a spin model formulation. This
description allows one to promptly identify the emerging correlated states in
different scenarios, such as an array of quantum emitters coupled to an
optical waveguide and an ensemble of Rydberg atoms. In both cases, the
formalism allows one to obtain the two-excitation bound state dispersion
relation exploiting a convenient description in the relative coordinate frame,
understand the effect of this dispersion relation on the correlation functions
of outgoing fields given an input pulse With this knowledge, and investigate
the transition from few- to many-photon behavior. In the many-body case, we
show how SIT or generalizations thereof emerge in all the systems studied.

While the two-excitation limit and many-excitation, semi-classical SIT limit
are possible to treat in semi-analytic fashion, in future work it would be
interesting to develop techniques to better understand the intermediate case,
where quantum effects in the many-photon pulse might still persist in the
output. For example, it might be interesting to see whether the spin model
might be more amenable to field theoretical techniques. Alternatively, it
might be feasible to use time-independent MPS techniques to target and
investigate the multi-excitation bound eigenstates themselves, in order to
better understand their nature within the atomic medium.

## Acknowledgments

The authors thank Sahand Mahmoodian, Anders S. Sørensen, Angelo Carollo,
Francesco Ciccarello and Cosimo Rusconi for valuable discussions. G.C.
acknowledge that results incorporated in this standard have received funding
from the European Union Horizon 2020 research and innovation programme under
the Marie Sklodowska-Curie grant agreement No 882536 for the project QUANLUX.
D.E.C. acknowledges support from the European Unions Horizon 2020 research and
innovation programme, under European Research Council grant agreement No
101002107 (NEWSPIN); the Government of Spain (Europa Excelencia program
EUR2020-112155, Severo Ochoa program CEX2019-000910-S, and MICINN Plan
Nacional Grant PGC2018-096844-B-I00), Generalitat de Catalunya through the
CERCA program, AGAUR Project No. 2017-SGR-1334, Fundació Privada Cellex,
Fundació Mir-Puig, and Secretaria d’Universitats i Recerca del Departament
d’Empresa i Coneixement de la Generalitat de Catalunya, co-funded by the
European Union Regional Development Fund within the ERDF Operational Program
of Catalunya (project QuantumCat, ref. 001-P-001644).

##  Appendix A Matrix product state simulation

Figure 8: Spectrum $E_{K}$ of two-excitation eigenstates for the three-level
atomic scheme presented in Sec. IV.1, as a function of the center of mass
momentum $K$ and for different detunings $\delta$ of the control field. The
continuous yellow line in (a)-(d) indicates the bound state obtained for the
three level atom, which differs from the expected TLA bound state shown by the
yellow dashed line in (a). In (e)-(f) the red continuous line indicates the
bound state that emerges around the effective Stark-shifted TLA transition, as
discussed in the main text. Its dispersion indeed approaches the one obtained
with the rescaled effective model given in Eq. (41) (red dashed line in (f)).
The dashed red lines in (a)-(d) indicate the region where the bound states
associated with the effective Stark-shifted TLA transition eventually arise
for large $|\delta|$. In all plots we choose $C_{6}/(d)^{6}=0$ and
$\Omega=\Gamma$.

Solving the full system dynamics described by the master equation (4) for
large atomic arrays, $N\gtrsim 20$, is extremely challenging with brute force
numerical techniques. To overcome this problem we adopt an MPS representation
for the atomic array, which allows us to highly reduce the degrees of freedom
needed to efficiently simulate the system. In particular, we make use of a
quantum trajectories algorithm where the state of the system evolves under the
effective Hamiltonian (5), along with stochastic quantum jumps as described in
the details of Ref. MPSJames . The main idea of the MPS representation
consists in reshaping a generic quantum state
$|\phi\rangle=\sum_{i_{1},..i_{N}}\psi_{i_{1},i_{2},..i_{N}}|i_{1},i_{2},..i_{N}\rangle$
into a matrix product state of the form:

$|\phi\rangle=\sum_{i_{1},..i_{N}}A_{i_{1}}A_{i_{2}}...A_{i_{N}}|i_{1},i_{2},..i_{N}\rangle$
(33)

where, for each specific set of physical indices $\\{i_{1},i_{2},..i_{N}\\}$,
the product of the $A_{i_{j}}$ matrices gives back the state coefficient
$\psi_{i_{1},i_{2},..i_{N}}$. Each matrix $A_{i_{j}}$ has dimension
$D_{j-1}\times D_{j}$ and finite-edge boundary conditions are assumed by
imposing $D_{1}=1$ and $D_{N}=1$. The bond dimension $D_{j}$ is a crucial
parameter because is directly connected to the entanglement entropy of the
system. This implies that problems characterized by a slow entanglement
entropy growth, as the ones considered in this paper, can be efficiently
described by MPS ansatz with small bond dimension Verstraete ; Schollw . To
compute the time evolution of the system we derive a matrix product operator
(MPO) representation for the effective Hamiltonian and jump operators as
described in Ref. MPSJames . The evolution of the initial ground state is then
computed by using a Runge-Kutta method.

One difficulty for the simulation arises in the case of Rydberg atoms because
the Rydberg interaction term,
$\sum_{n<m}V(|x_{n}-x_{m}|)\hat{\sigma}^{n}_{sg}\hat{\sigma}^{n}_{gs}\hat{\sigma}^{m}_{sg}\hat{\sigma}^{m}_{gs}$
with $V(|x_{n}-x_{m}|)=C_{6}/|x_{n}-x_{m}|$, does not have an exact MPO
representation due to its power law nature. To overcome this issue we
approximate the power law potential as a series of exponentials
$V(x_{r})\approx\sum_{n}\alpha_{n}\lambda_{n}^{x_{r}}$, as described in Refs.
Pirvu ; Crosswhite ; fro  and explicitly developed for atom-waveguide
interactions in Ref. james_ryd . It is sufficient to truncate the maximum
strength of the Rydberg interaction at the radius where the potential assumes
a value an order of magnitude bigger than other characteristic scales, given
by $\bar{\Gamma}$ and $\sqrt{N_{b}}n_{\rm ph}\bar{\Gamma}$ for the few and
many body regimes respectively.

In addition to the Hamiltonian term, an MPO representation can be derived also
for the output field and its associated correlation operators. Both the time
evolution and the computation of the observables at each time step are
evaluated by applying an MPO to an MPS. This operation progressively increases
the MPS bond dimension but, in order to keep the computation efficient, the
bond dimension $D_{j}\leq D_{\rm max}$ can be truncated after each step.

In section IV.2, we employ a slightly different technique, by directly
representing the density matrix in MPS form and solving the master equation
without making use of the quantum trajectories algorithm, as explained in Ref.
james_ryd . This choice, exploited for the plots of Fig. 7, is well motivated
by the relatively low maximum bond dimension required for the density matrix.

##  Appendix B Spin model approach for the Rydberg media

In this section we use the spin model to derive the single- and two-excitation
spectrum for the three-level Rydberg atom scheme presented in Sec. IV.1.

###  B.1 Single excitation and EIT

For a single excitation the Rydberg interaction term given in (28) does not
play a role and it is convenient to rewrite the spin Hamiltonian (5) in
wavevector space $k$ where it reads:

$H_{\rm
eff}=\tilde{J}_{k}\sum_{k}\hat{\sigma}_{eg}^{k}\hat{\sigma}_{ge}^{k}-\delta\sum_{k}\hat{\sigma}^{k}_{se}\hat{\sigma}^{k}_{es}+\Omega\sum_{k}\left(\hat{\sigma}_{es}^{k}+\hat{\sigma}_{se}^{k}\right),$
(34)

where $\tilde{J}_{k}=-\frac{\Gamma}{2}\cot{(kd/2)}-i\Gamma_{0}/2$ coincides,
except for the spontaneous emission term $\sim\Gamma_{0}$, with the dispersion
for a chiral atomic array derived in Eq. (15). Hamiltonian (34) can be exactly
diagonalized by considering a linear combination of Rydberg and excited state
excitations:
$|\psi_{k}\rangle=\alpha_{k}|e_{k}\rangle+\beta_{k}|s_{k}\rangle$. This leads
to the following coupled equations:

$\begin{split}&\left(E_{k}-\tilde{J}_{k}\right)\alpha_{k}=\Omega\beta_{k}\\\
&\left(\delta+E_{k}\right)\beta_{k}=\Omega\alpha_{k}.\end{split}$ (35)

These equations provide the upper (U) and lower (L) dispersion branches

$E^{(U/L)}_{k}=\frac{\tilde{J}_{k}-\delta}{2}\pm\frac{1}{2}\sqrt{(\tilde{J}_{k}-\delta)^{2}+4(\Omega^{2}+\delta\tilde{J}_{k})},$
(36)

associated to the states:

$|\psi^{(U/L)}_{k}\rangle=\Psi^{\dagger(U/L)}_{k}|0\rangle=\frac{(\delta+E^{(U/L)}_{k})\hat{\sigma}^{k}_{ge}+\Omega\hat{\sigma}^{k}_{gr}}{\sqrt{|\delta+E^{(U/L)}_{k}|^{2}+\Omega^{2}}}|0\rangle.$
(37)

The combination of these two branches gives the dispersion relation shown in
Fig. 6(b), which is characterized by the occurrence of an EIT transmission
band centered at the EIT resonance $E^{\rm EIT}=-\delta$ (dashed horizontal
line). Around this resonance the dispersion (36) reduces to

$E_{k}\simeq-\delta+\frac{\Omega^{2}}{\Gamma}kd+\frac{\Omega^{2}}{\Gamma^{2}}(\delta-i\Gamma_{0}/2)(kd)^{2}+...,$
(38)

and the state (37) becomes mainly a Rydberg-like excitation with
$\alpha_{k}\sim 0$. The modes occurring around this resonance do not suffer
from spontaneous emission into free space ($\Gamma_{0}$) due to the negligible
excited state population, and these “dark state” polaritons allow the medium
to be transparent to light in this frequency range lukinEIT . These dark state
polaritons propagate in the medium with EIT group velocity $v_{g}=\partial
E_{k}/\partial k|_{k=0}=\Omega^{2}d/\Gamma$ and effective mass $m_{\rm
eff}=\hbar/(2\partial E^{2}_{k}/\partial
k^{2}|_{k=0})=\hbar\Omega^{2}/(2v_{g}^{2}\delta)$.

###  B.2 Two excitation spectrum

The two-excitation subspace is spanned by the basis set
$\\{|e_{n}e_{m}\rangle,|s_{n}s_{m}\rangle,|e_{n}s_{m}\rangle,|e_{m}s_{n}\rangle\\}$
with $n>m$. Similarly as done for the two-level atom array in Sec. III.2, we
can re-parametrize the eigenstates in the center of mass, $x_{\rm
cm}=(x_{m}+x_{n})/2$, and relative coordinate, $x_{r}=|x_{n}-x_{m}|$, assuming
a plane wave ansatz along $x_{\rm cm}$ of the form

$\begin{split}&|\psi^{(2)}\rangle=\sum_{x_{\rm cm}}e^{iKx_{\rm
cm}}\left(f_{1}(x_{r})\hat{\sigma}^{(x_{\rm
cm}-x_{r}/2)}_{eg}\hat{\sigma}^{(x_{\rm cm}+x_{r}/2)}_{eg}\right.\\\
&\left.+f_{2}(x_{r})\hat{\sigma}^{(x_{\rm
cm}-x_{r}/2)}_{sg}\hat{\sigma}^{(x_{\rm cm}+x_{r}/2)}_{sg}\right.\\\
&\left.+f_{3}(x_{r})\hat{\sigma}^{(x_{\rm
cm}-x_{r}/2)}_{eg}\hat{\sigma}^{(x_{\rm cm}+x_{r}/2)}_{sg}\right.\\\
&\left.f_{4}(x_{r})\hat{\sigma}^{(x_{\rm
cm}-x_{r}/2)}_{sg}\hat{\sigma}^{(x_{\rm
cm}+x_{r}/2)}_{eg}\right)|0\rangle,\end{split}$ (39)

where $f_{l}(x_{r})$ are generic functions of the relative coordinate. The
effective single-excitation problem in the relative coordinate is more
complicated than the case of TLA, as the excitation can occupy one of four
distinct sectors characterized by the creation operators
$\hat{S}^{\dagger}_{ee}$, $\hat{S}^{\dagger}_{ss}$, $\hat{S}^{\dagger}_{se}$
and $\hat{S}^{\dagger}_{es}$. With this notation the Hamiltonian for the
relative coordinate reduces to

$\begin{split}&H^{K}=-\sum_{r>0}\left[\left(2\delta-V(x_{r})\right)\hat{S}^{\dagger
r}_{ss}\hat{S}^{r}_{ss}+\delta\hat{S}^{\dagger
r}_{se}\hat{S}^{r}_{se}+\delta\hat{S}^{\dagger
r}_{es}\hat{S}^{r}_{es}\right]\\\
&-i\frac{\Gamma}{2}\sum_{r,r^{\prime}>0}\sum_{\epsilon=\pm
1}\left[e^{-i\frac{K}{2}|x_{r}+\epsilon
x_{r^{\prime}}|}-e^{i\frac{K}{2}|x_{r}+\epsilon
x_{r^{\prime}}|}\right]\hat{S}^{\dagger r}_{ee}\hat{S}^{r^{\prime}}_{ee}\\\
&-i\frac{\Gamma}{2}\sum_{r>r^{\prime}}\left[e^{-i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{es}\hat{S}^{r^{\prime}}_{es}-e^{i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{se}\hat{S}^{r^{\prime}}_{se}\right]\\\
&-i\frac{\Gamma}{2}\sum_{r<r^{\prime}}\left[e^{-i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{se}\hat{S}^{r^{\prime}}_{se}-e^{i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{es}\hat{S}^{r^{\prime}}_{es}\right]\\\
&-i\frac{\Gamma}{2}\sum_{r,r^{\prime}>0}\left[e^{-i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{es}\hat{S}^{r^{\prime}}_{se}-e^{i\frac{K}{2}|x_{r}-x_{r^{\prime}}|}\hat{S}^{\dagger
r}_{se}\hat{S}^{r^{\prime}}_{es}\right]\\\
&+\Omega\sum_{r>0}\left[\hat{S}^{\dagger
r}_{es}\hat{S}^{r^{\prime}}_{ee}+\hat{S}^{\dagger
r}_{se}\hat{S}^{r^{\prime}}_{ee}+\hat{S}^{\dagger
r}_{es}\hat{S}^{r^{\prime}}_{ss}+\hat{S}^{\dagger
r}_{se}\hat{S}^{r^{\prime}}_{ss}+\rm H.c.\right],\end{split}$ (40)

which depends parametrically on the center-of-mass momentum $K$. Eq. (40) can
be diagonalized numerically.

####  B.2.1 Two excitation spectrum for $C_{6}=0$

In Fig. 6(c) of the main text we presented the two-excitation spectrum of the
Rydberg medium spectrum, obtained by diagonalizing the effective two-level
model (29). Here we want to show how this spectrum arises, within the full
three-level model, starting from the resonant case, $\delta=0$ and
progressively moving off resonance $\delta\gg\Gamma$. In Fig. 8 we plot the
spectrum obtained by the diagonalization of Eq. (40) for the paradigmatic case
where there is no Rydberg interaction, $C_{6}=0$. In Fig. 8 (a) for $\delta=0$
we observe a bound state in the gap (continuous yellow line) with a dispersion
that differs from the one of the TLA case (dashed yellow line) due to the
three level nature of the medium. In the main text we discussed the off
resonance regime where an effective (rescaled) TLA atom description can be
used. The transition towards this regime is shown in Fig. 8(b)-(d), where the
density of states inside the region highlighted by the red dashed lines
decreases for increasing $\delta$. Moving further off-resonance (panels
(e)-(f)), we observe that a new gap arises in this region, which is captured
by the effective TLA Hamiltonian (29). This effective model, in absence of
Rydberg interactions, owns a bound state solution (red dashed line in Fig. 8
(f)) that follows the dispersion:

$\omega(K)=2\frac{\Omega^{2}}{\Delta_{p}}-2\Gamma\left(\frac{\Omega}{\Delta_{p}}\right)^{2}\cot(Kd/2)$
(41)

and well approximates the one obtained by the full diagonalization of Eq. (40)
(continuous red line).

Figure 9: (a)-(c) Spectrum $E_{K}$ of two-excitation eigenstates versus
center-of-mass momentum $K$ for $\delta=-8\Gamma$, $\Omega=\Gamma$ and
different strengths of the Rydberg interaction: $C_{6}/d^{6}=0.1\Gamma$ in
(a), $C_{6}/d^{6}=10\Gamma$ in (b) and $C_{6}/d^{6}=1000\Gamma$ in (c). Here,
we focus on the region of the spectrum near the effective, Stark-shifted TLA
transition. The grey line in (a) indicates the bound state dispersion relation
in absence of Rydberg interaction due to the TLA nonlinearity, while the
Rydberg interaction induced bound states are given in red. The black lines
indicate the boundaries of the continuum states. (d)-(f) Rydberg population
distribution as a function of the relative distance index $r=x_{r}/d$, for the
corresponding strengths of the Rydberg interaction shown in (a)-(c). The bound
states populations plotted in the panels correspond to specific values of $K$
and are highlighted in (a)-(c) by the black dots.

####  B.2.2 Two excitation spectrum for $C_{6}\neq 0$

We now investigate the spectrum of the two-excitation subspace in presence of
Rydberg interactions. Different strengths of the Rydberg interaction affect
the bound state dispersion relation as illustrated in Fig. 9. In particular,
in Fig. 9(a), we see how a weak interaction is sufficient to drastically
change the bound state dispersion (red) compared to the one previously
obtained for an array of two-level atoms (grey). Increasing the interaction
strength $C_{6}$, multiple bound states progressively arise with a dispersion
that becomes steeper as they move further away from the continuum of states.
The increasing interaction strength also affects the population distribution
in the relative coordinate. In particular, in Figs. 9(d)-(f) one sees that the
distribution $P_{rr}(r)$ of two Rydberg excitations exhibits a progressive
growth in the separation between them. This localization distance is
approximately given by the Rydberg blockade radius,
$r_{b}=(C_{6}/\bar{\Gamma})^{1/6}$, defined in the main text for the frequency
regime close to the Stark-shifted resonance.

##  Appendix C Exciting multiple bound states in the Rydberg media

Figure 10: (a) Same eigenvalue spectrum as the one shown in Fig. 6(c) of the
main text. (b) Output intensity $I_{R}(t)$ and two-photon correlation function
$G^{(2)}(t)$, due to a coherent Gaussian input pulse propagating through an
array of $N=40$ atoms with $\sigma\bar{\Gamma}=2$ and $n_{ph}=1.0$. The
expected time delays arising from excitation and propagation of bound states,
$\tau_{l}=Nd/v^{l}_{g}(\omega_{in})$, are indicated by vertical dashed lines
and are computed numerically from the group velocity of the bound state
branches $l$ plotted in (a). The grey vertical line instead indicates the
initial time of the input pulse set to $\bar{\Gamma}t_{0}=10$. The simulation
has been performed with an MPS based quantum trajectories algorithm involving
$N_{t}=5000$ trajectories and maximum bond dimension $D_{\rm max}=40$.

Here, we show how the multiple two-excitation bound states occurring in the
Rydberg system can be excited given an input pulse. To simulate the dynamics
we make use of the effective TLA model (29) and employ the MPS simulation as
previously discussed. The bound states shown in the spectrum of Fig. 6(c)
(repeated in Fig. 10(a) for convenience) can be excited by sending a coherent
Gaussian input pulse at the frequency $E_{\rm st}$ through the atomic array.
In the ideal lossless regime, $\Gamma_{0}=0$, the time-dependent output
intensity $I_{R}(t)$ and the equal time two-photon correlation function
$G^{(2)}(t)$ exhibit isolated peaks that can be associated to the different
bound states of Fig. 10(a) through their acquired time delay
$\tau_{l}=Nd/v^{l}_{g}(\omega_{in})$, where $l$ is the label of the bound
state branches (circled numbers in Fig. 10). In order to resolve the multiple
bound states in the output field, we have considered in Fig. 10(b) a large
array of $N=40$ atoms. This allows that the delay between two different peaks
is bigger than the sum of their widths, a condition that can be recast in the
form $Nd\gg r_{b}(v_{g}^{l}+v_{g}^{m})/|v_{g}^{l}-v_{g}^{m}|$.

In a large array, moderate values of external spontaneous emission
$\Gamma_{0}/\Gamma\sim 0.1$ can already damp and smear out the features in the
output field. As discussed in Sec. IV.2 in the main text, the effect of
external dissipation decreases by going toward the many-body regime, where the
stimulated emission rate into the waveguide, $\sim\Gamma n_{\rm ph}$, can
overcome the free space emission for large photon number input pulses. On the
other hand, this few-excitation dynamics might be observable in circuit QED
platforms where Rydberg-like interactions are implemented in arrays of
artificial atoms coupled to a transmission line. It would be interesting to
further investigate how different types of correlated light might be realized
through the engineering of direct qubit-qubit interactions.

## References

  * (1) G. Agrawal, Nonlinear Fiber Optics, 5th edn Academic (2012).
  * (2) Y. Lai and H. A. Haus, Quantum Theory of Solitons in Optical Fibers. I. Time-Dependent Hartree Approximation, Phys. Rev. A 40, 844 (1989).
  * (3) Y. Lai and H. A. Haus, Quantum Theory of Solitons in Optical Fibers. II. Exact Solution, Phys. Rev. A 40, 854 (1989).
  * (4) O. Firstenberg, C. S. Adams, and S. Hofferberth, Nonlinear quantum optics mediated by Rydberg interactions, Phys. B: At. Mol. Opt. Phys. 49, 152003 (2016).
  * (5) C. Murray and T. Pohl, Quantum and Nonlinear Optics in Strongly Interacting Atomic Ensembles, Adv. in Atomic, Molecular, and Optical Physics 65, 321-372 (2016).
  * (6) O. Firstenberg, T. Peyronel, Q-Y. Liang, A. V. Gorshkov, M. D. Lukin, and V. Vuletic, Attractive photons in a quantum nonlinear medium, Nature 502, 71-75 (2013).
  * (7) Q-Y. Liang, A. V. Venkatramani, S. H. Cantu, T. L. Nicholson, M. J. Gullans, A. V. Gorshkov, J. D. Thompson, C. Chin, M. D. Lukin, and V. Vuletic, Observation of three-photon bound states in a quantum nonlinear medium, Science 359, 783 (2018).
  * (8) D. Roy, C. M. Wilson, and O. Firstenberg, Colloquium: Strongly interacting photons in one-dimensional continuum, Rev. Mod. Phys. 89, 021001 (2017).
  * (9) P. Lodahl, S. Mahmoodian, and S. Stobbe, Interfacing single photons and single quantum dots with photonic nanostructures, Rev. Mod. Phys. 87, 347 (2015).
  * (10) D. E. Chang, J. S. Douglas, A. González-Tudela, C. L. Hung, and H. J. Kimble, Colloquium: Quantum matter built from nanoscopic lattices of atoms and photons, Rev. Mod. Phys. 90, 031002 (2018).
  * (11) J.-T. Shen and S. Fan, Strongly Correlated Two-Photon Transport in One-Dimensional Waveguide Coupled to A Two-Level System, Phys. Rev. Lett. 98, 153003 (2007).
  * (12) J.-T. Shen and S. Fan, Strongly correlated multiparticle transport in one dimension through a quantum impurity, Phys. Rev. A 76, 062709 (2007).
  * (13) H. Zheng, D. Gauthier, and H.U. Baranger, Strongly-Correlated Photons Generated by Coupling a Three-or Four-level System to a Waveguide, Phys. Rev. A 85, 043832 (2012).
  * (14) A. S. Prasad, J. Hinney, S. Mahmoodian, K. Hammerer, S. Rind, P. Schneeweiss, A. S. Sörensen, J. Volz, and A. Rauschenbeutel, Correlating photons using the collective nonlinear response of atoms weakly coupled to an optical mode, Nat. Photonics 14, 719 (2020).
  * (15) S. Mahmoodian, G. Calajó, D. E. Chang, K. Hammerer, and A. S. Sörensen, Dynamics of many-body photon bound states in chiral waveguide QED, Phys. Rev. X 10, 031011 (2020).
  * (16) H. Zheng, D. Gauthier, and H.U. Baranger, Cavity-Free Photon Blockade Induced by Many-Body Bound States, Phys. Rev. Lett. 107, 223601 (2011).
  * (17) D. E. Chang, V. Vuletic, and M. D. Lukin, Quantum nonlinear optics - photon by photon, Nat. Photonics 8, 685 (2015).
  * (18) H. A. Bethe, Strongly correlated multiparticle transport in one dimension through a quantum impurity, Z. Phys. 71, 205 (1931).
  * (19) S. Mahmoodian, M. Cepulkovskis, S. Das, P.Lodahl, K. Hammerer, and A. S. Sörensen, Strongly Correlated Photon Transport in Waveguide Quantum Electrodynamics with Weakly Coupled Emitters, Phys. Rev. Lett. 121, 143601 (2018).
  * (20) P. Bienias, S. Choi, O. Firstenberg, M. F. Maghrebi, M. Gullans, M. D. Lukin, A. V. Gorshkov, and H. P. Büchler Scattering resonances and bound states for strongly interacting Rydberg polaritons, Phys. Rev. A 90, 053804 (2014).
  * (21) M. J. Gullans, S. Diehl, S. T. Rittenhouse, B. P. Ruzic, J. P. D Incao, P. Julienne, A. V. Gorshkov, and J. M. Taylor, Efimov States of Strongly Interacting Photons, Phys. Rev. Lett. 119, 233601 (2017).
  * (22) M. F. Maghrebi, M. J. Gullans, P. Bienias, S. Choi, I. Martin, O. Firstenberg, and M. D. Lukin, H. P. Büchler, and A. V. Gorshkov, Coulomb Bound States of Strongly Interacting Photons, Phys. Rev. Lett. 115, 123601 (2015).
  * (23) A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, Atom-light interactions in quasi-1D nanostructures: a Green s function perspective, Phys. Rev. A 95, 033818 (2017).
  * (24) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, D. E. Chang, Exponential Improvement in Photon Storage Fidelities Using Subradiance and Selective Radiance in Atomic Arrays, Phys. Rev. X 7, 031024 (2017).
  * (25) S. J. Masson and A. Asenjo-Garcia, Atomic-waveguide quantum electrodynamics, Phys. Rev. Research 2, 043213 (2020).
  * (26) A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, Waveguide quantum electrodynamics: collective radiance and photon-photon correlations, arXiv:2103.06824 (2021).
  * (27) M. T. Manzoni, D. E. Chang, and J. S. Douglas, Simulating quantum light propagation through atomic ensembles using matrix product states, Nat. Commun. 8, 1743 (2017).
  * (28) T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I. Cirac, and D. E. Chang, Quantum dynamics of propagating photons with strong interactions: a generalized input-output formalism, New J. Phys. 17, 113001 (2015).
  * (29) S. Prasad and R. J. Glauber, Coherent scattering by a spherical medium of resonant atoms, Phys. Rev. A 83, 063821 (2011).
  * (30) K. Lalumiere, B. C. Sanders, A. F. Van Loo, A. Fedorov, A. Wallraff, and A. Blais, Input-output theory for waveguide QED with an ensemble of inhomogeneous atoms, Phys. Rev. A 88, 043806 (2013).
  * (31) A. Albrecht, L. Henriet, A. Asenjo-Garcia, P. B. Dieterle, O. Painter, and D. E. Chang, Subradiant states of quantum bits coupled to a one-dimensional waveguide, New J. Phys. 21, 025003 (2019).
  * (32) L. Henriet, J. S. Douglas, D. E. Chang, and A. Albrecht Critical open-system dynamics in a one-dimensional optical lattice clock, Phys. Rev. A 99 (2), 023802 (2019).
  * (33) L. Ostermann, C. Meignant, C. Genes, and H. Ritsch, Super-and subradiance of clock atoms in multimode optical waveguides, New J. Phys. 21 (2), 025004 (2019).
  * (34) Y. Zhang, C. Yu, and K. Mølmer, Theory of subradiant states of a one-dimensional two-level atom chain, Phys. Rev. Lett. 122( 20), 203605 (2019).
  * (35) Y. Zhang and K. Mølmer, Subradiant Emission from Regular Atomic Arrays: Universal Scaling of Decay Rates from the Generalized Bloch Theorem, Phys. Rev. Lett. 125, 253601 (2020).
  * (36) J. Kumlin, K. Kleinbeck, N. Stiesdal, H. Busche, S. Hofferberth, and H. P. Büchler, Non-exponential decay of a collective excitation in an atomic ensemble coupled to a one-dimensional waveguide, Phys. Rev. A 102, 063703 (2020).
  * (37) D. F. Kornovan, N. V. Corzo, J. Laurat, A. S. Sheremet, Extremely subradiant states in a periodic one-dimensional atomic array, Phys. Rev. X 100, 063832 (2019).
  * (38) Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Inelastic scattering of photon pairs in qubit arrays with subradiant states, Phys. Rev. Lett. 123, 253601 (2020).
  * (39) J. Zhong and A. N. Poddubny, Classification of three-photon states in waveguide quantum electrodynamics, Phys. Rev. A 103, 023720 (2021).
  * (40) J. A. Needham, I. Lesanovsky, and B. Olmos, Subradiance-protected excitation transport, New J. Phys. 21, 073061 (2019).
  * (41) A. A. Makarov and V. I. Yudson, Quantum Engineering of Superdark Excited States in Arrays of Atoms, Phys. Rev. A 102, 053712 (2020).
  * (42) J. Zhong, N. A. Olekhno, Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Photon-mediated localization in two-level qubit arrays, Phys. Rev. Lett. 124, 093604 (2020).
  * (43) A. V. Poshakinskiy, J. Zhong, Y. Ke, N. A. Olekhno, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Quantum Hall phase emerging in an array of atoms interacting with photons, arXiv:2003.08257 (2020).
  * (44) A. V. Poshakinskiy, J. Zhong, and A. N. Poddubny, Quantum chaos driven by long-range waveguide-mediated interactions, Phys. Rev. Lett. 126, 203602 (2021).
  * (45) Y. Zhang, C. Yu, K. Mølmer, Subradiant Dimer Excitations of Emitter Chains Coupled to a 1D Waveguide, Phys. Rev. Research 2, 013173 (2020).
  * (46) A.N. Poddubny, Quasi-flat band enables subradiant two-photon bound states, Phys. Rev. A 101, 043845 (2020).
  * (47) P. Bienias, J. Douglas, A. Paris-Mandoki, P. Titum, I. Mirgorodskiy, C. Tresp, E. Zeuthen, M. J. Gullans, M. Manzoni, S. Hoerberth, D. Chang, and A. Gorshkov, Photon propagation through dissipative rydberg media at large input rates, Phys. Rev. Research 2, 033049 (2020).
  * (48) K. Stannigel, P. Rabl, and P. Zoller, Driven-dissipative preparation of entangled states in cascaded quantum-optical networks, New J. Phys. 14, 063014 (2012).
  * (49) H. Pichler, T. Ramos, A. J. Daley, and P. Zoller, Quantum Optics of chiral spin networks, Phys. Rev. A 91, 042116 (2015).
  * (50) P. Lodahl, S. Mahmoodian, S. Stobbe, P. Schneeweiss, J. Volz, A. Rauschenbeutel, H. Pichler, and P. Zoller, Chiral Quantum Optics, Nature 541 473 (2017).
  * (51) M. Wortis, Solitons and magnon bound states in ferromagnetic Heisenberg chains, Phys. Rev. 132, 85 (1963).
  * (52) T. Schneider, Solitons and magnon bound states in ferromagnetic Heisenberg chains, Phys. Rev. B 24, 5327 (1981).
  * (53) D. C. Mattis, The few-body problem on a lattice, Rev. Mod. Phys 58, 361 (1986).
  * (54) Fukuhara, T., Schaus, P., Endres, S. Hild, M. Cheneau, I. Bloch, and C. Gross, Solitons and magnon bound states in ferromagnetic Heisenberg chains, Nature 502, 76-79 (2013).
  * (55) F. Letscher and D. Petrosyan, Mobile bound states of Rydberg excitations in a lattice, Phys. Rev. A 97, 043415 (2018).
  * (56) C. D. Parmee and N. R. Cooper, Decay rates and energies of free magnons and bound states in dissipative XXZ chains, Phys. Rev. A 99, 063615 (2019).
  * (57) Y. Chougale, J. Talukdar, T. Ramos, and R. Nath, Dynamics of Rydberg excitations and quantum correlations in an atomic array coupled to a photonic crystal waveguide, Phys. Rev. A 97, 022816 (2020).
  * (58) F. Verstraete, V. Murg, and J.I. Cirac, Matrix Product States, Projected Entangled Pair States, and Variational Renormalization Group Methods for Quantum Spin Systems, Adv. Phys. 57, 143 (2008).
  * (59) U. Schollwöck, The Density-Matrix Renormalization Group in the Age of Matrix Product States, Ann. Phys. (Amsterdam) 326, 96 (2011).
  * (60) D. Reitz, C. Sayrin, R. Mitsch, P. Schneeweiss, and A. Rauschenbeutel, Coherence Properties of Nanofiber-Trapped Cesium Atoms, Phys. Rev. Lett. 110, 243603 (2013).
  * (61) R. Mitsch, C. Sayrin, B. Albrecht, P. Schneeweiss, and A. Rauschenbeutel, Quantum state-controlled directional spontaneous emission of photons into a nanophotonic waveguide, Nat. Commun. 5, 5713 (2014).
  * (62) P. Solano, P. Barberis-Blostein, F. K. Fatemi, L. A. Orozco, and S. L. Rolston, Super-radiance reveals infinite-range dipole interactions through a nanofiber, Nat. Commun. 8, 1857 (2017).
  * (63) M. Arcari, I. Söllner, A. Javadi, S. Lindskov Hansen, S. Mahmoodian, J. Liu, H. Thyrrestrup, E. H. Lee, J. D. Song, S. Stobbe, and P. Lodahl, Near-Unity Coupling Efficiency of a Quantum Emitter to a Photonic Crystal Waveguide, Phys. Rev. Lett. 113, 093603 (2014).
  * (64) I. Söllner, S. Mahmoodian, S. L. Hansen, L. Midolo, A. Javadi, G. Kiršanske, T. Pregnolato, H.El-Ella, E. H. Lee, J. D. Song, S.Stobbe, and P. Lodahl, Deterministic photon-emitter coupling in chiral photonic circuits, Nat. Nanotechnol. 10, 775 (2015).
  * (65) A. Goban, C.-L. Hung, S. P. Yu, J. D. Hood, J. A. Muniz, J. H. Lee, M. J. Martin, A. C. McClung, K. S. Choi, D. E. Chang, O. Painter, and H. J. Kimble, Atom-light interactions in photonic crystals, Nat. Commun. 5, 3808 (2014).
  * (66) J. D. Hood, A. Goban, A. Asenjo-Garcia, M. Lu, S. P. Yu, D. E. Chang, and H. J. Kimble, Atom-atom interactions around the band edge of a photonic crystal waveguide, PNAS 113, 10507 (2016).
  * (67) N. M. Sundaresan, R. Lundgren, G. Zhu, A. V. Gorshkov, and A. A. Houck, Interacting qubit-photon bound states with superconducting circuits, Phys. Rev. X 9, 011021 (2019).
  * (68) M. Mirhosseini, E. Kim, V. S. Ferreira, M. Kalaee, A. Sipahigil, A. J. Keller, and O.Painter, Superconducting metamaterials for waveguide quantum electrodynamics, Nat. Commun. 9, 1-7 (2018).
  * (69) M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang, and O. Painter, Cavity quantum electrodynamics with atom-like mirrors, Nature 569, 692-697 (2019).
  * (70) J. D. Brehm, A. N. Poddubny, A. Stehli, T. Wolz, H. Rotzinger, and A. V. Ustinov, Waveguide Bandgap Engineering with an Array of Superconducting Qubits, npj Quantum Mater., 6, 10 (2021).
  * (71) I. S. Besedin, M. A. Gorlach, N. N. Abramov, I. Tsitsilin, I. N. Moskalenko, A. A. Dobronosova, D. O. Moskalev, A. R. Matanin, N. S. Smirnov, I. A. Rodionov, A. N. Poddubny, and A. V. Ustinov, Topological photon pairs in a superconducting quantum metamaterial, Phys. Rev. B 103, 224520 (2021).
  * (72) M. Scigliuzzo, G. Calajo, F. Ciccarello, D.Perez Lozano, A. Bengtsson, P. Scarlino, A. Wallraff, D. Chang, P. Delsing, and S. Gasparinetti, Extensible quantum simulation architecture based on atom-photon bound states in an array of high-impedance resonators, arXiv:2107.06852 (2021).
  * (73) P. W. Anderson, Absence of Diffusion in Certain Random Lattices, Phys. Rev. 109, 1492 (1958).
  * (74) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entangle- ment, Rev. Mod. Phys. 91, 021001 (2019).
  * (75) N. Fayard, L. Henriet, A. Asenjo-Garcia, and D. Chang, Many-body localization in waveguide QED, arXiv:2101.01645 (2021).
  * (76) D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J. Kimble, Cavity QED with atomic mirrors, New J. Phys. 14, 063003 (2012).
  * (77) T.Shi, D. E. Chang, and J. I. Cirac, Multiphoton-scattering theory and generalized master equations, Phys. Rev. A 92, 053834 (2015).
  * (78) S. L. McCall and E. L. Hahn, Self-Induced Transparency by Pulsed Coherent Light, Phys. Rev. Lett. 18, 908 (1967).
  * (79) S. L. McCall and E. L. Hahn, Self-Induced Transparency, Phys. Rev. 183, 457 (1969).
  * (80) R. K. Bullough, P. J. Caudrey, J. C. Eilbeck, and J. D. Gibbon, A General Theory of Self-Induced Transparency, Opto-electronics 6, 121 (1974).
  * (81) P. O. Guimond, B. Vermersch, M. L. Juan, A. Sharafiev, G. Kirchmair, and P. Zoller, A Unidirectional On-Chip Photonic Interface for Superconducting Circuits, npj Quantum Inf. 6, 32 (2020).
  * (82) N. Gheeraert, S. Kono, and Y. Nakamura, Programmable directional emitter and receiver of itinerant microwave photons in a waveguide, Phys. Rev. A 102, 053720 (2020).
  * (83) M. Fleischhauer and M. D. Lukin, Quantum memory for photons: Dark-state polaritons, Phys. Rev. A 65, 022314 (2002).
  * (84) T. Peyronel, O. Firstenberg, Q.Y. Liang, S. Hofferberth, A. V. Gorshkov, T. Pohl, M. D. Lukin, and V. Vuletic, Quantum nonlinear optics with single photons enabled by strongly interacting atoms, Nature 48, 57-60 (2012).
  * (85) Z. Bai, C. S. Adams, G. Huang, and W. Li, Self-Induced Transparency in Warm and Strongly Interacting Rydberg Gases, Phys. Rev. Lett. 125, 263605 (2020).
  * (86) E. Shahmoon, P. Grisins, H. P. Stimming, I. Mazets, and G. Kurizki Highly nonlocal optical nonlinearities in atoms trapped near a waveguide, Optica 3, 031017 725 (2016).
  * (87) T. Caneva, J. S. Douglas, and D. E. Chang, Photon Molecules in Atomic Gases Trapped Near Photonic Crystal Waveguides, Phys. Rev.X 6, 031017 (2017).
  * (88) N. Stiesdal, J. Kumlin, K. Kleinbeck, P. Lunt, C. Braun, A. Paris-Mandoki, C. Tresp, H. P. Büchler, and S. Hofferberth, Observation of Three-Body Correlations for Photons Coupled to a Rydberg Superatom, Phys. Rev. Lett. 121, 103601 (2018).
  * (89) N. Stiesdal, H. Busche, J. Kumlin, K. Kleinbeck, H. P. Büchler, and S. Hofferberth, Observation of collective decay dynamics of a single Rydberg superatom, Phys. Rev. Research 2, 043339 (2020).
  * (90) B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, Matrix product operator representations, New J. Phys. 12, 025012 (2010).
  * (91) G. M. Crosswhite, A. C. Doherty, and G. Vidal, Applying matrix product operators to model systems with long-range interactions, Phys. Rev. B 78, 035116 (2008).
  * (92) F. Fröwis, V. Nebendahl, and W. Dür, Tensor operators: Constructions and applications for long-range interaction systems, Phys. Rev. B 81, 062337 (2010).