Optical properties of an atomic ensemble coupled to a band edge of a photonic crystal waveguide

We study the optical properties of an ensemble of two-level atoms coupled to a 1D photonic crystal waveguide (PCW), which mediates long-range coherent dipole-dipole interactions between the atoms. We show that the long-range interactions can dramatically alter the linear and nonlinear optical behavior, as compared to a typical atomic ensemble. In particular, in the linear regime, we find that the transmission spectrum reveals multiple transmission dips, whose properties we show how to characterize. In the many-photon regime the system response can be highly non-linear, and under certain circumstances the ensemble can behave like a single two-level system, which is only capable of absorbing and emitting a single excitation at a time. Our results are of direct relevance to atom-PCW experiments that should soon be realizable.


I. INTRODUCTION
Photonic crystal waveguides (PCWs) -dielectric media with a periodically varying refractive index -have attracted significant interest recently as a platform for realizing novel quantum light-matter interfaces. The ability to engineer PCW properties permits control of the electromagnetic environment experienced by nearby atoms, which may be leveraged to achieve strongly enhanced atom-photon coupling efficiencies [1][2][3][4][5] compared to more traditional approaches [6][7][8], as well as for the exploration of new regimes of quantum optics.
An exciting example of the latter [9][10][11][12][13][14][15] is the ability to engineer long-range coherent interactions between atoms. The emission properties of an atom coupled to a PCW differ markedly from those in free space when the resonant transition frequency of the atom, ω a , is inside a band gap, i.e. a frequency domain where guided modes are absent, as depicted in Fig. 1a for the blue-colored band. Rather than radiating away, the field emitted into this set of modes is localized around the atom by Bragg reflection, giving rise to an atom-photon bound state [16][17][18], shown schematically in Fig. 1b. The photonic component of this state constitutes an effective cavity mode that can mediate coupling to other atoms, resulting in effective dipole-dipole interactions which acquire the spatial features of the cavity mode. These interactions can in turn be exploited to realize novel regimes of quantum many-body physics and nonlinear optics [12,14,19].
Motivated by remarkable advances to interface atoms or solid-state emitters with PCWs [2,4,15,20], we study the fundamental optical properties of a system of two-level atoms with long-range interactions, finding rich behavior that differs markedly from an ensemble of independent atoms. For example, as is evident from Fig. 1c, even the linear transmission spectra can be dramatically modified as compared to a typical atomic ensemble. We will describe the origin of the multiple transmission dips, and show that the resonances associated with some of them exhibit large optical nonlinearities that yield strong photon anti-bunching in the scattered light. Our results are a step towards developing a comprehensive understanding of the linear and nonlinear optical properties of such systems of interacting atoms.

