Dynamics of many-body photon bound states in chiral waveguide QED

We theoretically study the few- and many-body dynamics of photons in chiral waveguides. In particular, we examine pulse propagation through a system of $N$ two-level systems chirally coupled to a waveguide. We show that the system supports correlated multi-photon bound states, which have a well-defined photon number $n$ and propagate through the system with a group delay scaling as $1/n^2$. This has the interesting consequence that an incident coherent state pulse breaks up into different bound state components during propagation, which can become spatially separated at the output by photon number in a sufficiently long system. For sufficiently many photons and sufficiently short systems, we show that linear combinations of $n$-body bound states recover the well-known phenomenon of mean-field solitons in self-induced transparency. Our work thus covers the entire spectrum from few-photon quantum propagation to genuine quantum many-body (atom and photon) phenomena, and ultimately the quantum-to-classical transition. Finally, we demonstrate that the bound states can undergo elastic scattering with additional photons. Together, our results demonstrate that photon bound states are truly distinct physical objects emerging from the most elementary light-matter interaction between photons and two-level emitters. This opens the door to studying quantum many-body physics and soliton physics with photons in chiral waveguide QED.


I. INTRODUCTION
Generating quantum many-body states of light remains one of the outstanding challenges of modern quantum optics [1]. Such many-body states of light are of fundamental physical interest as they arise from nonequilibrium systems with strong interactions between light and matter. On the other hand, they also promise to form novel resources for quantum technologies, for example, in quantum enhanced metrology [2]. The main obstacle in the pursuit of generating such many-body states has been the difficulty in developing one-dimensional systems with a sufficiently strong nonlinear response at the few-photon scale [3]. Recently however, significant progress has been made in creating an ideal light-matter interface between atoms or artificial emitters coupled to a one-dimensional continuum of photons at optical [4][5][6] and microwave frequencies [7][8][9]. Such an interface creates a highly nonlinear medium as photons propagating in a waveguide interact deterministically with atoms. Systems of this kind have been thus far used to propose or demonstrate the generation of states of photons with strong two-or three-body correlations [10][11][12][13][14][15][16][17][18][19][20].
Studies of photon correlations in these systems typically consider steady-state driving and photon correlations are subsequently measured in relative coordinates. On the other hand, here we show that pulse propagation through two-level systems (TLSs) strongly coupled to a waveguide leads to very distinct temporal features which reveal the underlying dynamics and has the potential to generate temporally ordered many-body states of light.
In particular we theoretically consider the propagation of pulses of coherent and Fock states of light through a waveguide to which N TLSs are chirally coupled. We show that photons undergo a Wigner delay [21,22] -a delay of the pulse center due the optical excitation transferring to the atom and back to the waveguide -which is dependent on the number of photons. Therefore, incident pulses break up into a state that is temporally ordered by its photon correlations. For example, as shown in Fig. 1(a), an incident coherent field can produce a pulse with three-photon correlations followed by two-photon correlations states followed by uncorrelated photons. The underlying physics that causes this time delay is examined by considering the many-body photonic scattering eigenstates of a TLS chirally coupled to a waveguide.
Of central importance are the class of bound eigenstates, where two or more photons propagate together. We show that the photon bound states propagate past each atom with a photon-number-dependent delay of τ n = 4/(n 2 Γ),where Γ is the decay rate into the waveguide, as can be understood in terms of absorption and stimulated emission of a photon as shown in Fig. 1(b). Since the interaction is chiral, the delay is also proportional to the number of emitters N in the waveguide. We also show that the pulse distortion encountered by an n-photon bound state scales as n −6 . Pulses of higher order bound states therefore propagate with negligible distortion, i.e. much like a soliton. Indeed we show that the classical mean-field limit of a soliton solution emerges from a combination of photon bound states with n 1. Going beyond mean field, we show that with a mesoscopic number of photons n but a large number of atoms  1. (a) N two-level atoms (blue circles) are chirally coupled to a waveguide with decay rate Γ and driven by an input Gaussian pulse, which can be a coherent or Fock state. The light pulse propagates with a correlation-number-dependent group velocity leading to an output state where one-, twoand three-photon bound states are spatially separated. (b) Schematic of the bound states propagation. When an nphoton bound state is scattered by an atom, it re-emits the absorbed photon with a stimulated emission rate Γn, which also coincides with the inverse of its width. Since only a single photon out of n is delayed by an amount 4/(Γn), the pulse preserves its shape but is delayed by τn = 4/(Γn 2 ). N 1, the soliton breaks apart into a temporally ordered many-body quantum state of light. In this limit mean-field theory fails and our full quantum many-body treatment provides a simple description of photon propagation in waveguide quantum electrodynamics (wQED).
In mean-field theories solitons are known to be highly stable object which are unaffected by external perturbations. Here we show that similar properties exist for fewphoton bound states. We outline how one can conduct scattering experiments between photon bound states and individual photons. In the considered scattering experiments the bound state is deflected by the interaction but is otherwise unperturbed by it. Together the results obtained here demonstrates that the bound states should be considered as truly distinct physical entities emerging from the underlying light-matter interactions.
This manuscript is arranged as follows: in Section II we introduce the model for chiral wQED and outline various theoretical approaches for computing the dynamics through II A mean-field theory, II B the photon scattering eigenstates, and II C the matrix-product states (MPS). In Section III we compute how the input pulse propagates through the medium and compute the representation in terms of photon bound states. This is followed by Section IV which shows that the many-body photon bound states can be used to construct the mean-field soliton solutions obtained in self-induced transparency. Finally, in Section V we show that photon bound states can undergo elastic scattering with individual photons modifying the delay of the bound state but otherwise leaving it unaltered, much like classical solitons.

