Simulating quantum light propagation through atomic ensembles using matrix product states

A powerful method to interface quantum light with matter is to propagate the light through an ensemble of atoms. Recently, a number of such interfaces have emerged, most prominently Rydberg ensembles, that enable strong nonlinear interactions between propagating photons. A largely open problem is whether these systems produce exotic many-body states of light and developing new tools to study propagation in the large photon number limit is highly desirable. Here we provide a method based on a “spin model” that maps quasi one-dimensional (1D) light propagation to the dynamics of an open 1D interacting spin system, where all photon correlations are obtained from those of the spins. The spin dynamics in turn are numerically solved using the toolbox of matrix product states. We apply this formalism to investigate vacuum induced transparency, wherein the different photon number components of a pulse propagate with number-dependent group velocity and separate at output. Numerical simulation of light propagation through 1D atomic systems in the many-body limit rapidly saturates hardware capabilities. Here, the authors tackle the problem by mapping the dynamics to an open 1D interacting spin system and solving it using the matrix product state ansatz.

A tomic ensembles are a very successful platform used to couple input light to atomic degrees of freedom, allowing the development of new quantum technologies. Historically, the weak optical nonlinearities associated with atomic ensembles have allowed most of the processes of interest, such as quantum memories for light 1 , to be describable within a limited realm of classical linear optics or Gaussian quantum states 2 . More recently, it has become possible to engineer strong interactions between photons in atomic ensembles and thereby realize highly non-Gaussian states. Under weak field inputs, for example, phenomena such as photon blockade 3 or two-photon bound states 4 in atomic Rydberg gases [5][6][7] have been experimentally demonstrated. A number of other systems, such as photonic waveguides coupled to atoms [8][9][10][11][12][13][14] or "artificial" atoms such as superconducting qubits and quantum dots [15][16][17][18] also show potential to realize similar physics.
A major question of interest is what occurs in such systems at higher field inputs. In particular, it is expected that strong interactions might lead to interesting many-body phenomena involving photons, such as photon crystallization (illustrated schematically in Fig. 1a). To address this question theoretically seems challenging: the systems are out of equilibrium, being driven by an external laser source; are open, where spontaneous decay of atoms leads to losses; and have long range interactions between atoms mediated by the exchange of photons. Some progress has been made in limiting regimes, where for example effective theories can emerge under certain approximations 11,12,[19][20][21][22][23][24][25][26] . While these effective theories provide useful insights, it would also be highly desirable to develop numerical methods with minimal approximations to verify these models and investigate regimes where the approximations break down, potentially revealing new physics as a result.
Currently numerical techniques are quite limited. For example, the standard approach is to describe the atom-light dynamics using Maxwell-Bloch equations 27,28 , which are solved by discretizing the atomic and photonic fields 3 , i.e., by modeling the continuum in which fields propagate by a finite number M of boxes, as depicted in Fig. 1b. To describe a state with up to n photons in each box then requires a Hilbert space with dimension of at least (n + 1) M (plus any atomic degrees of freedom associated with each box). In practice, this has limited numerical simulations to two 3,4,20,29,30 or three 23 total excitations over the entire system, and excluding the full effects of dissipation (in particular, either quantum jumps in wave function evolution, or population recycling terms in density matrix evolution).
Here, we show instead that simulating the high-photon number limit is possible by mapping the propagation problem onto the dynamics of a one-dimensional (1D) open "spin" model, which can be solved using the powerful matrix product state (MPS) ansatz 31,32 . The essence of the spin model is that in light propagation through an ensemble, the only independent degrees of freedom are the atomic internal states ("spins"), where the light fields mediate interactions between the atoms. A common theoretical approach is then to integrate out the fields, reducing the description to a problem of N interacting spins (N being the number of atoms) [33][34][35][36][37][38][39] . Furthermore, many light propagation experiments are quasi-one dimensional, where the input and output light are in a single transverse mode 2 . In this case, a simpler model capturing the same physics is a 1D spin chain coupled by the modes of a 1D photonic waveguide, as shown in Fig. 1c.
Our model further takes advantage of the fact that typically the number of spins in the chain does not appear independently, but rather macroscopic observables such as optical depth depend on the product of atom number and the atom-photon interaction probability. By tuning this probability to be large, a relatively small (~10 2 ) number of atoms is sufficient to model most atomic ensemble experiments. We can then use the MPS technique from condensed matter physics to numerically simulate the system, which is well adapted to treating chains of hundreds of spins. This technique depends on the fact that many of the quantum states we encounter in reality do not have large amounts of entanglement and are confined to a small portion of the in-principle exponential Hilbert space, allowing for a more efficient state representation. While the amount of entanglement present in atom-light interfaces is generally an unstudied problem, we give heuristic arguments below why we expect the MPS ansatz to work for such systems. As a benchmark of our model, we use it to simulate pulse propagation in the case of vacuum induced transparency 40,41 , which is one of the few cases where manyphoton propagation is qualitatively understood 30,42 .