II. MODEL AND METHODS
We model the simple case of an ensemble of two-level atoms with ground and excited states |g , |e , coupled to two guided modes of a 1D PCW, with dispersion relations (frequency ω vs. Bloch wavevector k) as depicted schematically in Fig. 1a. A realistic platform for implementing such a model is the 'alligator' PCW described in Ref. [2]. The atomic transition frequency, ω a , is in a band gap of one set of modes, which in practice may have, e.g., transverse electric (TE) polarization. As described above, the absence of propagating states at ω a prevents an excited atom from radiating into the TE modes. However, an exponentially decaying photonic cloud of length L can form around the atom, facilitating coherent excitation exchange with a proximal atom in its ground state with a strength V . The quantities V, L depend on the curvature of the dispersion relation at the lower band edge and the frequency separation δ between the band edge and atomic transition, as discussed further in Ref. [13]. We assume the detuning between ω a and the upper TE band edge to be sufficiently large that this band may be neglected.
While the above mechanism can produce novel long-range interactions, it does not allow the atoms to be efficiently probed using the same TE modes, since it is not possible to send in propagating fields at frequency ω a . To this end, a second band of modes, which may correspond to transverse magnetic (TM) polarization, is used as a conventional waveguide to excite the atoms and collect the resulting scattered fields. Our model is then a minimal description that produces long-range interactions, and allows their effect to be probed. In particular, the level structure is considerably simplified as compared to previous schemes to realize optical nonlinearities in PCWs [12,14]. The coupling strength of a single atom to the TM band is characterized by an emission rate Γ 1D , and the associated wavevector k a is shown in Fig. 1a. We assume the atoms are trapped along the axis of the crystal in a lattice of period d, equal to the period of the crystal itself. The dynamics of the full system may be described by an effective non-Hermitian Hamiltonian for the atoms alone [13,21]: Here ∆ = ω L −ω a denotes the detuning of the probe field (incident from the left with wavevector k L , driving frequency ω L , and Rabi frequency Ω) from the atomic resonance frequency, Γ ′ is the decay rate into free space and all other modes, and z j is the position of the j th atom along the waveguide. In all simulations we set Γ 1D /Γ ′ = 0.3, however we will also provide more general conclusions independent of the specific parameter choices. The operator σ j µν = |µ j ν j |, where µ, ν = g, e are energy eigenstates of the atom j. We assume that the atoms are trapped at the anti-nodes, where the interaction with the TE modes is maximized, and results in an integer phase factor θ jk = (z j + z k )/d (as z j is an integer multiple of d).
To obtain the fields scattered into the waveguide, we use the generalized input-output methods described in Ref. [22]. To preserve the norm of the wave function, dynamics under the non-Hermitian Hamiltonian of Eq. (1) must in principle be supplemented with quantum jumps; however, we will primarily be interested in field correlations under weak driving (Ω/Γ ′ ≪ 1), in which case it can be shown that quantum jumps may be neglected. We time-evolve the initial atomic wavevector |g ⊗n under Eq. (1); the transmitted (T) and reflected (R) fields are then reconstructed using the input-output relations where the transmitted (reflected) field is evaluated at a position z beyond the last (first) atom of the array (henceforth we drop the subscript out and the argument z). Physically, Eqs. (2)-(3) state that the output field properties are completely encoded in the input field and the correlations of the atomic scatterers. The vacuum input noise operators F T,R may be neglected on our correlations of interest. The steady-state transmittance and two-photon correlation function of the transmitted field are given by and similarly for the reflected field, where |ψ is the steady-state wavevector. For stronger driving fields, we solve the full atomic density matrix equations, thus accounting for quantum jumps. Fig. 1c shows a representative single-shot (atoms fixed in position) linear transmission spectrum for one configuration of n = 10 atoms randomly placed in a lattice of N = 200 sites. We choose the lattice constant d such that k a d = π/2, which for unit filling factors (i.e. for n/N = 1) minimizes back-reflection from the array, and thus most closely corresponds to a free-space atomic ensemble. However, different choices of k a d (excluding those very close to integer multiples of π) do not qualitatively change the results.

III. LINEAR OPTICS
For the case V = 0 (red curve in Fig. 1c), the transmission spectrum is identical to that of a free-space ensemble, i.e. a broadened Lorentzian centered at the single-atom resonance frequency. The resonant (∆ = 0) transmittance is given by T ∼ e −D , where D ≈ 2nΓ 1D /Γ ′ is the optical depth. The blue curve shows the spectrum in the presence of the bandgap (BG) interaction, with strength V /Γ ′ = 4 and length L/d = 100. Clearly, the system now has new resonance frequencies that may be resolved spectroscopically.
To understand the structure of the linear spectrum, note that in the limit L → ∞, the BG interaction term in Eq. can be lifted to yield additional observable resonances, as shown in Fig. 1c, while the maximum resonance energy, ω max , is correspondingly reduced. However, as we explain below, the properties of this maximum resonance can still be well quantified.
Current techniques for interfacing atoms with PCWs do not allow precise control of the atomic positions z j , and moreover, due to finite trap lifetimes, single-shot spectroscopy is not possible. We therefore seek to understand the average optical properties, taking (unless stated otherwise) the mean of the appropriate quantities over 1000 random spatial configurations. Fig. 1d shows the average linear transmission spectra for the same conditions as Fig. 1c. While each realization has a unique value of ω max , as we explain in Appendices A 1 & A 2, for L/d ≫ 1 the average ω max scales linearly with the number of atoms in the system, as shown in Fig. 2a, with a slope that is dependent on the ratio L/N d. Indeed, fixing N d and defining an effective interaction strength V eff (L) such that ω max (L) = V + (n− 1)V eff (L), a simple estimate gives (see Appendix A 1) As shown in Fig. 2b, this simple model agrees well with the numerically obtained solution. As we explain further in Appendix A 1, when n ≪ N , V eff (L) may then be used to determine the frequency ω max for a PCW with known V , L, N , and d. It follows that the number of atoms may be inferred from the frequency of the highest-energy transmission minima. This may be a useful experimental characterization tool, given that the resonant optical depth ceases to be directly related to the total atom number in the presence of BG interactions.
Another important feature of the maximum resonance ω max to characterize is the associated transmittance T dip , as illustrated in Fig. 1d. Intuitively, this should be a function of the probability of creating the collective excitation associated with the maximum eigenvalue. In the limit L → ∞, it may easily be shown that this excitation has the form We can then quantify the probability of exciting this state by considering its overlap with the initial excitation |ψ in provided by the input field, where |ψ in = (1/ √ n) j e ikLzj |e j . Fig. 2c shows (as we explain fully in Appendix A 3) how, as a function of the number of atoms n, the mean overlap φ max |ψ in decays approximately as 1/ √ n for moderate filling factors, and tends to zero for very large filling factors. We find that φ max |ψ in ∼ √ κ n , where κ n = (N −n)/( √ 2nN ), gives a very good estimate of the overlap in all regimes.
As a result of the decreasing overlap, the transmittance T dip increases as a function of n, as shown in Fig. 2d (dots) for two different values of L. As explained fully in Appendix A 3, we can also construct an effective model that treats the system as having only two degrees of freedom -an 'ensemble' level |E with the bare single-atom resonance frequency, and the maximum resonance |1 at frequency ω max , with coupling strengths and decay rates as shown in the level diagram of Fig. 3a. The crosses in Fig. 2d show the effective model predictions for T dip for the two values of L, agreeing well with the full model in both cases. Moreover, for the case L → ∞ the effective model may be solved exactly to give T dip ∼ Γ ′ 2 /(Γ ′ + nκ n Γ 1D ) 2 , which is analogous to the resonant transmission T = Γ ′ 2 /(Γ ′ + Γ 1D ) 2 through a single atom coupled to the probe band [21].