II. MODEL
We consider a system of N two-level systems (TLS) chirally coupled coupled to a linearly dispersive onedimensional bath of photons. The Hamiltonian for this system ( = 1) iŝ where all integrals are over , σ ± i are Pauli operators for the ith TLS,â(x) (â † (x)) is a photon annihilation (creation) operator at position x, the group velocity is set to unity, and the energy has been renormalized to that of the TLS. The positions of the TLSs x i are arbitrary and only change the overall phase of the outgoing field. In this Hamiltonian the first term captures the free-propagation of photons in the waveguide with linear dispersion, while the next terms take into account the interaction between the photons and the emitters. In the absence of the interaction terms the photonic eigenstates of the Hamiltonian are plane waves with wavevector k and frequency ω = k. Equation (1) constitutes the typical scenario for chiral wQED [23]. In this manuscript, we are interested in describing the propagation of a multiphoton input field through the system. This goal can be achieved by exactly computing the scattering eigenstates of Hamiltonian (1), as illustrated in Subsection II B.
In addition to computing the eigenstates of (1), we introduce an equivalent formulation that is well-suited for using the MPS technique that we are going to introduce later in Subsection II C. This approach is based on the observation that the transmitted field depends directly on the emitters evolution. This can be seen by formally integrating the Heisenberg equation for the field operatorâ(x, t), that provides, within Born-Markov approximation, the following generalized input-output relation for the transmitted field [24,25] where we have definedâ out (t) =â(x + N , t) as the output field measured right after the last atom. Note that, within this approach, we assume the input field E in (t) to be a classical coherent field on resonance with the atomic transitions. With these assumptions the emitter dynamics driven by the input field is known to be described by a purely dissipative chiral master equation (ME) of the form [26,27] (see supplementary material (SM) for more information): Here, is the effective Hamiltonian, which provides the non-Hermitian collective evolution of the emitters, while the term gives the coupling of the emitters to the input field.
The combination of Eqs. (2) and (3) provide a full description of the photon propagation through the chiral medium. In particular the spin dynamics can be efficiently solved by making use of an MPS ansatz [28,29] as recently described in [25]. As will be shown in the following, this approach will allow us to fully explore the limit of many photons and large atomic arrays, a scenario that is challenging to simulate with standard numerical techniques.
A. Mean-field Theory and Self-induced Transparency Before considering the full many-body dynamics of the Hamiltonian (1), we consider the system within the mean-field limit. We present this mean-field limit to contrast its predictions with the full many-body theory presented below.
The first treatment of (1) within mean-field theory dates back to the work on self-induced transparency (SIT) [30][31][32]. In these early experiments, gasses of twolevel atoms were excited by short intense laser pulses. Although the atoms are not ideally coupled to a single waveguide mode in such systems, the laser pulses are sufficiently short that decay channels to modes other than the laser mode can be neglected. Furthermore, both the weak coupling of the atoms to this mode and the high intensity of the laser means that one can consider the atoms as a spin continuum under the mean-field approximation where quantum correlations between the atoms and the light field can be neglected. Under these approximations, the equations of motion give the SIT equations Here, a( the expectation values of the spin operators in the continuum limit. These nonlinear equations of motion have SIT soliton solutions. Following the treatment in [32], for a resonant pulse the field can be taken to be real, and one can map the equations of motion onto a nonlinear pendulum equation that can be solved exactly. This gives the fundamental soliton solution for the field wheren is the number of photons in the field,(or more precisely, the total energy in the original SIT work), and V =n 2 Γ/(4ν) is the relative velocity in the SIT medium which depends on the field strengthn, the decay rate Γ, and the linear density of the gas ν. An important feature of the solitonic solution (6) is that the integrated Rabi frequency Ω = 2 √ Γ dta(x, t), which is proportional to the area under the pulse, is fixed by the relationship between the pulse amplitude and pulse width, and always evaluates to be 2π. This corresponds to a full Rabi cycle of complete excitation and subsequent deexcitation. The SIT soliton therefore can be physically interpreted as a rapid excitation and de-excitation of the atoms, which suppresses spontaneous emission of the excited state and thus makes the medium transparent.
Inspired by the apparent dependence of velocity on photon number, it is interesting to ask whether this property extends to the few-photon limit, thus enabling, e.g., photon-number separation at the output. A full quantum treatment is necessary to answer this question, which we turn to in the following sections.

B. Many-body Scattering Eigenstates
In contrast to the mean-field treatment we now consider the full many-body eigenstates. Since the Hamiltonian (1) preserves the combined number of atomic and photonic excitations, eigenstates with different numbers of excitations decouple. By computing the eigenstates in the one-and two-excitation subspaces, one can generalize the result to an arbitrary number of excitations. This technique is often referred to as Bethe's ansatz [33] and is used to diagonalize a class of one-dimensional manybody Hamiltonians [34]. In particular, it has previously been used to diagonalize the Hamiltonian in Eq. (1) [35]. Since we are interested in the state that emerges after interaction with the TLSs, we are interested in the scattering eigenstates. These are the photon states that interact with the TLSs and emerge unchanged apart from an overall transmission coefficient. The n-body scattering eigenstates have the form where C k,n,S is a normalization constant which varies with wavevector k, excitation number n, and the type of eigenstates S, where S labels different states as explained below;â † (x) =â † (x 1 )â † (x 2 ) . . .â † (x n ); d n x = dx 1 dx 2 . . . dx n ; and ↔ indicates summing over all n! permutations of x i to symmetrize the wavefunction. The energy of the eigenstates is E = n i k i . Upon scattering off all N emitters, the eigenstates are multiplied by the eigenvalue t N k = n j=1 t N kj , where t k = (k − iΓ/2)/(k + iΓ/2). Since, by assumptions, the system does not contain any dissipation and is chiral, all the transmission coefficients have |t k | = 1. This means that the transmission coefficients simply multiply the eigenstate by a phase. For example a single photon on resonance undergoes a π-phase shift with t 0 = −1. Importantly, the phase that is imparted on the eigenstate varies with k, i.e., the TLSs introduce dispersion to the system. As we soon show, different types of eigenstates also accumulate different phases. We note that the output states is here described using the position coordinates x. This is equivalent to using a time variable t in (2) as we have set the group velocity to unity.
In addition to different eigenstates for different values of n, there are also different possible types of eigenstates within different excitation-number manifolds. In the single excitation subspace there is only one type of eigenstate and they are characterized by real wavenumbers k and are fully extended in space. For n = 2, the wavenumbers k 1 and k 2 can either both be real which gives rise to a fully extended solution of the form (7). They can also assume complex values k 1 = E/2 + iΓ/2 and k 2 = E/2 − iΓ/2 which are called strings [34]. These values give additional valid solutions. Since the wavevectors can form complex conjugate pairs, the eigenstates become localized in a relative coordinate while being fully extended within the center-of-mass coordinate. This localization is associated with the formation of photon bound states that manifest in the bunching of two photons which travel together during their propagation [36,37].
For n > 2 the m-body string (m ≤ n) has the wavevector k j = K/m − i [m + 1 − 2j] Γ/2 with j = 1, 2, . . . , m where K is the total energy of the m-string. In total, the n-photon manifold has p(n) string combinations where p(n) is the number of partitions of n. For example, for three photons, there is a completely extended scattering eigenstate, a completely bound eigenstate, and a hybrid state with two bound photons and one extended photon. The different string combinations give the different types of scattering eigenstates S by substituting the complex wavevectors into (7). The transmission coefficient is also obtained by substituting the complex wavevectors into the expression for t k .
For n photons, the n-string state is a fully localized bound state with energy E, where C n,B = Γ n−1 (n − 1)!/(2πn). The transmission coefficient for the n-photon bound state is then, Importantly, the phase of the transmission coefficient varies with n, i.e., the system has a photon-numberdependent dispersion. With all the eigenstates at hand, the scattering matrix for interacting with all N emitters in the n-photon manifold can be formally written aŝ where the sum over S is a sum over the different string combinations of the n-photon manifold. We note that the eigenstates are orthogonal, thus the scattering matrix for N emitters simply requires taking the eigenvalue to the N -th power. Since the number of string combinations increases as p(n), the number of terms in the sum increases exponentially for large n [35,38]. In this manuscript we compute the full output states for up to n = 3 using this formalism. We note that a formalism exists where one does not have to sum over string combinations [39]. Nevertheless the form used here is particularly insightful as it gives direct access to the number-dependent transmission coefficient which, as we show, plays a central role in understanding the many-body pulse propagation.

C. MPS Ansatz
In order to study the dynamics for stronger input pulses (n > 3) than the one computed with the Smatrix formalism, we solve equations (2) and (3) using an MPS ansatz. Specifically, the system evolution can be solved either by directly solving the ME (3) [40] (method used for Fig. 2(d) and Fig. 3) or by using a quantum trajectories algorithm where the state of the system evolves under the effective Hamiltonian (4) and stochastically experiences quantum jumps [25] (method used for Fig. 2(f)). In both cases an MPS representation is applied either to the quantum state or to a vectorized form of the density matrix. Here for simplicity we limit the discussion to the former while the latter is discussed in the supplementary material. The MPS ansatz consists of reshaping the generic quantum state where, for each specific set of physical indices {i 1 ,ī 2 , ..ī N }, the product of the A ij matrices gives back the state coefficient ψ i1,ī2,..ī N . Each matrix A ij has dimension D j−1 × D j known as the bond dimension and finite edge boundary conditions are assumed by imposing  For instance, if D j = 1 for all j, the matrices A ij are scalars and the state (11) reduces to a product state with no entanglement. For arbitrary states the bond dimension grows exponentially with the number of particles. The advantage of the MPS ansatz is that, in many physical scenarios as the one considered here, the entanglement can grow slowly with the system size allowing an efficient description of the state in terms of a smaller bond dimension [29]. An important figure of merit of the efficiency of the MPS ansatz is given by D max , the maximum bond dimension, used to represent the system during the entire time evolution (see SM for more information). We will make use of this quantity in Sec. IV to quantify the amount of many-body correlations present in the system.

III. MANY-BODY PULSE PROPAGATION
We are interested in studying multi-photon propagation through the chirally coupled array. Here, we consider coherent and Fock input states with mode creation operatorâ † in = dxE(x)â † (x), and we evaluate the transmitted field with the two methods described in the previous section. In particular, for the exact solution we compute the transmitted photon state using the eigenstates for up to three photons while for higher excitations we make use of the MPS ansatz. In Fig. 2 we consider the propagation of a Gaussian photonic mode with am- Here the x coordinate is chosen to be in the reference frame of the pulse propagating in the absence of emitters. In Fig. 2(b) we see that the two-photon bound state comes out earlier than the extended state. The two-photon bound state thus clearly propagates with a faster velocity than the extended state. The bound state also undergoes significantly less broadening and distortion. For the three-photon transport in Fig. 2(c) again the extended state is distorted and delayed, while the three-photon bound state has less significantly distortion and delay. In the three-photon manifold there is also a string combination that forms a hybrid state where one photon is completely extended while the other two are bound. The evolution of this state is as follows determined by the individual components of the state: the two bound photons propagate in a similar manner to the two-photon bound state in Fig. 2(b), while the extended photon propagates like a single photon. This separation of the propagation of the bound and extended parts of this state can be shown explicitly in the long pulse limit σΓ 1 (see SM). This is significant because it implies that in order to understand the pulse evolution one does not have to understand the behaviour of all the different string combinations S. Rather, one can focus on simply understanding the behaviour of the bound states. On the other hand, for short pulses this separation is not completely true because one needs to include the effect of interactions between the components, see Sec. V.

A. Evolution of bound states
To better quantify the difference in propagation observed above, let us compute the pulse evolution in the center-of-mass coordinate. This can be done by using the form of the n-photon bound state and its transmission coefficient given in Eqs. (8) and (9). Within the n-photon manifold the projection of an input Gaussian state and the bound state is where g(x J ) = i<j |x i − x j | written in Jacobi coordinates. The center-of-mass evolution is clearly determined by the second integral. This now has the standard form of Gaussian pulse propagation through a linear dispersive medium. Defining t N E,n ≡ e iN φ(E) , the first to third derivatives of φ(E) give the delay, broadening, and distortion that the pulse undergoes per emitter respectively. The delay per emitter is τ n (k 0 ) = Γ/(k 2 0 + n 2 Γ 2 /4). This is largest for a resonant pulse τ n = 4/(n 2 Γ). (13) This gives the Wigner delay imparted on an n-photon bound state by a single emitter in wQED: the photons propagate with a number-dependent velocity. By taking higher-order derivatives we also compute the pulse broad- , which is zero on resonance. The third-order pulse distortion term on resonance is d = −32/(n 6 Γ 3 ). The pulse distortion is thus drastically reduced for higher-order bound states. This indicates that many-photon bound states suffer negligible pulse distortion while propagating through the array of nonlinear and dispersive atoms.
In order to verify that indeed the physics of the bound states dominates the wave-packet evolution we also compute the evolution of a coherent input pulse as shown in Figure 2(d). Here, the pulse width Γσ = 3 √ 2 is the same as in Figs. 2(a)-(c), while the average photon number in the pulse isn = 0.5. We compute the output both by truncating the coherent state to three photons, and solving exactly using Eq. (10), or by solving Eqs. (2) and (3) with the MPS algorithm. The evolution of the bound states is seen in the position of the peaks in the power distribution â † (x)â(x) as well as in the difference between the power and the n-th order correlation functions n . The localized nature of the bound states in the relative coordinates is shown in Fig. 2(e) where we computed the two-point second order correlation function . While exact analytical calculations become unfeasible forn 1, the validity of our arguments and the importance of the bound states can still be seen in MPS simulations. For example, in Fig. 2(f), we calculate the transmitted power for an input coherent staten = 8 and N = 60. Here, for this particular system length, the photon-numberdependent Wigner delay clearly manifests itself as separate peaks for up to six-photon bound states.
This low-distortion propagation of the bound states can be intuitively explained by returning to the simple schematic shown in Fig. 1(b). When a multi-photon bound state propagates trough the atomic array one of the photons in the packet can be absorbed and re-emitted by the atom. This process occurs on a time scale ruled by the inverse of the photon-number dependent stimulated emission rate that coincides with the bound-state packet width, ∆t n ∼ 1/(Γn), allowing the pulse to preserve its shape. This continuous absorption and re-emission of photons during the bound-state propagation leads to a time delay of one out of n photons by an amount 4/Γn, leading to the group delay in Eq. (13).
Finally, while these results obtained in the quantum regime were initially inspired by the classical regime of SIT, we now establish a concrete relationship between the SIT solitons and quantum bound states.

IV. CONNECTION TO THE SIT SOLITON
In the previous sections we have shown how the Hamiltonian (1) leads to the SIT solitonic solutions in the mean-field limit, and that the full quantum mechanical treatment of this Hamiltonian predicts correlationordered photon propagation. In this section we aim to bridge the gap between these two regimes: first we show that indeed the many-body theory reduces to the meanfield result in the limit of large photon number. Secondly, we derive the quantum corrections to the mean-field results which become relevant when both the number of photons n and the number of emitters N are large. Finally, we push the numerical simulations to the manyphoton limit to verify the theoretical predictions.
Let us consider a wave packet composed of a linear combination of many-body bound states. This is expressed with the ansatz We later show that this is the expected form of the state for a high-power coherent input with a large spectral width. Such a state can be probed by measuring either the field ψ|â(x)|ψ or the m-th normally ordered observable, which we consider in the center-of-mass coordinate ψ|(â † (x)) m (â(x)) m |ψ . We compute these observables in the limit where the average photon number is largē n 1, the pulses are spectrally broad σ 1/Γ and the order of the correlation function is much less than the photon number m n. These give (see SM for full calculations), wheren is the average number of photons in the pulse. These observables reproduce the fundamental soliton solution of SIT [31,32,41] given in Eq. (6). Self-induced transparency is thus the classical limit of the photon bound state when the photon number becomes large, or conversely photon bound state are simply the quantum limit of a soliton, a quantum soliton. Just like the SIT solitons, the integrated Rabi frequency of the pulse is Ω = 2 √ Γ dx â(x) = 2π, i.e., the intense pulse of light rapidly excites and de-excites the emitters resulting in a 2π Rabi oscillation. Since we consider high photon numbers n 1, the effects of pulse distortion due to higherorder dispersion are negligible and interaction with the emitters simply maps x → x + 4N/(n 2 Γ).

A. Beyond mean-field theory
As a lowest-order approximation, the variation in photon number n making up the pulse (14) can be ignored, and the state will simply propagate at a reduced speed dictated by the mean photon number. This leads to the SIT soliton prediction, that simply maps x → x + 4N/(n 2 Γ). However, for a coherent state, the uncertainty in the photon number scales as ∆n ∼ √n , which leads to a gradual broadening of the pulse due to the different photon number components accumulating different time delays. This can become significant for a sufficiently large number of emitters N . This broadening is not captured by the mean-field theory. The breakdown of the mean-field theory therefore occurs when the difference in the delay of then and then + ∆n photon bound states becomes on the order of the of the width of thē n photon bound state. This gives (τn+∆n−τn) ∆t(n) ∼ N/n 3/2 . This means that when N/n 3/2 1 the mean-field theory breaks down even for large input power. This provides the boundary between mean-field theory and genuine quantum many-body dynamics.
Whenn 1 and N n 3/2 , the Fock-state-dependent delay can be included in the expression for the field and the power to give with equivalent expressions for higher-order correlation functions. We note that similar expressions exist for the bound state of the nonlinear Schrödinger equation with an attractive interaction [42,43]. Equations (17) provide a simple description of the observables of a quantum many-body state of light. In order to investigate the full transition from the multi-photon bound-state propagation to the formation of the SIT solitons we make use again of the master equation simulation. In Fig. 3    On the other hand it is important to emphasize the difference in physical effects. While the SIT solitons can be fully described by a mean-field semi-classical treatment, the formation of distinct bound-state peaks is characterized by a highly-correlated state of light and represents the break down of the mean-field solution due to quantum effects.
Within the MPS ansatz, one natural way to characterize the amount of correlations in the system is to allow the maximum truncated bond dimension D max to vary, and to record the value D th max (N a ,n) at which the truncation error exceeds some acceptable threshold value (see SM). In Fig. 3, we show this quantity as a function of the solitonic input pulse strengthn. We see that the breakdown of the mean-field description occurs in the regime where bound state formation occurs and amount of correlation in the system is high (large values of D max ), while in the limit of largen the bond dimension tends to shrink approaching the mean-field limit (D th max = 1).

V. SOLITON INTERACTIONS
So far we have characterized the propagation of light through ideal media and show that it can be understood in terms of the photon bound state. To fully characterize and understand these objects it is important to also investigate their interactions and robustness to disturbances. To this end, we now make a preliminary investigation of a scattering experiment between a single photon and a two-photon bound state. We consider the input state |in = C d 3 x a † (x)|0 φ(x 1 , x 2 , x 3 ) with C being a constant and i.e., a product state composed of a two-photon bound state and a single photon which are centered at a 1 and a 2 respectively with a 1 < a 2 . Figure 4 shows the evolution of the pulse delay for this state as it propagates through the ensemble, i.e. for different N . Here since a 1 < a 2 and the Wigner delay is larger for single photons τ 1 > τ 2 , the bound-state photons catch up to and overtake the single photon. In this process the photons interact when the two parts overlap. The interaction causes a change in the Wigner delays which is seen as kinks in the lines in Figs. 4(a)-(b). The region of interaction is highlighted by the third-order correlations shown in Fig.4(c). After the interaction N 10 the lines in Fig. 4 (a) continue with the same slopes as before the interaction, signifying that there is still a two-photon bound state and an unbound photon. The collision between the bound state and the free photon is thus elastic and the bound state is stable against external influence. Potentially, the interaction identified here could be used to realize passive photonphoton gates for quantum computing.

VI. CONCLUSION
In conclusion we have shown that chiral wQED platforms provide a highly nonlinear medium suited for exploring nonlinear optics at the quantum level. From the most elementary light-matter interaction, the interaction between photons and two-level systems, emerges correlated photonic states. These bound states are truly distinct physical objects with their own dispersion relation and are stable against external influence. Our work provides a clear recipe for how these features can potentially be observed in experiments. In the limit where the number of photons is high the bound state approach approach the well-known soliton solution of SIT. Our full quantum description on the other hand covers the entire spectrum from few-photon quantum propagation to genuine quantum many-body (atom and photon) phenomena, and ultimately the quantum-to-classical transition. In particular the analysis highlights how the mean-field solution with weak quantum correlations break down through a region of maximal quantum correlations, until finally resulting in a state with weaker correlations. This phenomena can be used to generate quantum many-body states of light with strong position-dependent photon-photon correlations. With strong light-matter interaction available in current circuit QED platforms, the physics presented here can be observed with current setups provided they are modified to exhibit a chiral interaction. Supplementary Material: Dynamics of many-body photon bound states in chiral waveguide QED

S1. ONE-, TWO-AND THREE-PHOTON TRANSPORT USING SCATTERING EIGENSTATES
In this section we detail how we compute the ouput photon states when N emitters are impinged upon by one, two, or three photons. We also combine these scattering calculations to consider scattering of a coherent state with up to three photons.

A. Single-Photon Transport
In the few-photon transport calculations we we consider Gaussian input states. The normalized single-photon input state has the form with its Fourier representation The creation operatorsâ † (k) are eigenstates of the single-photon scattering operator. After N scattering events, we haveâ † (k) → t N kâ † (k), with the transmission coefficient Throughout this document we consider an ideal 1D reservoir without losses, and we also set v g = 1. The scattered state is thus This is computed numerically. Most of the integrals that appear in the few-photon transport calculations are Fourier integrals, which can be efficiently computed using Fast Fourier Transforms.

B. Two-Photon Transport
The two-photon scattering matrix has two types of eigenstates, extended states |W E,∆ and bound states |B E [36]. The scattering matrix has the form where all integrals range from −∞ to ∞. The two-photon scattering eigenstates can be obtained from Bethe's ansatz in the main text (7). In a position-space representation they are Here, and and x c = x1+x2 2 , x = x 1 − x 2 , E = k + p is a two-photon detuning, and ∆ = k−p 2 is a difference in photon energies, where k and p are photon detunings.
The input state for two photons is Projecting this on the two-photon states gives where Erfc = 1 − Erf, where Erf is the error function, and D F is Dawson's Function. The output state is then written in terms of two integrals given by the bound-state and extended-state contributions and (S14) The integrals in these two functions are computed numerically.

C. Three-photon transport
Using the same rules rules we construct the three-photon scattering matrix for N emitterŝ where The eigenstates are with the normalized real-space representations, (S18) The three-photon Gaussian input state is Both the projection of the input state on the scattering eigenstates and the subsequent integrals in (S15) are performed numerically.

D. Output states and observables
We now have one-, two-, and three-photon output states which we express as (S20) These are also used to compute the output for coherent states We are interested in computing the normally ordered expectation values and correlations â † (x)â(x) , . For a coherent state with up to three photons, these give the expressions (S22)

S2. SEPARABILITY OF THREE-PHOTON HYBRID STATES
The behaviour of the Hybrid state in Fig. 2(c) of the main text is significant. Its shape indicates that the bound two-photon part of this state propagates with the same effective velocity as the actual two-photon bound state, while the extended photon in the hybrid state evolves like a single photon.
In this section we show that, for spectrally narrow input states, the evolution of the Hybrid state follows that of an independent two-photon bound state and an independent single photon. We consider the input state |in 3 given in (S19) that is on resonance k 0 = 0. The hybrid state has the form given in (S18). The contribution of the hybrid state to the full three-photon output state is We now focus on the inner product H K,k |in . Noting that the input state is separable and therefore already fully symmetric we can relabel the 3! terms in the Hybrid state to make them identical. We end up with Using the two-body coordinates x c = x 1 + x 2 /2 and x = x 1 − x 2 allows writing the integral as Unfortunately all the integrals cannot be computed exactly. We now perform the integrals over the x coordinate by considering Γσ 1. The exponentially decaying functions decay with Γ, while the Gaussian functions have a width ∼ 1/σ. When taking Γσ 1, we can make the approximation e Γ/2x e −x 2 /(2σ 2 ) ∼ e Γ/2x for x < 0. We can then write the integral over x as where a = x 3 − x c . We then compute the integrals The integral in (S25) then becomes where we have used I 1 + I 2 = −4/Γe −Γ|a| and I 1 − I 2 = 4/Γ sgn (a)(1 − e −Γ|a| ). Unfortunately we are unable to compute I(K, k) exactly. Instead we compute the projection of this function on the K-axis, i.e.Ĩ(K) = dkI(K, k). We show that this function is much narrower than the linewidth. Consequently, conservation of energy E = K + k will imply that its projection on the k-axisĪ(k) = dKI(K, k) is also narrower than the linewidth. After computing the integrals and some manipulation, in the limit Γσ 1, we obtain where Erfi is the imaginary error function (which is actually real-valued). In the limit K → 0 and Γσ 1 the function is On the other hand, in the limit Kσ 1 and Γσ 1, using the large argument approximation e −x 2 /4 Erfi (x/2) = 2/( √ πx) + O(1/x 3 ) we have The Gaussian part becomes negligible here and the function therefore scales as 1/K 2 σ 2 which is much less than one when K → Γ. Within the linewidth of the emitter, the function is largest when Kσ 1, where it is well-approximated by a Gaussian function. Since the tails of this function scales as 1/K 2 σ 2 , they only serve to slightly broaden the Gaussian function. This implies that (S30) is a reasonable approximation for the function.
This means that the Gaussian part of the function is dominant over the other parts. This part of the function I(K, k) comes from the integrals I 0 and I 3 and can actually be computed analytically when Γσ 1. It is given by This then gives This projection has two very useful features. The Gaussian part is separable in K and k and this part is highly localized about the origin. In fact the function multiplying the Gaussians can be approximated by a constant in K and k when Γσ 1. Substituting the projection H K,k |in into the (S23) gives the hybrid wavefunction The integrals over K and k are separable and can be computed individually. They are simply inverse Fourier transforms of the Gaussian functions multiplied by their transmission coefficients. The two-photon bound part of the hybrid state therefore evolves according to the two-photon bound state transmission coefficient, while the single photon part evolves according to the single-photon transmission coefficient. One can therefore interpret the hybrid eigenstates as a twophoton bound state and a single photon which are weakly interacting. The probability of interaction then becomes negligible in the long pulse limit.

S3. ELECTRIC FIELD OF BOUND STATES AND RELATION TO SELF-INDUCED TRANSPARENCY
In this section we compute the electric field of a linear combination of bound states. In particular, we show that in the classical limit, where the sum is composed of a sum of many-body bound states, the field resembles that of a self-induced transparency soliton.
We begin by considering an arbitrary state composed of a sum of bound states The electric field is then proportional to The bulk of this calculation involves computing the matrix element in the integrand. Using the expression for the n-photon bound state in (8) of the main text, the matrix element is where K n = E − nE n−1 . The bound-state wave functions are highly symmetric, which can be used to greatly restrict the domain of integration. This symmetry will be exploited to reduce the exponentials with negative absolute value arguments to exponentials without absolute values, but with arguments that are always negative. Noting that the integral is invariant under permutation of any x i ↔ x j , we can restrict the domain of integration to x 1 > x 2 > . . . x n−1 and correspondingly multiply by (n − 1)!. This gives There is still an exponential factor with an absolute value in its exponent. This is where θ(x) is Heaviside's step function. In general, this contains 2 n−1 terms when fully expanded, but most of these terms integrate to zero because the integral is over x 1 > x 2 > . . . > x n−1 . In fact, there are only n terms that contribute, where the notation θ( . This then gives n many-body integrals each over n − 1 variables, The j-th integral is (S43) The n − 1 integrals thus split into two sets of integrals: the first j − 1 integrals and the second n − j. The first j − 1 integrals have a form giving where α k = Γ [n − 2 − 2(j − 1) + 1/2]−iK n /n and the sums over α k can be easily evaluated. The latter n−j integrals must be evaluated in reverse order. They givẽ where we have defined K n = K n /n. We note that in I j (x) the exponents containing Γ all cancel each other.
Putting everything together the matrix element can be written Using the properties of the complex Gamma function Γ (z) one can express The entire expression becomes which leads to If the Gaussian function is slowly varying with respect to the sech function, one can expand it about K n = 0 giving 1 + O(K n ). Taking this limit and evaluating the Fourier integral gives Defining the integration Rabi frequency as Ω = 2 √ Γ dx â(x) , one obtains The normalization of the wavefunction implies 1 = n dE|c n (E)| 2 = n |c n | 2 √ nπ σ and therefore the integrated Rabi frequency in the large n limit approaches Ω = 2π. (S58) From (S56), in the limit of a large average photon number and a narrow photon distribution, the electric field can also be approximated as wheren is the average photon number.

S4. NORMALLY-ORDERED CORRELATION FUNCTION OF BOUND STATES AND MEAN-FIELD THEORY
Here we compute the m-th-order normally ordered correlation function of the bound states (n ≥ m). Following the previous section one obtains an integral As before, one must reduce the integral to the domain x 1 > x 2 > . . . > x n−m and then expand the product and reduce it to its n − m + 1 non-zero terms. Computing all the integrals using the techniques of the previous section, one obtains where now K = E − E .
A. Reduction to mean-field theory (n m) We now show that if the particle number n is much large than the order of the correlation function m, the outcome of the measurement probes the mean fields and the correlation function factorizes. One can start by making the approximation Since the bound states are eigenstates, scattering off N emitters simply maps c n (E) → t N E,n c n (E). As in the main text, one can write t N E,n = e iN φn(E) and Taylor expand φ n (E). Since we have found the pulse distortion terms decrease rapidly as n increases we simply consider the first order term in the expansion of φ n (E). Following the calculations starting from (S51) one obtains for zero detuning (k 0 = 0) which gives (17) from the main text. A similar process can be followed to obtain the normally ordered correlation functions.

S6. EFFECTIVE SPIN MODEL AND MPS ANSATZ
In this appendix we present an alternative route to the S-matrix scattering formalism used in the main text to describe the light propagation through an array of chirally coupled emitters. The full system dynamics can be indeed described by the combination of a driven-dissipative master equation for the emitters and a generalized input-output relation that allowed to reconstruct the electric field [24,25]. The dynamics of the emitter can than be efficiently solved by using a matrix product states (MPS) ansatz [25,28,29,40] which allow us to push the simulation to larger atomic arrays and higher input power.

A. Effective spin model
Let us consider the most general case (extension of the cascade scenario described by Eq.(1) of the main text) of N TLA with equal frequency ω 0 = ck 0 coupled asymmetrically to a waveguide with decay rates Γ R and Γ L associated respectively to the emission of right and left propagating photons. The addition of the two decay rates gives the total emission rate of the TLA into the waveguide Γ = Γ L + Γ R . In a realistic implementation the emitters can also radiate into other external channels, the decay rate associated to this emission is indicated with Γ 0 . The right propagating input pulse is treated as a classical coherent field, E in (t, x) = E in (t)e ikinx and for the rest of the discussion we will assume it to be on resonance with the atomic transition, i.e. k in = k 0 . The coupling of the emitters to the input field is given by the Hamiltonian H drive = − j √ Γ R E in (t, x j )σ + j + H.c. whereσ + j = |g e| j is the j-th emitter annihilation operator.
Under Born-Markov approximation the emitters dynamics is known to be described by a chiral master equation(ME) of the general form [27]:ρ where is the effective Hamiltonian which provides the non-Hermitian collective evolution of the emitters. The recycling term assures the conservation of the density operator trace and it arises by the quantum jumps that the emitters experience after the emission of a photon into the waveguide or into free space. For Γ L = Γ R Eq. (S72) coincides with the standard waveguide QED master equation [45,46] where both coherent and dissipative interactions between the emitters occur depending on the atoms position. In the limit of a perfect chiral waveguide considered in the main text, i.e. Γ L = Γ 0 = 0 and Γ = Γ R , the master equation Eq. (S72) reduces to the well known cascade master equation [26,27] characterized by an unidirectional purely dissipative interaction between the emitters. In this limit only the spatial order of the atoms and not the specific positions matters and the excitation is transferred only to the right without any back action. Once solved the emitters dynamics, the transmitted and reflected (for the semi-chiral case) electric field is obtained by the following generalized input-output relations [24,25]: The combination of Eq. (S72) and Eqs. (S75)(S76) allows to correctly describes the photon propagation trough the chiral medium.

B. MPS ansatz
In order to observe the main results obtained in the main text, e.g. photon number dependent bound states separation, it is necessary to consider substantial input pulse amplitudes and large atomic arrays, a scenario that is challenging to simulate with standard numerical techniques. To treat this many body problem we use a recently developed algorithm which involve an MPS representation to efficiently describe the effective spin dynamics [25]. The system evolution can be solved with two different approaches, either by directly solving the ME (S72) [40](method used for Fig. 3 of the main text) or by using a quantum trajectories algorithm where the state of the system evolves under the effective Hamiltonian (S73) and it stochastically experiences the occurring of quantum jumps [25] (method used for Fig. 2(f) of the main text). In both case an MPS representation is applied either to the quantum state or to the the density matrix once mapped to a vector. The MPS ansatz consists in reshaping the generic quantum state |φ/ρ = i1,..i N ψ i1,i2,..i N |i 1 , i 2 , ..i N into a matrix product state of the form: where, for each specific set of physical indices {ī 1 ,ī 2 , ..ī N }, the product of the Aī j matrices gives back the state coefficient ψī 1 ,ī2,..ī N . Each matrix Aī j has dimension D j−1 × D j known as bond dimension and finite-edge boundary conditions are assumed by imposing D 1 = 1 and D N = 1. The bond dimension reflects the entanglement entropy. For arbitrary states the bond dimensions grow exponentially with the system size. The advantage of the MPS ansatz it is that, in many physical scenarios, as the one considered here, the entanglement can grow slowly with the system size allowing an efficient description of the state in terms of a smaller bond dimension [29].
In order to compute the evolution of the system it is possible, in both methods, to derive a matrix product operator (MPO) representation for either the Liouvillian or the effective Hamiltonian and jump operators. The derivation of this representation for our semi-chiral case follows straightforwardly from the bi-directional case presented in [25,40]. The evolution of the system is then computed by applying a linear expansion of the master equation (or of the time evolution operator) 1 + dtL(H eff ) or a Runge-Kutta method. In addition a MPO representation can be derived also for the observables and it allows to evaluate the expectation values. After each application of an MPO to an MPS the MPS bond dimension increases. The bond dimension is than truncated after each step in order to keep an efficient description of the state. We indicate as D max the maximum bond dimension used to represent the system during an entire time evolution. In all the plots considered in this work we used D max = 150 for the ME simulations and D max = 40 for the quantum trajectories algorithm.

S7. EFFECT OF OTHER DECAY CHANNELS
In the main text we considered the ideal scenario of an array of atoms perfectly chirally coupled to a waveguide. In real implementations this is not always the case and the atoms can either emit into the other direction or into free space. Let us first assume that the atoms can emit also into the direction opposite to the pulse propagation, i.e. Γ L = 0, while the emission to other external channels is suppressed Γ 0 = 0. In this case the master equation (S72) becomes dependent on the atomic positions. To better understand this point we consider the two paradigmatic cases of atoms equally spaced either by a distance k 0 d = k 0 |x i − x j | = (n + 1)π/2 or k 0 d = k 0 |x i − x j | = nπ with n = 1, 2, ... We focus on these two configurations because in the bi-directional limit they are the ones providing purely coherent or dissipative interactions respectively. In Fig.S2(a)-(b) we plot the transmitted intensity for the two cases and we compare it to the one obtained in the perfectly chiral scenario. We observe that the emission into the opposite channel affects the intensity in a totally different manner depending on the atomic distance. In the first configuration, Fig.S2(a), it leads to an overall spreading of the bound states peaks for increasing values of Γ L . In the second, Fig.S2(b), the bound state separation not only seems to be robust to the deviation from a perfect chiral emission, but it even leads to a better resolution of the bound state peaks at the price of a weaker pulse in transmission. This effect is caused by an additional time-delay of the BS peaks, induced by partial reflection, which affects mostly the few photon bound states making the bound states separation even more evident. The many-photons BS, on the contrary, not only experience less delay than their few-photons counterpart, but they are also more robust against distortion tending to keep their solitonic behaviour. These features are more evident in Fig. S2(c) where we plot the reflected field intensity. Here the majority of the reflection occurs at the interface and it causes the damping in transmission. On the other end, when the pulse enter in the medium, we observe a major reflection for the BS with few photons. This phenomena make sense in light of the connection to the SIT limit at large photon pulses. In this limit, the chirality is not an essential requirement to observe the formation of solitons. The rule of chirality is instead crucial to bring the solitonic behaviour to the level of few photons.
In addition this observation suggests that the k 0 d = k 0 |x i − x j | = nπ spacing is the optimal configuration to observe the effect in an experiment. On the other hand, in many realistic implementation, it is difficult to have control on the position of the emitters. For this reason it makes sense consider the effect of a semi-chiral emission by making an average on random positions. This scenario is similar to consider a decay of the atoms in an external channel with rate Γ 0 . The effect of this external decay is plotted in Fig. S2(d) and shows that in general the effect is visible for values Γ 0 ≤ 0.2. Again we expect the overall effect to be robust for many-photon BS where the major requirement in this limit relies on the trapping of a 1D long array of atoms.