Results
1D spin model of light propagation. While there are some phenomena in atomic ensembles that are truly three-dimensional, such as radiation trapping 43 and collective emission at high densities [44][45][46] , within the context of generating many-body states of light, the problems of interest largely involve quasi onedimensional propagation [19][20][21][22][23][24]26 . Indeed a typical experimental design is to input light in a single transverse mode and detect the light in the same mode after it traverses the ensemble. The standard approach to describe light propagation in such a system Fig. 1 Photon propagation in atomic ensembles. a Propagation through atomic ensembles with strong inter-atomic interactions can lead the quantum properties of the light to be strongly modified, represented here as crystalline-like output of photons from the medium. b In a typical approach to treat such quantum light propagation numerically, space is broken up into M discrete units centered at positions z j , with corresponding local Fock basis states 0 zj ; 1 zj ; for the photons. The equations of motion for the field may be solved on this space, where the growth of the Hilbert space with the number of Fock components on each site typically restricts calculations to a maximum of two or three photons in the entire system. c Here instead we model the quasi-1D propagation problem by a 1D waveguide coupled to N atoms. In this case the degrees of freedom associated with the light field are integrated out, to produce an effective model consisting of an interacting "spin" chain (with the spins being associated with atomic internal states). Output fields are then directly calculated from the spin dynamics using an input-output relation and the Hilbert space is reduced to the atomic one, which can be efficiently treated using matrix product states. As discussed in the Results section, the atoms are coupled to the waveguide modes with strength Γ 1D and for the simple case of two-level atoms the excited state |e〉 may also decay into modes outside of the waveguide at rate Γ′ is to use Maxwell-Bloch equations 27,28 in their one-dimensional, paraxial form 2-4, 20, 47-49 . There, the electric field operator can be decomposed as the sum of a forward propagating mode that travels in the direction of the input, taken here to be the positive z direction, and a backward propagating one, E(z, t) = E + (z, t) + E − (z, t). These have propagation equations, determined by the atomic polarization density operator P ge (z). For concreteness, we assume that the probe field couples to a single dipole-allowed transition from atomic ground state |g〉 to excited state |e〉, however, it is straightforward to modify these equations to account for additional atomic levels and transitions, driving fields, and interactions. The atomic polarization density, on the other hand, is driven by the field, and obeys an optical Bloch equation where ω eg is the atomic transition frequency, P ee,gg are the excited and ground state populations, and F describes the quantum noise associated with decay rate Γ′. Here we have introduced the coupling rate Γ 1D of an individual atom to the one-dimensional input mode. In principle, this rate can vary with z depending on the details of this mode, but for notational simplicity we assume here that it is constant. In this standard formulation of the Maxwell-Bloch equations, it should be noted that the interaction of the atoms with the remaining continuum of three-dimensional modes is reduced to an independent emission rate Γ′, meant to approximately capture scattering of photons out of the transverse mode of interest. The question of when this approximation breaks down is quite complicated and rich 20, 50, 51 and will not be discussed here; in any case, Eqs. (1,2) are widely accepted as the standard model for quasi-1D light propagation through atomic ensembles. Eqs. (1,2) represent an open, interacting quantum field theory, for which a general solution is unknown. The complexity is reduced in ensembles that lack strong non-linearities, where for example one can linearize the atomic system, such that the resulting joint quantum state of matter and light is Gaussian 2 . However, in the highly non-linear ensembles that are interesting for many-body physics we are typically restricted to solving numerically the Maxwell-Bloch equations by discretizing the fields, which, as mentioned in the introduction, is not feasible for more than a few photonic excitations.
Here instead, we take advantage of the fact that the Maxwell-Bloch equations presented above can also formally arise from a simple toy model of atoms coupled to a 1D waveguide 36,[52][53][54][55] . In particular, one can consider a model of atoms coupled to a bidirectional waveguide of infinite bandwidth, and frequencyindependent propagation speed c and coupling strength. In that case, the propagation equations of the forward and backward going modes are exactly those in Eq. (1), where P ge ðzÞ ¼ P N j¼1 σ j ge δ z À z j À Á for N atoms with positions z j . These equations can then be formally integrated giving a solution for the field as the sum of the input field E in ðz; tÞ and the field radiated by the atoms [36][37][38][39] , In the above equation the propagation of the field from the atomic position z j to z leads to a phase factor determined by the wavevector k 0 = ω eg /c. On the other hand, the time delay in freespace propagation is neglected and the field sees the response of an atom at another point instantly. That is, any pulse delay that arises is due to the atomic dispersion itself. This is justified when the free-space propagation time over the length of the system L is much smaller than the time scale characterising the atomic evolution, e.g., when L=c ( 1=Γ 0 , a condition easily satisfied in atomic ensemble experiments 36 . In limiting cases, this approximation can be further validated by solving for dynamics exactly and seeing that the results are the same 56 . Removing time retardation provides a drastic simplification of the problem as the equations of motion of the atoms and fields are now all local in time. Indeed, inserting Eq. (3) into the Heisenberg equation for the atomic coherences, one finds that the dynamics of the atoms depend only on the input field and on the state of the other atoms at the same time. Part of the system dynamics can then be derived from an atomic interaction Hamiltonian 36,54 , which describes the process of emission by an excited atom at z j into the waveguide, and the subsequent absorption of that photon by a ground-state atom at z l . These effective atom-atom interactions also lead to collective spontaneous emission, described in the master equation by the Lindblad operator In addition, we can add a phenomenological independent decay rate Γ′ as in the Maxwell-Bloch equations, which corresponds to scattering out of the quasi-1D input mode. This is described by the locally acting Lindblad operator The coupling of the atoms to the input field is given by þ H:c:Þ. In the following we will consider the case most common in experiments of a coherent state input, where E in can be treated as a classical field 57 and we neglect the associated quantum noise term, as this does not contribute to the normally ordered correlation functions of the output field that we will be interested in. The output field itself is the field measured past the position of the last atom, E out ðtÞ ¼ E z þ N ; t ð Þ given by Eq. (3), which is completely determined by the solution of the driven spin system and the input.
In the model above, the coupling of the atoms to the waveguide and the positions of the atoms must be chosen carefully to reproduce phenomena associated with free-space ensembles. In particular, as we discuss below, our numerical calculations are facilitated by choosing ratios of Γ 1D =Γ ′ $ 1. It is known that for a weak resonant input field, a single two-level atom can produce an appreciable reflectance of Γ 2 1D = Γ 1D þ Γ ′ À Á 2 52, 53 . The reflectance can be further enhanced if multiple atoms are placed on a lattice with lattice constant defined by k 0 a = π, in which case the reflectance from individual atoms constructively interferes 54,58,59 . While it is possible to observe similar effects in atomic ensembles 60,61 , this situation is atypical and will not be discussed further here. To reproduce the typical case in atomic ensembles where reflection is negligible, we choose a waveguide spacing of k 0 a = π/2, in which case reflection from different atoms in the lattice destructively interferes. In this configuration, the 1D waveguide model reproduces one of the key features of an atomic ensemble, that of decay of the transmitted field with increasing optical depth. If we consider the transmittance T ¼ hE y out E out i=jE in j 2 , then for a resonant weak coherent state input we find in the 1D waveguide model T = exp (−OD), where the optical depth is OD = 2NΓ 1D /Γ′ for Γ 1D ≲ Γ ′12 . Since OD ≲ 10 2 in realistic atomic ensembles of~10 6 weakly coupled atoms, by artificially choosing Γ 1D $ Γ ′ , the same optical depth is achieved with just tens or hundreds of atoms. At the same time, the essential properties of most quantum nonlinear optical phenomena are believed to depend only on optical depth 47,62,63 , or on a limited number of other parameters where atom number does not appear independently (such as the optical depth per blockaded region in a Rydberg gas 6,26,29 ). By matching these parameters using a much smaller number of atoms, we can then model the physics of interest in 3D atomic ensembles. The possibility that artificial effects (such as saturation) arise from low atom number can be eliminated by numerically checking that observables converge with increasing N while, e.g., decreasing Γ 1D in proportion to keep the key parameters fixed.
While our model can be used to reproduce the macroscopic observables of light propagation in a traditional atomic ensemble, we also note that it quantitatively captures the microscopic details of experiments where atoms or other quantum emitters couple to 1D channels. This includes atoms coupled to nano-fibers (Γ 1D =Γ ′ $ 0:05) 13 or photonic crystals (Γ 1D =Γ ′ $ 1) 14 , or "artificial" atoms such as superconducting qubits or quantum dots coupled to waveguides (Γ 1D =Γ ′ ) 1) [15][16][17][18] . In these cases, our model is valid when the spacing between the atoms is of the order of the wavelength of the light or larger, for smaller atomic distances additional effects can occur 51,64 .
Simulations using matrix product states. Using the 1D spin model described above significantly reduces the size of the Hilbert space required to simulate the light propagation problem, but the dimension still grows exponentially with atom number. This growth can be avoided in the case where the input field is sufficiently weak that the Hilbert space can be truncated to a maximum number of total excitations likely to be found in the system 12,36 . In the more general case, where many-photon effects are important, the full Hilbert space may be treated numerically for around 10-20 atoms depending on the size of the single-atom Hilbert space dimension d. Going beyond this requires some reduction of the Hilbert space and here we choose to use matrix product states, which have been successfully used in condensed matter to model a wide variety of 1D interacting spin systems 31,32 .
The key idea behind MPS is to write the quantum state of the spin chain in a local representation where only a tractable number of basis states from the full Hilbert space is retained. In the case of time evolution, these basis states are updated dynamically in order to have optimum overlap with the true state wave function. In particular, the wave function of a manybody system ψ j i ¼ P σ1; ;σN ψ σ 1 ;σ 2 ; ;σ N σ 1 ; σ 2 ; ; σ N j ican be represented by reshaping the N-dimensional tensor ψ σ1;σ2; ;σN into a matrix product state of the form where σ j represent the local d-dimensional Hilbert space of the atoms, e.g., σ j 2 e j i; g j i f gfor two-level atoms. Each site j in the spin chain has a corresponding set of d matrices, A σj , and by taking the product of these matrices for some combination of σ j 's we then recover the coefficient ψ σ 1 ;σ 2 ; ;σ N . The matrices have dimensions D j−1 × D j for the jth site (D 0 = D N+1 = 1), which are referred to as the bond dimensions of each matrix. We also define D = max j {D j } as the maximum bond dimension of the state ψ j i MPS . This representation is completely general, and as such the bond dimensions grow exponentially in size for arbitrary quantum states. In certain circumstances, however, the bond dimension D needed to approximate a state well might grow more slowly with N due to limited entanglement entropy, which enables MPS to serve as an efficient representation.
For example, this forms the underlying reason for the efficiency of density-matrix renormalization group algorithms for computing ground states of 1D systems with short-range interactions 65 . A priori, for our system involving the dynamics of an open system with long-range interactions, we know of no previous work that makes definitive statements about the scaling of D. However, we can provide some intuitive arguments that MPS should work well (at least without additional interactions added to the system). First, we note that although the dipole-dipole interaction term in Eq. (4) appears peculiar, being infinite-range and non-uniform, it conserves excitation number. For a single excitation, it simply encodes a (well-behaved) linear optical dispersion relation that propagates a pulse from one end of the atomic system to the other 12 , and thus does not add entanglement to the system. While the spin nature in principle makes the atoms nonlinear, thus far in atomic ensemble experiments the strength of nonlinearity arising purely from atomic saturation remains very small at the level of single photons, and thus one can hypothesize that only a small portion of the Hilbert space is explored.
Once extra interactions are added, at the moment the scaling of D must be investigated on a case-by-case basis. However, generically one expects that the system has a memory time corresponding roughly to the propagation time of a pulse through the length of the system. Thus, if the system is driven continuously, it should generally reach a steady state over this time and there will not be an indefinite growth of entanglement. In the case of pulsed input, the number of photons in the pulse limits entanglement and it is possible to establish upper bounds on the bond dimension required as we discuss in the Supplementary Note 1. For arbitrary n-photon wave functions the bound may still scale exponentially in the number of photons N n/2 , however in the case of vacuum induced transparency that we investigate as a benchmark the scaling is instead approximately quadratic in the number of photons.
In our MPS treatment of the spin model we adopt a quantum jump approach to model the time dynamics of the master equation 66 , which has been successfully applied to many-body dissipative systems 67,68 . As we describe in more detail in the Methods, this method decomposes the master equation evolution into an ensemble of quantum trajectories, which are formed by deterministically evolving pure states under an effective Hamiltonian H eff and stochastically applying quantum jumps to the system. As an aside, however, we note that MPS-based techniques for evolution of density matrices have also been developed [69][70][71][72] . Whether and when such techniques out-perform quantum jump methods for our problem is likely a subtle question, which will be explored in more detail in future work.
To numerically simulate the spin model there are then four essential manipulations of the MPS as illustrated in Fig. 2 and described in greater detail in the Methods. The first is (a) deterministic evolution of the MPS over a small discrete time step δt. Here an approximation of the time evolution operator e ÀiH eff δt % 1 À iH eff δt is applied to the MPS representation of the state. This is achieved by expressing the operator as a matrix product operator (MPO), a generalization of MPS to operators.
After this deterministic evolution the MPS then undergoes (b) stochastic quantum jumps that account for dissipation, realised again by applying MPOs to the MPS, where this time the MPOs correspond to the quantum jump operators O l . In each case, after applying an MPO (with corresponding bond dimension D W defined in the same way as for an MPS) to an MPS, the MPS bond dimension increases. The state must then be (c) compressed to constrain the growth of its representation in time. At any time step we may then (d) calculate observables, such as the output field, given the MPS representation of a state and the MPO corresponding to the observable. The steps are then repeated to obtain the full time evolution.
Vacuum induced transparency. The model introduced above gives a powerful and flexible algorithm for simulating the interaction of light with atomic ensembles in the multi-photon limit.
To demonstrate the utility of this approach we now investigate the phenomenon of vacuum induced transparency (VIT) 40 . This example also serves to benchmark our method, as exact solutions for non-trivial multiphoton behavior are not available, while in the case of VIT at least the qualitative nature of the system dynamics is understood.
VIT is closely related to the effect of electromagnetically induced transparency (EIT) 1 , which occurs in three-level atomic media. In a two-level medium incoming probe light that couples resonantly to an atomic transition |g〉-|e〉 is absorbed and scattered by the atoms into other directions, leading for example to the strong attenuation in the linear transmittance for high optical depth, T = exp(−OD). In EIT, an additional metastable level |s〉 (Fig. 3a) is also coupled to the excited state by a classical control field with Rabi frequency Ω, allowing probe photons to couple with spin-wave excitations from state |g〉 to |s〉, forming so-called "dark-state polaritons". The coupling to the spin wave leads to a strongly reduced group velocity relative to free space (v g = Ω 2 a/(2Γ 1D ) for a waveguide with spin chain lattice constant a 12 ), while the absence of population in |e〉 enables a pulse to propagate with minimal attenuation.
In VIT the control field is replaced by strong coupling of the atoms to a resonant cavity mode as shown in Fig. 3b 40,73 , which is described by the Hamiltonian H cav ¼ g P N j¼1 ðσ j es b þ H:c:Þ=2 in the case of uniform coupling g to a cavity mode with annihilation operator b. Here even when the cavity is empty the atomic medium can become transparent as vacuum Rabi oscillations transfer population from state |e〉 to |s〉 41  propagation of light in the system then takes on the nature of the non-linear coupling of the atoms to the cavity. Specifically, the formation of a spin wave from n probe photons is accompanied by the excitation of the same number of cavity photons, which produce an effective control field strength of ffiffiffi n p g. Since in EIT the group velocity of the light is determined by the control field, where v g ∝ Ω 2 , the group velocity in VIT becomes number dependent v n $ ng 2 a= 2Γ 1D ð Þ 30,42 . Fock states |n〉 input into the system are then expected to propagate at v n .
On the other hand, a coherent state |α〉 that has average number of photons |α| 2 is a superposition of Fock states, where n photons are present with probability e À α j j 2 α j j 2n =n!. Input into the VIT medium, these components are then expected to spatially separate due to their different propagation velocities, given sufficient optical depth. The output intensity can then be calculated naively by simply delaying the input Fock components by a time τ n = L/v n , where L is the length of the atomic medium. The output intensity in time resulting from such a toy model is shown in Fig. 3c, for a coherent state input pulse with average number 〈n pulse 〉 = 1. We have taken the system length to be L = 100a and the single photon velocity v 1 = 4aΓ′, which results for example from taking g = 4Γ′ and Γ 1D = 2Γ′ in which case the system's optical depth is 400. We note that the experimental conditions needed to observe photon number separation in VIT are difficult to achieve 41 , and thus our parameters are chosen to observe the desired effect, rather than correspond to a given experiment.
A plot similar to Fig. 3c was given in a previous theoretical treatment of VIT 42 , as at that time it was unknown how to calculate observables in the presence of losses and spatiotemporal effects, such as occurring from pulse entry and exit from the atomic medium. More recently, VIT has also been studied numerically in the weak-field limit using the space discretization technique schematically illustrated in Fig. 1b 30 . In the weak field limit, only the single photon manifold contributes to the output intensity and the higher number components are only visible in higher order correlation functions like g (2) . This also means that quantum jumps have a negligible effect on the system dynamics, and they were neglected in the calculations. In more general circumstances, using MPS simulations, we will show that the effects of quantum jumps and pulse distortion can have a significant effect on the output field.
For concreteness, we take input pulses with central frequency ω p and Gaussian envelope , which have an average photon number of n pulse ¼ α j j 2 $ 1. The average photon number chosen is not due to any intrinsic limitation coming from the MPS method itself, but rather because in VIT the spatial separation is largest for the Fock components with low photon number (Fig. 3c) and with |α| 2 = 1 the single photon and two photon components of the coherent state give an equal contribution to intensity emphasizing this effect. In this case, number states with three or more photons make up 8% of the input state and constitute 26% of the input intensity due to their high photon number.
To treat VIT, we include in the spin model formalism the atomic part H 0 of the total effective Hamiltonian (Eq. (8) in Methods section), Here Δ = ω p − ω eg is the detuning of the probe light from the |e〉-|g〉 transition frequency, δ c = ω p − ω c − ω sg is the VIT twophoton detuning and κ is the decay rate of the cavity mode. In what follows we assume both the probe and cavity are resonant with their respective transitions, so that Δ = δ c = 0. Dissipation via the various loss channels is then included through quantum jump operators, where the collective emission into the waveguide is described by O ± as detailed in the Methods. The jump operator corresponding to cavity decay is O c ¼ ffiffi ffi κ p b and we assume that the atomic excited state can decay via free-space spontaneous emission into either state |g〉 or |s〉 (taking these decay rates to be equal for simplicity), leading to 2N jump operators O j;ge ¼ ffiffiffiffiffiffiffiffiffi The cavity mode is represented in our MPS treatment by an additional site in our spin chain, which can support up to n c bosonic excitations. In the simulations we present here we have taken n c = 10 and observe no difference in observables if n c is increased.
In Fig. 4a, b, we show the time-dependent output pulse intensity I out ðtÞ ¼ hE y out ðtÞE out ðtÞi calculated from an MPS simulation of 100 atoms and an input pulse with |α| 2 = 1. We also show the zero-delay second-order correlation function I ð2Þ out ðt; tÞ ¼ hE y out ðtÞE y out ðtÞE out ðtÞE out ðtÞi. In the output intensity two main peaks are observed, where the first peak in time (tΓ ′ $ 23) is due to photon number components with two or more photons, while the last peak (tΓ ′ $ 36) is associated with the slow propagation and exit of the single-photon component. That the most delayed part contains only single photons can be confirmed by looking at the second order correlation function which is only non-zero in the first part of the pulse. In Fig. 4b, we see good agreement between the features of the numerical pulse shape and the expected group velocity for each part of the pulse (compare with Fig. 3c), where the vertical black dashed lines represent the expected times for the peaks of the Fock state components, that is, with delays τ n .
Compared with the ideal picture in Fig. 3c, where a clean separation is seen between one and two photons, one can see that the full simulation produces a much larger intensity between the one-and two-photon peaks. We now show how the trajectories from the MPS simulations can be further filtered and analyzed, to gain insight about the underlying physics. In particular, we find that quantum jumps play a key role in blurring the separation between the different number components in the output, even for the very good system parameters that we have chosen (OD = 400, g=κ $ 130). An intuitive picture of how the blurring occurs can be gained by considering two photons that enter the medium, and initially propagate at a velocity v 2 = 2v 1 . During evolution, this state may decay via spontaneous emission into free space and leave behind a single photon propagating in the medium, at which point the group velocity is slowed to v 1 . This change in group velocity can happen at any point in the system and leads to single photons that arrive at the output earlier than expected if just a single-photon Fock state was input into the system, destroying the perfect separation of the single photon output from the two photon component. We can quantify this behavior by analyzing the quantum jumps that happen in our simulations, where due to the choice of physical jump operators discussed in the Methods, the total number of jumps in a given trajectory corresponds to the number of photons emitted from the system. Furthermore, the type of jumps (and thus the emission channel) can be explicitly tracked, between free-space loss, cavity loss, or detection in the waveguide output. In Fig. 5a, we show a histogram of the average number of jumps into the output waveguide channel versus time for the 20000 trajectories used to produce Fig. 4. This provides an alternative way (compare to Fig. 4) to calculate the intensity, as would be done in an experiment where detector counts are averaged over many identical realizations.
Moreover, we can classify the jumps according to whether they come from trajectories where 1, 2 or 3+ photons are emitted into the waveguide (as indicated by the different bar colors in Fig. 5a). As we see in the plot, the higher the number detected in the waveguide, the earlier in time the jumps happen, in agreement with the simple theoretical model and with the calculations of I out (t) and I ð2Þ out ðt; tÞ, discussed above. We can also select only the trajectories where a single photon is detected at the waveguide output, and further separate those trajectories into two distinct cases: (i) when that is the only jump event (indicating a single photon was input and successfully propagated through the system), and (ii) where a multi-photon state was input, and all but one photon decayed into other channels. The histogram according to this classification in time is shown in Fig. 5b, where we see that the tail of faster arriving single photons, seen to the left of the main peak, results from the decay of number states with two of more photons, and the resulting mixing of propagation velocities.
Alternatively, we can use the jump statistics from a coherent state input to identify the intensity resulting from a Fock state input. Since the VIT system does not support any long lived excitations (compared with the simulated time scale), the total number of photon jumps (into any channel) out of the system for any one trajectory is equal to the number of the photons that entered the system for that trajectory. By post-selection on the total number of jumps we can then find the intensity that results from a Fock state input as shown in Fig. 5c. Here we see the same effect of jumps as noted above but observed in a different way. In particular, while we categorized the trajectories in Fig. 5a, b by the number of photons that survive and are output, in Fig. 5c we classify them by the number that are input. For Fock state inputs of two or more photons, the output intensities show tails of longer than expected delay times, again as a result of photon loss and the mixing of propagation speeds.
These altered delay times are not only due to quantum jumps however, they can also result from distortion of the multi-photon wavepacket as it enters the medium 30 . This distortion happens as the input pulse crosses the boundary of the atomic ensemble, as we illustrate for a two-photon wave function in Fig. 6a. For example, if the two photon wave function has a Gaussian pulse shape, the two photons can arrive at the boundary of the atomic ensemble at different times. The first photon that enters then travels at v 1 until the time that the second photon enters and the group velocity becomes 2v 1 . A similar process occurs when the photons exit the medium. In this case the further the photons are separated in the original pulse, the larger the delay of the photons. This process distorts the two-photon input Gaussian into a heart shaped output and higher photon number manifolds into higher dimensional hearts. In Fig. 6b we show how this behavior can be observed in the two time correlation measurement of the output photons for an input coherent pulse at low input photon number. For higher photon number input the heart shape is distorted as higher photon number manifolds with larger group velocity smear out the distribution.