IV. NONLINEAR OPTICS
We now turn to the nonlinear optical properties of the system, beginning with the second-order correlation function of the reflected field, g R , when the system is driven at its maximum single-excitation energy, ω max . The reflected field arises purely due to scattering from the ensemble, whereas the transmitted field is an interference between the scattered and incident fields, see. Eq. (2). Thus, g R yields information about the ability of the system to scatter two photons simultaneously. The limiting case of a single two-level atom produces perfect anti-bunching, g (2) (0) = 0, indicating that it can only absorb and re-emit single excitations at a time.
For a given value of V , we compute the correlation function g   Fig. 4a. Evidently, the stronger the interaction, the more prevalent strong anti-bunching effects become. In particular, as Fig. 4b shows, for an interaction strength V /Γ ′ = 6, the proportion of configurations with g R (0) < 0.1 approaches 90%, suggesting that in the vicinity of the maximum resonance ω max , the ensemble typically behaves as an effective two-level system.
The origin of the anti-bunching in the L → ∞ limit may easily be understood: by diagonalization of the bandgap interaction (proportional to V in Eq. (1)), we find the maximum single-photon eigenenergy to be ω (1) max = nV , and the maximum two-photon eigenenergy to be ω Thus, it is clear that for large V the system approaches ideal two-level behavior. For finite L, provided that the linear spectrum is well characterized by V eff (L) (Eq. 6), the anharmonicity A(L) ∼ −2V eff (L)see Appendix B 2.
Motivated by the above observation, an interesting question is whether a global Rabi pulse applied to all of the atoms can produce only a single collective excitation in the ensemble (as well as the associated fidelity as a function of the various system and laser parameters). To answer this, we study the dynamics under Eq. (1) with a driving field Ω that is sufficiently large to produce Rabi oscillations, solving the full master equation for the atomic density matrix when the system is initialized in its ground state |g ⊗n . Fig. 5a shows a representative sample of the time evolution of the single-excitation manifold population, p 1 , for a system of 6 atoms driven at its maximum resonance frequency and for different driving strengths, with V /Γ ′ = 10. Clearly, depending on the driving strength, p 1 can be very large. Intuitively, to maximize p 1 , Ω must be sufficiently large to overcome dissipation, yet sufficiently small to minimize subsequent population of the doubly-excited manifold: we observe this trade-off in Fig. 5a, where the choice Ω/Γ ′ = 5 gives a larger maximum value of p 1 than the cases Ω/Γ ′ = 1 and Ω/Γ ′ = 10. Fig. 5b shows the maximum value of p 1 during the evolution time, as a function of both Ω and the detuning ∆ max = ω L − ω max from the maximum resonance frequency. Since the anharmonicity is negative, to minimize the probability of two-photon absorption the optimal choice of ∆ max is positive. The optimal single-excitation probability p opt 1 ≈ 0.81 occurs for the parameters Ω opt /Γ ′ = 6.75, ∆ opt max /Γ ′ = 0.4. While Figs. 5a,b are calculated with full density matrix dynamics, we now introduce a simplified effective model, building on the linear model discussed above. This enables both a simple understanding of the dependence of the maximum single excitation probability p 1 versus system parameters, and predictions when the atom number becomes too large for exact density matrix solutions to be feasible. As depicted in Fig. 3b, the model treats the system as having two separate 'branches', corresponding to the ensemble of independent atoms (with m th excited state |E m ), and the spin-wave excitation of the maximum-energy resonances in the one-and two-photon manifolds, denoted by states |1 and |2 . Fig. 5c shows the effective model prediction for p 1 as a function of the external field properties, under the same conditions as for Fig. 5b. While the optimal laser parameters (Ω opt /Γ ′ = 6.75 and ∆ opt max /Γ ′ = 2.4) do not match those of the full model, the optimal single-excitation probability is accurate to within ∼ 1% (p opt 1 ≈ 0.80). We may then consider the fidelity of introducing only a single excitation for given physical resources n and V . Fig. 5d shows the maximum value of p 1 obtained for the case n = 6 as a function of V , with both models predicting that p 1 approaches unity for large V .
Using the effective model, we may analytically obtain an estimate of the error in introducing a single excitation, as a result of population in the second-excited manifold. As we explain in Appendix B 4, in the limit Ω ≪ √ nV eff , where the ensemble may be neglected, the dominant error comes from the finite anharmonicity of the maximum resonance, which yields a small population in state |2 of p |2 ∼ Ω 2 n /4V 2 , where Ω n = √ nκ n Ω. The total error is then wellapproximated by a sum of the error in creating perfect population inversion for an ideal two-level system (under the conditions Ω n , ∆ max , Γ 1D and Γ ′ ), plus the contribution from p |2 . While there is no simple closed-form expression, the total error can nonetheless be easily optimized numerically (see Appendix B 4).