Discussion
In summary, we have described a novel technique to numerically simulate quasi-1D quantum light propagation through atomic ensembles, which is based on the powerful toolbox of matrix product states. This technique is versatile and adaptable to many cases of theoretical and experimental interest (e.g., with regard to  Fig. 4. b Stacked bar graph for quantum jumps from trajectories where only a single photon is detected in the output of the waveguide. These jumps are then divided into jumps that are not accompanied by any other jump into other channels, and those that are. We see that the tail of photons detected earlier are due to trajectories where 2 or more photons entered the medium but all but one were lost into other channels. c Post-selection of trajectories to find evolution for Fock state input. By selecting only trajectories where there were a total of 1 (blue), 2 (red), or 3 (yellow) jumps into any channel we can reconstruct the intensity output for the corresponding input Fock state level structure, types of interactions, additional degrees of freedom, etc.). Similar to the important role that DMRG and MPS played in one-dimensional condensed matter systems, we envision that results gained from our numerical technique could be used to push forward the development of effective theories of strongly interacting systems of light 11,12,19,[22][23][24][25] , and conversely that such analytical work could be used to improve numerical algorithms. Beyond that, it would be also interesting to investigate further why MPS apparently works well in the context of our open, long-range interacting system, and under what conditions MPS might fail. This could yield a better understanding of the growth of entanglement, which naively seems like a potentially useful resource, but which has not been explored for such systems to our knowledge.
The ability to formally map atom-light interactions to quantum spin models is intriguing in general, and it would be valuable to explore whether other techniques for solving spin systems give further insights into atom-light interactions. Finally, it should be noted that this mapping essentially relies on the fact that the atomic response dominates the dispersion of near-resonant light fields, as compared to the dispersion of empty space. It would thus be interesting to investigate whether a similar effective theory could be derived for other strongly dispersive systems, such as exciton-polariton condensates [74][75][76] , to shed new light on interacting photon dynamics in those settings.