V. SUMMARY AND OUTLOOK
In summary, we have shown that an ensemble of two-level atoms with long-range interactions has rich optical properties in both the single-and many-photon regimes. We have considered the simplest conceptual model, primarily in order to understand and inform experiments that are within reach. Moving forward, it may be possible to exploit strong nonlinearities to produce a blockade effect, where within some region of the system only a single excitation is supported. If the full system is large enough to support several such regions, one would expect the scattered light to exhibit rich spatiotemporal correlations [19].
A. LINEAR OPTICS Figure 2a in the main text shows that for finite values of the characteristic interaction length L, the average energy of the highest single-excitation resonance ω max scales linearly with the number of atoms n in the system. This motivates the definition of an effective interaction strength V eff (L), through the equation ω max (L) = V + (n − 1)V eff (L). The first factor of V is the energy of a single atom, and its average interaction energy with each subsequent atom is V eff (L).

Effective interaction strength
Here we give a simple model for estimating V eff for a photonic crystal of given V , L, and a lattice of N sites separated by distance d.

The band-gap interaction Hamiltonian is
where, for simplicity, we neglect the phase factor resulting from the product of Bloch wavefunctions. When the scaling is linear, it is sufficient to consider the average interaction energy of two atoms, in which case the eigenstates of H BG in the single-photon manifold are simply the symmetric and anti-symmetric states |± = (1/ √ 2)(|g ± |e ), and we may easily determine the corresponding energies. Fixing the first atom to be in the middle of the lattice (i.e. in its average position), the energy of the symmetric state |+ (the state of higher energy), averaged over all possible separations, is where we assume L ≫ d. As shown in Fig. 2b of the main text, this simple model agrees well with the results from exact numerics, however we now consider the limits of its applicability.
The derivation of Eq. A2 makes use of the fact that we know the maximum-energy eigenstate to be |+ , in which the excitation is equally distributed between the two atoms. The assumption that we may use the same expression for the energy in the case of n atoms therefore assumes in turn that the n-atom eigenstate is of the form |ψ n = (1/ √ n) j |e j , which is strictly true only for the case L/d = ∞. The approximation using V eff over-estimates the slope of the linear relationship between ω max and n, as is shown in Fig. A1(a) for several different cases. The accuracy of the approximation therefore diminishes as n increases.
To gain further insight into the validity of the model, we consider the relative error ǫ(L) in the slope V eff , defined as where S(L) is the true slope of the maximum resonance frequency vs. n line, as computed using the full model. Fig.  A1(b) shows this error as a function of L, for the cases N = 100 and N = 200. As one would expect on the basis of the above discussion, ǫ(L) tends to zero as L → ∞. In both cases, the maximum error occurs for L/d ≈ N/2.