Methods
Quantum jump formalism. To find the time dynamics for the spin model we must evolve the master equation in time. Numerically this can be done directly by evolving the full density matrix ρ in time using standard techniques such as the Runge-Kutta algorithm. Alternatively, we can instead use the "quantum jump" approach to unravel the master equation into trajectories of evolving pure states 66,68 . Here we briefly review the quantum jump formalism, which we implement with MPS as discussed below.
We write the master equation for our 1D spin model in the form _ ρ ¼ ÀiðH eff ρ À ρH y eff Þ þ P l O l ρO y l , where O l are the "jump" operators associated with the dissipation resulting from emission into the waveguide and into free space, and H eff is a non-Hermitian effective Hamiltonian. This division of the master equation into jump terms and an effective Hamiltonian H eff is not unique and we attempt to do so here in a way that the jump operators have a physical significance. In particular, the emission of a photon into the forward going mode of the waveguide may interfere with the input light that is also traveling in the positive z direction (Eq. (3)), an interference that would be present in real detection of photons output from the waveguide. This interference can be taken into account in our jump operator, and as such we take the forward going jump operator to be ge as in ref. 36 ). The backward going jump operator is simpler given the lack of input field in that mode, j ge corresponding to the free space decay, giving a set of possible jumps With the jumps formulated in this way the effective Hamiltonian becomes In general H 0 can describe any additional atomic evolution; in the specific case of two level atoms coupled to a probe of frequency ω p we can write, in the frame rotating with in the input frequency, The quantum jump approach uses the above decomposition of the master equation to restate the evolution of the density operator as a sum of pure state evolutions called trajectories 66 , where the wave function evolution is divided into (a) deterministic evolution under H eff and (b) stochastic quantum jumps made by applying jump operators O l . Starting from a pure state |ψ(t)〉 at time t, the deterministic evolution over a time step δt gives ψðt þ δtÞ j i¼ e ÀiHeff δt ψðtÞ j i. However, during this evolution the norm of the state decreases to δp ¼ 1 À hψðtÞje iH y eff δt e ÀiHeff δt jψðtÞi, as the effect of the jump operators is neglected. The effect of these operators is instead accounted for stochastically, where after each deterministic evolution we generate a random number r between 0 and 1. If r > δp the system remains in state |ψ(t + δt)〉. Otherwise, the state makes a random quantum jump to ψðt þ δtÞ j i¼ O l ψðtÞ j iwith probability δp l ¼ δt hψðtÞjO y l O l jψðtÞi. The state is then normalized and the process repeats for the next time step and each sequence of evolutions gives a quantum trajectory. Any observable can be obtained by averaging its value over many trajectories. Furthermore, as we choose our quantum jumps to relate to physical processes, the distribution of the jumps can be thought of as corresponding to actual photon detection in an experiment.
Time evolution with MPS. As discussed in the Results section and depicted in Fig. 2, to evolve the MPS of the spin model in time and measure the expectation value of different observables we perform the following four steps.
(a) Deterministic time evolution. To evolve the state |ψ(t)〉 in time we need to apply the operator e ÀiHeff δt to the MPS representation. This is achieved by applying a matrix product operator (MPO) to the state, where just as a state can be decomposed into an MPS, any operator W can be expressed in a local representation as Γ′ 2 Fig. 6 VIT pulse distortion. a Distortion process as a two-photon wavefunction ψ(z 1 , z 2 ) enters the atomic medium. The initial Gaussian distribution of the photon positions z 1 and z 2 is shown as a circle. The twodimensional space of the photon pair is divided into regions where only one photon is inside the medium and has group velocity v 1 , indicated by the dashed lines, and when both photons are inside the medium having velocity v 2 = 2v 1 , the square box. Photon pairs with greater separation spend more time in the regions where only one photon is inside the medium, delaying them compared to pairs with z 1 = z 2 , leading to a characteristic heart shaped pattern. b Two-time correlation function I ð2Þ out t 1 ; t 2 ð Þ for the output field of the VIT system, after excitation with a coherent Gaussian input pulse for various average input photon number: |α| 2 = 0.01, 0.25, 1.0, and 2.0. At weak input field the correlation function is purely due to the two photon component and shows a clean heart shape. As the number of input photons increases, higher photon number components contribute, which travel faster through the medium distorting the pattern and pulling it forward in time. The system parameters were for an optical depth of OD = 60, with N = 30, Γ 1D = Γ′, Δ = δ c = 0, g = 4, Γ′, κ = 0.03Γ′, σ t = 4/Γ′ and t 0 = 6/Γ′. We used bond dimension D = 30 and time step δt = 0.01/Γ′ where convergence was observed for all observables of interest (see Methods) Here W σ ′ j ;σj are a set of matrices at site j, where the matrices now have two physical indices σ ′ j ; σ j due to W being an operator. An MPO may be "applied" to an MPS via a tensor contraction over the physical indices σ j of the MPS and MPO, as shown in Fig. 2a. This generates a new MPS with higher bond dimension, as the bond dimension of the MPO, D W , multiplies the bond dimension of the original MPS, and for the calculation to be tractable D W must be small. Such a compact form is not known for the operator e ÀiHeff δt ; however, the first order approximation e ÀiHeff δt % 1 À iH eff δt has a compact MPO form if H eff does. This is the case for the 1D spin model where the MPO representation of the effective Hamiltonian has D W = 4. We can write σ j are matrices of operators given by for 1 < j < N, and end vectors and Here λ ¼ e ik0 a , I j is the identity operator for site j and the H j loc contain all the local terms in H eff . The MPO of the linear expansion of the time evolution operator 1 − iH eff δt can be obtained from this MPO without increasing the bond dimension. It is enough to replace W 1 with W t:e: to obtain the desired MPO. Using a small time step δt we can then advance the wave function in time by applying this MPO. (b) Quantum jumps. After evolving a time δt, the state is either kept and renormalized, or a jump is applied. To apply the quantum jump formalism we then just require an MPO form of the jump operators that can be applied to the MPS at each time step, see Fig. 2b. The jump operators of the 1D spin model can be written in compact MPO form, where loss into the waveguide requires an MPO of bond dimension D W = 2. The jump operator corresponding to the emission of a photon in the +z output channel can be written as O + = Z 1 Z 2 .. Z N with for 1 < j < N, and end vectors and The MPO of O − is analogous, but without the external field term in Z 1 and with k 0 replaced by −k 0 . For the jumps associated with spontaneous emission into free space an MPO representation is not required as these jumps just require an operator to be applied locally to a single site. (c) State compression. After applying the time evolution operator or jump operators the size of the MPS increases as the bond dimension of the operator multiplies the bond dimension of the original state. Over time this would lead to exponential growth in the MPS size if not constrained. This increase in bond dimension can correspond to the true build up of entanglement, but may also correspond to the new state being expressed inefficiently in the MPS form. In the second case, a more efficient representation can be found and the bond dimension compressed to a smaller value, as in Fig. 2c. This can be done using singular value decompositions to find low rank approximations of the matrices A σj in the MPS representation, or by variationally exploring the space of MPS states with a fixed bond dimension that are closest to the original state 31,32 . The validity of such a compression can be evaluated by checking how strongly the parts of the state discarded in the compression contribute to the description. From this an error can be calculated and the bond dimension in the compression adjusted so the error remains small (see below). (d) Calculating observables. At any point in time observables such as the spin populations or output field may be calculated for a particular quantum trajectory by applying the appropriate operator associated with that observable in MPO form to the state. For example, to find the output intensity, hψðtÞjE y out ðtÞE out ðtÞjψðtÞi, one can express the individual elements as matrix product states or operators. The intensity for that trajectory can then be evaluated through a tensor contraction, as shown in Fig. 2d. This intensity is then averaged over all the quantum trajectories to find the expectation value I out ðtÞ ¼ hE y out ðtÞE out ðtÞi. Multi-time correlation functions such as I ð2Þ out ðt; t þ τÞ ¼ hE y out ðtÞE y out ðt þ τÞE out ðt þ τÞE out ðtÞi can also be found. This is done by propagating the state in time until time t and then applying the operator E out to the state. The state is evolved a further time τ and the operator applied again. The norm of the resulting states are then averaged over many such evolutions to find the two-time correlation.
VIT matrix product operators. The MPO of the VIT Hamiltonian can be obtained by extending the bare spin model case above. The cavity degree of freedom is associated with an additional site in the spin chain at position N + 1, in which case the VIT MPO of H eff is obtained by adding two columns and rows to the bare representation: however, for the results of the simulation to match reality, this compression must be done in a controlled manner to avoid discarding important information from the state. One straightforward way to check the validity of the simulations is then to do the same simulation for various maximum bond dimensions D and see that values of the observables of interest converge as D increases. In Fig. 7a, b we plot the output intensity and zero-delay second-order correlation function for different bond dimensions for the trajectory without jumps (similar results hold for the other trajectories) for VIT with parameters as given in Fig. 4. We see that the intensity has an excellent convergence already for D = 20, while the zero-delay second-order correlation function requires higher bond dimension D~50 for convergence. From this we conclude that to accurately model the higher number components of the pulse requires larger bond dimension, as the higher number components of the pulse have higher weight in the second order correlation function. Correspondingly, in our simulations we fix the bond dimension to D = 50.
Another way to ensure the simulations accurately model the physical reality is to monitor the error incurred in each compression step. The compression can be done variationally minimizing the distance between the larger MPS ψðtÞ j i D′ and its compressed version ψðtÞ j i D , or by a sequence of singular value decompositions of the bond connections between each site 32 . In the latter case, the large bond dimension D′ yields D′ singular values λ t,j,l at bond site j and time t, with 1 ≤ l ≤ D′. If we order the singular values to monotonically decrease with increasing l, we may then reduce the bond dimension by keeping only the singular values with l ≤ D. One measure of this compression error is the norm of the difference of the original state and the compressed state ϵ t ¼ ψðtÞ j i D À ψðtÞ j i D′ , which can be expressed in terms of the discarded singular values, where ϵ t ¼ 1 À Q NÀ1 j¼1 1 À ϵ t;j À Á with ϵ t;j ¼ P l>D λ 2 t;j;l . The error accumulated during the whole time evolution is ϵ tot ¼ 1 À Q t 1 À ϵ t ð Þ . Since all the terms are small one can approximate the products with sums and obtain giving a figure of merit for the overall quality of the time evolution. In Fig. 7c, we plot the accumulated error for different bond dimensions. We find that starting from D = 20 the accumulated error decays as a power law ϵ tot $ D Àα , with α ≈ 2.9.
Data and code availability. The data presented in the figures of this manuscript are available from the corresponding author upon request. The code used for the MPS simulation of the spin model is available at https://github.com/jdnz/ MatrixProductStates.jl.