Breakdown of linear scaling of ωmax with n
For interaction lengths L/d ≫ 1 the scaling of ω max with n is linear, and moreover in the limit n ≪ N , ω max is well characterized by V eff . In contrast, for very short range interactions (L/d ∼ 1) the coupling between distant atoms can become negligible, resulting in a total energy that scales sub-linearly in the number of atoms n. As a simple example, consider the case of n = 4 atoms in a lattice of N = 200 sites, with an interaction length L/d ∼ O(1). One possible configuration of atomic positions in this system is z j = (1,2,199,200), where the atoms form two disconnected 'blocks'. In this instance, the maximum resonance energy is proportional to that of two interacting atoms, rather than four. Figure A2 shows how the mean maximum resonance energy varies with n for the case L/d = 1, for different total system sizes. The energy is largest (smallest) in the case N = 60 (N = 1000) since the average separation between atoms is shortest (longest). In all cases, the rate of increase in energy diminishes as the system filling increases, in line with the fact that atoms that are increasingly distant couple increasingly weakly.  Fig. A1(a).

Overlap of highest resonance with initial excitation
The TM-mode probe field creates an initial excitation in the system with a phase profile determined by the positions of the atoms. We write this initial state as |ψ in = 1 √ n j e ikLzj |e j , where |e j denotes an excitation at the j th atom, with all others in the ground state. For the choice k L d = π/2 this gives |ψ in = 1 √ n j (i) θj |e j , with θ j = z j /d. Meanwhile, in the limit L → ∞, it may easily be shown that the maximum-energy resonance of the BG interaction has the form |φ max = 1 √ n j (−1) θj |e j , i.e. where each dipole oscillates π out of phase with its nearest neighbour. The overlap κ of the two states is then Since the atomic positions are non-deterministic for any given realization, the latter equality shows that the overlap is proportional to a sum of random phases along the real and imaginary directions. By analogy to random walks, statistically -i.e in taking many samples and averaging -the magnitude of this sum scales as ∼ n/ √ 2, from which it follows that the mean overlap κ scales as 1/ √ 2n. Only in the special case k L d ≈ π does the summation of phases become coherent, where the n dependence cancels and κ → 1.
The above reasoning breaks down as the number of atoms approaches the total number of lattice sites, since then the sum contains an equal number of terms of opposite phases and thus tends to zero. For example, for n = N = 200, there are 50 terms of each of the phases 1, i, −1, −i, which when summed give a vanishing overlap. As shown in Figure   2c in the main text, we find that the simple formula √ κ n = (N − n)/ √ 2nN agrees very well with the numerically obtained scaling for all n.
The above reasoning also assumes that L → ∞, however, as Fig. A3 shows for a lattice of N = 200 sites, the same scaling behavior for κ holds for interactions with L/d > ∼ 1. Below this limit, the BG interaction is no longer dominant and thus no longer determines the maximum-energy resonance.

Effective model for transmittance at maximum resonance
Here we bring together the key features of the maximum-energy resonance that we have found in the linear regime, constructing an effective model that aims to reproduce the results of the full model of Eq. 1 in the main text. Specifically, we consider the transmittance T dip as illustrated in Fig. 1d of the main text. We define the detuning of the probe field from the maximum-energy resonance to be ∆ max = ω L − ω max . The ingredients for the model, assuming L/d to be sufficiently large, are: • An ensemble of n independent atoms at the frequency −(∆ max + nV eff (L)), denoted by a single level |E , whose coupling to the input field is enhanced by a factor of √ n, and whose decay into the guided modes is enhanced by a factor of n. This approximately captures the response of a normal atomic ensemble driven far from resonance.
• A single degree of freedom |1 at frequency ∆ max = 0 representing the maximum-energy resonance. In order to determine the coupling of the input field to this state, we note that the matrix element for coupling of the field to the initial atomic excitation |ψ in is ∼ √ nΩ, while the subsequent overlap with the maximum-energy resonance is given by ∼ √ κ n , as described in the previous section. Overall, therefore, the matrix element coupling the input field to the maximum resonance is √ nκ n Ω.
• For L/d → ∞ the frequency of the maximum resonance is independent of the atomic positions, whereas for finite L/d each configuration has a different value of ω max . As a result, in the averaged spectrum there is less attenuation at a given frequency in the latter case than the former, and hence the transmittance is higher. For the case N = 200, L/d = 50 presented in the main text, numerical analysis shows the spread of ω max to be well approximated by a Gaussian distribution, so that for given values of L and n we can obtain the standard deviation σ(n, L), which we then input this to the model via a random number η with the appropriate statistics. Finally, we average over many samples of different η. Figure 3a in the main text schematically illustrates the level diagram for the model, with the various frequencies, couplings, and decay rates. The Hamiltonian is Here, as described above, η is a random number of zero mean and standard deviation σ(n, L). We emphasize that we do not rigorously derive this effective Hamiltonian from the full Hamiltonian of Eq. 1 in the main text, but rather it represents a minimal model meant to capture the most important characteristics of the system. The output operator for the transmitted field is where the noise operator F describes quantum jump processes, which may be neglected for weak driving (Ω/Γ ′ ≪ 1). Using this model we obtain the results of Fig. 2d (crosses) in the main text, taking 1000 samples of different η in the case L/d = 50.
In the case L → ∞, we can also use the effective model to analytically obtain a simple formula for the transmittance T dip . Assuming √ nV eff ≫ Ω, nΓ 1D , the ensemble |E may be neglected and the problem is reduced to that of a two-level system comprising the levels |g , |1 . Solving the corresponding optical Bloch equations in the steady state (SS), we find that the transmittance is given by (A7)

Spectra for random number of atoms
In addition to averaging over atomic configurations for a fixed number of atoms, we can also take into account the effects of fluctuations in the number of atoms in the system. For a given average number of atoms n we use a Poisson distribution to give a random number m of atoms, generate a random configuration of positions of these m atoms, then obtain the transmission spectrum; we then repeat many times and take the average. As figure A4 shows, the implication for the maximum resonance is that it can now occur over an even larger range of energies, so that in the averaged transmission spectra its signature is flattened out.  Figure 4 in the main text shows that there is a very high probability of observing strong photon anti-bunching in the reflected-field correlation function g (2) R (0) for V > Γ ′ when the system is driven at its highest single-photon resonance frequency. It is natural to then consider whether the same is true when the driving is instead at the frequencies of other resonances in the single-photon manifold. To investigate this, we diagonalize the interaction Hamiltonian H BG in the single-excitation manifold to obtain a set of (in general non-degenerate) resonances, which we index by m = 1, ..., n from highest to lowest energy (i.e. m = 1 is the maximum-energy eigenstate). We then weakly drive the system at each frequency ω m and compute the correlation function g

Anharmonicity for finite L
In the limit L → ∞ we argued in the main text that the anharmonicity -defined as A = ω max ) is the highest eigenvalue of the second (first) excited manifold -is given by A(∞) = −2V . For finiterange interactions and n ≪ N , as discussed in Appendices A 1 and A 2 above, ω (1) max is characterized by an effective interaction energy V eff (L). It is then natural to ask whether in the finite-range interaction regime, the anharmonicity can also be characterized by V eff (L): here we investigate this question numerically. Figure B2 shows the mean anharmonicity A (solid dots, with averaging over 1000 samples of atomic positions), normalized to the interaction strength V , as a function of L for two different system sizes, N = 50 and N = 200, with n = 20 atoms. In addition, we plot (dashed lines) the hypothesized anharmonicity −2V eff (L). In both cases we see that for large L/d > ∼ N/2, the approximation based on V eff (L) is very accurate, while for small L there is a significant difference between the two models, which clearly increases for larger systems.

Effective model for many-photon dynamics
Here we describe the effective model for the multi-photon behavior, which may be used to make approximate predictions when either the number of atoms or excitations becomes sufficiently large that solutions of the full master equation become infeasible. We build on the linear model described in Section A 4, assuming now that the system may be described in terms of two bosonic-type modes: (1) a mode corresponding to the 'ensemble' of n independent atoms, which can support up to n excitations, and; (2) a spin-wave corresponding to the maximum resonance |1 of the first-excited manifold, and its doubly-excited state |2 in the two-photon manifold. The Hamiltonian is given by Here, σ µν = |µ ν| for energy eigenstates |µ and |ν , |E m denotes the m th excited state of the ensemble, with corresponding energy ω m = −m(∆ max + nV eff ), and |E 0 = |g , and all other quantities are as defined in the effective model of Appendix A 4. Dissipation is described by the Lindblad operator where Γ n = Γ ′ + nΓ 1D is the decay rate of the ensemble, and γ n = Γ ′ + nκ n Γ 1D is the decay rate of excitations in the spin wave. We now illustrate the predictive power of the effective model, numerically solving the master equatioṅ ρ = −i[H eff , ρ] + L eff [ρ] in two example cases, and comparing with the results of the full model of Eq. 1 in the main text. Furthermore, in Appendix B 4 we use the effective model to obtain simple analytical results in the regime where the system is driven close to the resonance frequency ω max . Fig. B3(a) shows the predictions of both the effective and full models for the maximum first-excited state population as a function of the number of atoms n, optimized over the detuning ∆ max , for a driving strength Ω/Γ ′ = 1. Here, in order to obtain the solution for the full model, we first verify that for small n < ∼ 6, the dynamics for this particular choice of Ω are very well approximated by truncating the Hilbert space at the second excitation manifold (i.e. by discarding states with 3 or more excitations). Since for larger values of n the maximum-energy resonance is driven more weakly (due to the reduced overlap with the input field), and hence the probability of introducing additional excitations is suppressed, we then obtain solutions for up to n = 20 with the truncated Hilbert space.
Rather than averaging the dynamics over many configurations, in this case we simply search through 1 million samples of random atomic positions and, for each n, choose the configuration whose maximum single-photon resonance |φ max has an overlap with the input excitation |ψ in that is closest to the overlap √ κ n used in the effective model. In other words, for each n we seek the configuration of atoms that minimizes the quantity D = | √ κ n − φ max |ψ in |.
We see from the figure that the two models agree reasonably well in their predictions for the maximum fidelity of introducing only a single excitation, with an average error (i.e. the average absolute difference between the values predicted by the models, normalized to the value from the full model) of approximately 7.5%. Since in general it is not possible to find a configuration of n atoms such that the difference in overlaps D is identically zero, we expect there to be some degree of discrepancy between the two cases. Figure B3(b) shows D for each of the cases in Fig.  B3(a), illustrating that the agreement between the two models is best when D is smallest. Having verified the validity of our effective model, we now apply it to the situation described in the main text. In particular, we apply a global Rabi pulse to all atoms at a frequency near the maximum resonance, and consider the fidelity of producing only a single excitation. Figure B4 shows the prediction for the optimal single-excitation probability p 1 as a function of interaction strength V and number of atoms n, optimized over both driving strength Ω and detuning ∆ max . We see that as the interaction strength increases, the fidelity of introducing a single excitation is very high regardless of the number of atoms, since by increasing the strength of the driving one can overcome the weaker coupling of the input field to the spin wave.

Analysis of single-excitation probability for driving close to ωmax
In this section, we use the effective model to analytically quantify the errors that reduce the fidelity of introducing only a single excitation. We assume that the dominant errors are (1) dissipation via emission into free space and into the waveguide, and (2) population of the doubly-excited manifold. We begin by considering the latter type, and will show that under the appropriate conditions these errors become negligible, in which case the fidelity is then limited only by errors of type (1).
The system can be doubly-excited by acquiring population either in the state |2 , or in the second excited state |E 2 of the ensemble. In order to illustrate how the populations of these states scale with the system parameters, we assume the system to be lossless and solve the Schrödinger equation (SE) for the atomic wavefunction: the generalization to the case with dissipation is straightforward, and does not qualitatively affect the following conclusions. The SE reads i˙ c = H c, where the state vector c = (c g , c E1 , c 1 , c E2 , c 2 ) T , c i are the probability amplitudes of the energy eigenstates (labelled by i) of the model, and the Hamiltonian is given by Eq. B1.
We are interested in the regime where there is a high probability of exciting state |1 , so by taking Ω ≪ √ nV eff we can strongly suppress excitation of the ensemble. Solving for the amplitudes c E1 and c E2 , we find (taking -for simplicity -L/d → ∞, such that V eff → V )