Nondestructive photon counting in waveguide QED

Number-resolving single-photon detectors represent a key technology for a host of quantum optics protocols, but despite signiﬁcant efforts, state-of-the-art devices are limited to few photons. In contrast, state-dependent atom counting in arrays can be done with extremely high ﬁdelity up to hundreds of atoms. We show that in waveguide QED, the problem of photon counting can be reduced to atom counting, by entangling the photonic state with an atomic array in the collective number basis. This is possible as the incoming photons couple to collective atomic states and can be achieved by engineering a second decay channel of an excited atom to a metastable state. Our scheme is robust to disorder and ﬁnite Purcell factors, and its ﬁdelity increases with the atom number. Analyzing the state of the re-emitted photons, we further show that if the initial atomic state is a symmetric Dicke state, dissipation engineering can be used to implement a nondestructive photon-number measurement, in which the incident state is scattered into the waveguide unchanged. Our results generalize to related platforms, including superconducting qubits.


I. INTRODUCTION
Single-photon detectors have a long history, with a plethora of technologies available [1].Applications in quantum optics, such as quantum state preparation, quantum metrology [2], entanglement distribution [3], and quantum computing [4,5] have placed a renewed focus on single-photon detectors capable of resolving the number of incoming photons.Perhaps most promising are superconducting transition-edge sensors, which have been demonstrated to achieve a (perphoton) detection efficiency (percentage of detected photons) of η ≈ 95% [6] and to distinguish up to seven photons, with a negligible dark count rate (clicks in the absence of incoming photons).They are based on the principle that near the critical temperature, the resistance of a superconductor is very sensitive to temperature changes, down to the level of single-photon energies.While very impressive, the device is limited to optical photons, destroys the photonic state, and is difficult to scale.
One strategy that in principle also allows for nondestructive measurements is based on quantum memory.An itinerant photon may be caught by an atom in a cavity [7] or by a cavity with tunable coupling [8], which however requires time-dependent tuning of the atom-cavity or cavity-waveguide coupling in accordance with the photon wave packet shape.More generally, electromagnetically induced transparency (EIT) can be used to slow down a light pulse such that it fits within an atomic cloud [9,10].This allows the storage of arbitrary photon pulses within a length given by system parameters.The number of polaritons can in principle be read out by employing a cycling transition [11,12].However, this scheme also requires knowledge about the arrival of the photon wave packet.Furthermore, the combination of quantum memory and quantum nondemolition detection of one or few excitations in a 3D cloud of atoms makes this very challenging.EIT in a waveguide QED setting is one strategy to alleviate the imaging problem (see, e.g., Ref. [13]).
Here, we instead design a detector that does not require knowledge about the shape or arrival time of incoming wave packets and therefore does not require implementation of a time-dependent Hamiltonian.This can be achieved either through continuous measurement, for example through dispersive coupling [14][15][16][17][18][19], or if the detector permanently changes its state and is read out later, as in impedancematched -systems [20][21][22][23][24][25][26].The former class of detectors suffer from measurement backaction on the photonic state [27,28], whereas the latter does not produce "clicks" to indicate detection event, and therefore does not localize photons in the waveguide, which avoids measurement backaction.The advantage of this strategy is that no knowledge about the photon wave packet is needed, and that photons do not have to be destroyed to be detected, which opens the way toward nondestructive photon-number measurements.
In recent years, quantum emitters coupled to waveguides have emerged as a powerful experimental paradigm [29,30].Examples of such systems include cold atoms levitated near optical fibres [31,32] or photonic crystal waveguides [33], but also solid-state realizations such as quantum dots [34,35], superconducting qubits [36,37], or nitrogen-vacancy centres [38][39][40].Due to the strong confinement of light, even a single emitter can have a profound effect on light propagation, and many emitters show remarkable collective effects [41].At the same time, impressive experimental breakthroughs have FIG. 1. Sketch of the proposed setups.We consider arrays of quantum emitters coupled to (a) a semi-infinite waveguide terminated by a mirror at x = 0, (b) an infinite waveguide, and (d) to two waveguides.The decay e → g couples to the waveguide (blue, g ), but also to free-space modes (red, free ).Key ingredient is a tunable engineered decay from the excited |e to a metastable state |s [green, s ( )], implemented for example as a Raman transition as shown in (c).The photons may be emitted either into free space [(a) and (b), Sec.III)] or back into the waveguide for a nondestructive measurement [(d), Sec.IV)].enabled the control [42] and, in particular, readout of neutral atoms, 1 superconducting qubits [43,44], and trapped ions [45].Thus quantum emitters coupled to waveguides appear to be an ideal platform to produce [33,[46][47][48] and control [49,50] quantum light.Here we show that they also allow for detection of quantum light.
The key idea here is to engineer atoms such that for each incident photon in the waveguide, exactly one atom changes its internal state, such that a subsequent measurement of the atomic state yields the number of photons in the scattered wave packet.We do this by identifying conditions such that all photons are absorbed in one atomic transition (g → e, blue in Fig. 1) and dissipated in a different one (e → s, green).This way, the atomic array keeps a memory of the number of scattered photons [20][21][22][23][24][25].
In order to identify such conditions, we study a realistic model for photon absorption by atomic arrays in the presence of free-space decay (red) and disorder.In view of experimental realizations, we provide recipes for engineering an additional decay channel, and pay particular attention to spatial disorder, as this fundamentally modifies the eigenmodes of the system.Surprisingly, we find that it is possible to engineer dissipation to achieve full absorption independent of disorder.As a consequence of collective enhancement in the atom-waveguide coupling, detection efficiency and bandwidth grow with atom number, such that even with moderate Purcell factor (ratio of waveguide to free-space decay rate), detection 1 In an array of 160 atoms, a recent experiment achieved 99.94% readout fidelity, albeit at the expense of a constant atom loss rate that dominates this theoretical fidelity [86].Atoms trapped in tweezer arrays can circumvent this problem, offering similar fidelity without losing the atoms [87].
efficiencies η > 99% can be reached for intermediate numbers of atoms (N > 20).We also consider emitters coupled to chiral waveguides, which have distinct advantages due to the natural suppression of backscattering.Our detection scheme is scalable, as the number of atoms required to detect a certain number of photons with a fixed error grows only polynomially.
In a second part, we study a natural modification, in which the dissipated photons are emitted back into the waveguide [cf.Fig. 1(d)].If the resulting output wave packet coincides with the input wave packet, this realizes a quantum nondemolition (QND) measurement. 2In order to achieve high fidelities, this requires collective enhancement of both decay channels, which requires the atomic array to be prepared in a symmetric Dicke state between ground states.We show that such states can be prepared through coherent interaction of two arrays and subsequent conditioning.Ultimately, this can reach a scaling where for a given error probability, the number of atoms that can be prepared in this way scales exponentially with Purcell factor N ∼ exp(P), which we verify with numerical simulations.
The rest of this paper is structured as follows.In Sec.II, we present the scattering theory for atomic arrays on waveguides.Building on the principle of coherent perfect absorption, we show explicitly that perfect absorption can be obtained by tuning only dissipation, even in the presence of arbitrary disorder.We apply this principle to establish the feasibility of number-resolving detection by absorption in a mirror geometry (Sec.III A) and in an infinite waveguide (Sec.III B).For completeness, we also discuss how dissipation can be used in a chiral waveguide to obtain a number-resolving detector (Sec.III C).We extend these concepts to a QND measurement in Sec.IV.In Sec.V, we discuss the experimental implementations of our proposal, including engineered dissipation, atomic species, readout, detector reset, the effect of nonidealities, and bandwidth.We conclude in Sec.VI.

II. ABSORPTION IN ATOMIC ARRAYS
Below, we first formulate a generic input-output theory for photons scattering off atomic arrays (Sec.II A), following previous work [13,51,52].This allows us to generically identify the conditions under which all photons are absorbed in one transition and emitted in the other (Sec.II B).

A. Scattering theory
A generic Hamiltonian describing the system-waveguide interaction is given through [13,51,52] In Eq. ( 1), the label ν runs over different sets of waveguide modes, the index i denotes whether the waveguide mode couples to the g ↔ e transition (i = 1) or the s ↔ e transition (i = 2), the index n runs of the N atoms, the individual waveguide modes are labeled by their wave vector k, and the coupling of them to the atoms is given by the rate g (n) ν,k,i .This Hamiltonian describes a situation where the two transitions are coupled to different fields, which can either be waveguide modes or free-space modes.In an infinite waveguide [Fig.1(b)], g (n)  ν,k,1 = 2c g exp(ikνx n ) and ν ∈ ±, corresponding to left-and right-moving modes, whereas for the semi-infinite waveguide [Fig.1(a)], there is only one set of waveguide modes with coupling g (n)  k,1 = c g sin(kx n ).In these expressions, g is the decay rate to |g of an individual atom into the waveguide, and ω k = ck.Integrating out the bath modes yields quantum Langevin equations for the spin operators (see Appendix A), where the sum over input fields ν and atoms m is implied.The non-Hermitian Hamiltonians H eff,1 (H eff,2 ) describes both coherent interaction of the quantum emitters and decay into the waveguide induced by the coupling of the g ↔ e (s ↔ e) transition to the waveguide, and a in,ν,i are waveguide fields coupling to these transitions.This coupling is given by the generically nonsquare matrix L i .If coupled via an infinite waveguide, the ith transition couples to two input fields, a in,i = (a in,+,i , a in,−,i ), corresponding to right-and leftmoving photons, such that L i is a N × 2 matrix.This is in contrast to a semi-infinite waveguide, which only has on input field, such that L i is a vector.Independent free-space decay at rate free,i , can also be captured by adding a term H free,i = −i free,i 1/2 L 3 = √ free 1, and introducing N independent noise operators a in,free,i .Note that the assumption of independent decay into free-space does not necessarily hold and depends on how closely spaced the quantum emitters along the waveguide are.Such a situation may lead to a reduction in free-space decay [53], which would benefit the detector proposed here.Multiple baths coupling to the same transition appear as several terms taking either the form of the first or the second line on the right-hand side.
The Langevin equation for σ (n)  se can be obtained by exchanging both 1 ↔ 2 and g ↔ s everywhere in Eq. ( 2).The explicit form of H eff for the mirror geometry and the infinite waveguide are given in Eqs. ( 7) and (10) below.
If there are only few excitations compared to the number of atoms, one can linearize the Langevin equation using a Holstein-Primakoff transformation that sends σ ge → b [13,29,54,55].In this approximation, σ (n)  gs σ (m) se → δ mn σ (n) ge , such that only the diagonal terms in the second line of Eq. ( 2) survive H eff,2 → −i s 1/2 (green in Fig. 1), independent of whether the decay e → s corresponds to guided or nonguided modes.Combining this with free-space decay (shown in red in Fig. 1), we obtain uniform incoherent decay at rate = free,1 + s .We thus arrive at Here, we have neglected all input fields except the ones pertaining to the waveguide, which is valid if they are in vacuum.In order to describe destructive photon measurements, it is then sufficient to show that all incoming photons are absorbed and re-emitted into the bath modes coupling to the s ↔ e transition.Strictly speaking, decay to |s eliminates the atom from the dynamics, but this effect is neglected in the linearization.

B. Complete absorption
We now examine in general how tuning in Eq. ( 3) can lead to complete absorption of photons by an atomic array.For the rest of this section, we drop the indices from H eff and L, for sake of generality and simplicity.In order to find the linear scattering properties of the array, we use the input-output equations that relate the output field operators to the input fields a out (t ) = a in (t ) − L † b(t ) [54].Solving Eq. ( 3), we obtain the scattering matrix A detector that counts the number of photons in a specific input port (say, a in,+ ) needs to absorb all of them and dissipate them via the transition to |s .This is captured by our key figure of merit, the detection efficiency η, which is the product of the probability that a photon is not reflected 2 , and the probability that it is dissipated via the engineered channel (rate s ) rather than into free space, η = p abs s /( s + free ).Thus a high fidelity requires s free and p abs ≈ 1.In waveguide QED, this regime can be reached through the collective enhancement of the emitter-waveguide coupling as compared to the free-space decay.
Unity absorption is attained if one of the eigenvalues of the scattering matrix S(ω) is zero.This corresponds to a pole of the inverse scattering matrix A pole of S −1 arises whenever ω coincides with an eigenvalue of H † eff − i /2, which implies that the scattering matrix S has a zero if ω − i /2 coincides with an eigenvalue of H eff .Thus tuning ω and , this can always be achieved, for any H eff .Absorption based on this principle has been observed in a variety of systems [56][57][58], and been termed coherent perfect absorption [59,60].Note that to reach this conclusion we did not have to assume anything about the form of H eff (apart from linearity), which is the reason it works for arbitrary disordered systems.If these conditions are fulfilled, there exists an eigenvector e 0 , such that S(ω 0 )e 0 = 0. Generically, e 0 describes a linear combination of various input fields of the system, which implies that coherent perfect absorption is only a sufficient condition for unity absorption property when there is only one input field, such as in the mirror geometry.Nevertheless, in the infinite waveguide efficient detection is still attained for large atom numbers.
One can verify explicitly that S −1 is the inverse scattering matrix (here, Γ is a matrix, for generality) FIG. 2. Photon detection in mirror geometry.(a) Photon loss probability (p loss = 1 − η) on resonance as a function of atom number N for an array in the atomic mirror configuration (lattice spacing a = λ) in the presence of weak (black), intermediate (blue) and strong (green) spatial disorder (averaged over a normal distribution with standard deviation σ ≡ δx 2 = {0.01,0.08, 0.2}λ, where λ is the wavelength of light).The standard deviation of each curve is shown as lightly shaded area, averaged over 2000 realizations.(b) Same system as (a) , but with engineered dissipation set to the average expected largest dissipation rate s = Im[ μ σ,N ], which mitigates the saturation in detection efficiency.In red we show the result for a completely random array, but with fixed (characterizable) disorder.(c) Photon loss as a function of frequency calculated for a completely disordered array of N = 200 atoms when tuned to the first, second, third, and fourth eigenvalue (from dark to light blue).Purcell factor for all plots is P = g / free = 10. where The square brackets in Eq. ( 5) vanish, which can be checked explicitly for the two examples below.It can also be shown to hold generically, since if Γ = 0, the scattering matrix is unitary S −1 = S † , which is only true if the term in square brackets vanishes.
Physically, this can be interpreted as the fluctuationdissipation theorem, as the anti-Hermitian part of H eff specifies the damping, whereas L captures how strongly the modes are coupled to the input noise operator.

A. Mirror geometry
We now turn to the specific model that best illustrates these ideas, namely an array of atoms coupled to a waveguide terminated by a mirror as sketched in Fig. 1(a).In this geometry, the non-Hermitian Hamiltonian induced by integrating out the waveguide photons reads where g is single-atom decay rate for the transition e → g into the waveguide, x n the position of the nth atom, and k 0 is the wave vector of the emitted light (wavelength λ = 2π/k 0 ).The coupling of the atoms via the waveguide contains both a term due to photons traveling directly in between them, accumulating a phase k 0 |x m − x n |, and one mediated by photons being reflected from the mirror, which incurs a minus sign and a phase k 0 (x m + x n ).Since there is only one input and output field [cf.Fig. 1(a)], the matrix L is now a vector L n = g sin(k 0 x n ).It is instructive to see an example of how perfect absorption manifests in this setup.Placing the atoms in the atomic mirror configuration at positions x n = (1/4 + n)λ (due to the infinite-range interactions, the lattice need not have unity filling), the photonic field only couples to the symmetric collective atomic excitation B = n b n / √ N, which also is an eigenmode of the atomic array.All other modes are dark and do not participate in the dynamics.In terms of this collective mode, the governing equations reduce to the input-output equations for a one-sided cavity [61] with internal dissipation where we have introduced the total decay rate tot = N g + .As in our discussion above, = free + s comprises both free-space decay and the decay from e → s, which is assumed to be tunable.Solving Eqs.(8a) and (8b) in frequency space, we calculate the number of photons in the output field If decay is tuned such that tot = N g + = 2N g , there is perfect absorption on resonance (p abs = 1), with a bandwidth of 2N g , which in this single-mode picture is equivalent to impedance matching or critical coupling.In order to match the collective decay, the engineered decay s ∝ N, such that as the atom number is increased, the detection efficiency η = p abs s /( s + free ) can become arbitrarily close to 1.In Fig. 2(a), we include spatial disorder and show how the photon loss on resonance p loss ≡ 1 − η scales with atom number.Clearly, while the setup works very well for low spatial disorder (σ /λ < 1%), it suffers significantly from disorder.In the following we show how this is mitigated.Note that in Fig. 2 and indeed all plots in this paper we choose the Purcell factor P = 10, which close to the current state of the art [62].However, the collective enhancement of the atom-waveguide coupling means that the primary effect of having a lower Purcell factor is that more atoms have to be employed and therefore read out.
In the presence of disorder, the energies and decay rates of the eigenmodes of the atomic array are shifted, and many collective atomic modes couple to the input field.Thus the picture presented above breaks down.Yet, as we have shown following Eq.( 4), full absorption can be attained generically, independent of disorder, by tuning the engineered dissipation s to one of the eigenmodes.However, the eigenmodes in the presence of disorder are not known a priori, so this approach appears infeasible.Surprisingly, one can still vastly improve over naively setting s = N g as we did in Fig. 2(a), if the standard deviation σ ≡ δx 2 of atomic positions is known.For a given N, σ , one can then calculate the average largest eigenvalue μ N,σ and tune the engineered dissipation to its imaginary part s = Im[ μ σ,N ].Note that Re μ σ,N = 0 if every configuration is as likely as its reflection.This restores the favourable scaling of detection efficiency with N, as illustrated by the blue (σ = 0.1λ) and green (σ = λ) curves in Fig. 2(b).Most strikingly, this works even in the presence of disorder equal to the lattice spacing (green), which is essentially equivalent to a fully random configuration.The reason it works lies in the fact that the largest eigenvalue and thus the absorption bandwidth grows faster than the fluctuations of the largest eigenvalue around its mean.
If the disorder is fixed as a result of fabrication, such as in solid-state implementations, one can further improve the scaling by measuring the largest decay rate.In this case, s can be tuned exactly to the largest eigenvalue.This situation corresponds to the red curve in Fig. 2(b), calculated for completely random configurations, where in each case s was set to coincide with the imaginary part of the largest eigenvalue.This clearly leads to a better result than without exact tuning (green).The scaling is still worse than with low disorder (blue), because the largest eigenvalue grows slower in completely disordered configurations compared with mostly ordered ones.
So far, we have just discussed absorption and detection on resonance.Equally important is the detection bandwidth, given by the engineered decay rate.Since also the detection efficiency η depends on the ratio between engineered dissipation and total decay rate, best detection is achieved when tuning to the most dissipative eigenmode [cf.Fig. 2(c)].This is our choice for all plots.

B. Infinite waveguide
Let us now turn to an atomic array coupled to an infinite waveguide, which has the simpler effective Hamiltonian since there is only one path for a photon to travel from one atom to the next.As illustrated in Fig. 1(b), there are now two input and two output modes, a right-moving one (+) and a left-moving one (−) [cf.Eq. ( 3)].The atomic lowering operators couple to the input operators via the N × 2 matrix L n,ν = g exp(ik 0 νx n ), where ν ∈ {±1} labels right-and left-propagating modes.The scattering matrix reads where M is defined as above [Eq.( 6)].Since M is symmetric, transmission of right-and left-moving waves (the diagonal elements of S) are equal.
One can show that for atoms arranged in a periodic array, parity symmetry implies that reflection amplitude of rightand left-moving photons differs only by a phase.This is because JMJ = M under the action of the exchange matrix J, which consists of ones on the antidiagonal and otherwise zeros, and JL + = exp(iφ)L − for some φ.In this case, we can parametrize the scattering matrix as with eigenvalues μ = A ± B exp(iφ).In this system, coherent perfect absorption (one eigenvalue is zero) is equivalent to B = exp(iθ )A and exp(iφ + iθ ) = ±1.In the end, full absorption of a uni-directional wave packet may only be attained if S = 0. Thus parity symmetry implies that perfect absorption may only occur if the scattering matrix is zero, corresponding to an exceptional point.Interestingly, as the atom number N → ∞, the scattering matrix S → 0 for all arrays except the atomic-mirror configuration.In the latter, the scattering matrix can be shown to reduce to (for full absorption s = N g ) which gives perfect absorption for wave packets that are symmetric superpositions of left-and right-propagating modes, but not for wave packets incident from one direction, an effect seen before [21,63].This can be overcome through atomic lenses [64], a mirror as above [23], or by using any other lattice spacing [21].The latter can be deduced by estimating the magnitude of the elements of the scattering matrix through |t| N ∼ exp(−N 1−α ), having assumed that the largest eigenvalue scales as N α g .Only in the atomic mirror configuration is α = 1, otherwise α < 1.
We numerically check the behavior for a range of other spacings as well as fully disordered arrays and find similar behavior as long as the atomic positions differ sufficiently from the atomic mirror configuration.In Figs.3(b) and 3(c), we demonstrate this for a selection of array spacings (without including disorder).We choose k 0 = π/(2λ) for no particular reason, except that this appears to be a good choice.Including disorder, and employing the same technique of tuning to the average largest eigenvalue, we find similar scaling behavior as in the mirror geometry, illustrated in Fig. 3(a).The upshot is that arbitrary detection efficiencies can again be attained by increasing atom number, independent of disorder.

C. Chiral atom-waveguide coupling
Interestingly, the recently demonstrated platforms for chiral atom-waveguide coupling [65] are another architecture in which robust photon detection may be achieved.Such a coupling is realized in a range of situation, for example when the light field is strongly confined [66][67][68][69], when giant atoms are tuned to give a chiral coupling [70], or in topological systems [71,72].
By design, (almost) no backscattering occurs in these systems, there are no collective effects, and the spacing of the atoms is immaterial, making the analysis straightforward.If an atom coupled to right-and left-moving modes at different rates, transmission and reflection on resonance are captured This graph illustrates the fact that an array in the atomic-mirror configuration (a = λ) cannot serve as photon counter, whereas all other generic spacings work similarly well, a fact that can be understood analytically (cf.Sec.III B).(b) Probability for an undetected photon on resonance for an atomic array coupled to an infinite waveguide as a function of atom number.The blue line denotes the limit of a perfectly ordered array with spacing a = λ/4 (or 5λ/4 etc).Purcell factor is P = 10, the average for the other cases was performed over 2500 disorder realizations (standard deviation shown as lightly coloured area).(c) Scaling of largest eigenvalue in ordered arrays with varying spacing compared to a fully random one (red).While the atomic mirror configuration (a = λ) is clearly different, it is a fine-tuned exception, with all other generic arrays (ordered or disordered) behaving remarkably similar.The robustness of our scheme relies to a large degree on this universal eigenvalue scaling.
by [65] As before, = s + free .This allows us to calculate the peratom absorption probability (for the + mode) The probability that the photon is dissipated in the right channel is s / as before.In the limit of many atoms, the transmission probability vanishes, as all photons are either reflected or absorbed.The detection efficiency is therefore the ratio of absorbed photons times the probability they are dissipated to |s , which to first order in γ − /γ + is given by Note that even for moderate γ − /γ + , the second term can be reduced arbitrarily by increasing s , with the caveat that a larger number of emitters is needed before complete extinction is attained.In the absence of backscattering, this scheme is intrinsically robust against disorder.On top of that, the detection bandwidth depends only on the bandwidth of chirality and thus is-at least in principle-independent of s .This comes again with the caveat that photons far detuned from resonance on a scale of s can only be absorbed with a large number of emitters.

A. Outline
In what we have discussed so far, we have disregarded the photons emitted via the engineered decay.After the state of the atomic array is measured, the only information obtained concerns the number of photons in the pulse.Since for each photon absorbed, one is emitted via the engineered channel, the question is pertinent whether a situation can arise in which the emitted photonic state coincides with the input state.This requires the outgoing photons to be disentangled from the atoms, save for in the collective number basis, and that the photonic state is not distorted in any other way.If these conditions are fulfilled, the setup realizes a quantum nondemolition (QND) photon-number measurement.
It is immediately obvious that if the photons are emitted into free space and thus scattered in all directions, the outgoing photonic state is (a) useless, (b) still entangled with the atoms, and (c) distributed over many different modes.In principle, this is mitigated to some degree if the photons are emitted back into the waveguide.However, if each atom were to emit independently, the probability of any one photon getting lost is 1/(1 + P), where P is the Purcell factor, which severely limits the fidelity of the scattered wave packet for realistic Purcell factors P. Furthermore, a photon has a different phase depending on the atom it is emitted from, causing residual entanglement of the photons with the atomic state, except in the atomic-mirror configuration.Thus QND detection requires the mirror geometry [cf.Fig. 1(d)], which we study in detail below. 3t turns out that the same trick that can make absorption robust against free-space decay and disorder-collective decay-can also be used to protect the re-emission into the waveguide.This is achieved if the atoms are in a superposition of |g and |s instead of all in |g .As we demonstrate mathematically below, this yields collective enhancement of the atom-waveguide coupling for both decay channels.This mode of operation therefore holds one further big advantage over the non-QND operation, namely that the bandwidth of the overall detector is not limited to the single-atom bandwidth.If the final readout of the atomic state is to give information about the number of photons, the ground state superposition must initially posses a known number of atoms in |g .These requirements are fulfilled by Dicke states, fully symmetric states with a definite number of excitations where sg is a collective spin operator and |G is the state in which all atoms are in |g .

B. QND measurement with Dicke states
We now aim to describe how a number-resolving QND measurement works in practice.We first focus on analytical results and approximations to motivate and illustrate the idea, and then corroborate the results with numerics.In the atomic-mirror configuration (a = λ) the coupling via a semiinfinite waveguide (7) simplifies to H eff,1,mn = −i g /2.We also assume that the e → s transition couples to either the same waveguide at a slightly different frequency or to another waveguide field, such that H eff,2,mn = −i s /2.In this case, the Langevin equations ( 2) reduces to a single equation of motion for the collective spin operator and another one with g ↔ s and 1 ↔ 2, where αβ are collective spin operators.The input-output equations in this case read ) Note that neglecting the input fields corresponding to free-space decay in Eq. ( 18) is an approximation, as they take the system out of subspace of Dicke states.The reason is that when a photon decays into free-space, it destroys the coherence, which leaves the system in a mixed state.Therefore Eq. ( 18) fails to account properly for the dynamics after the first photon has been dissipated to free space.The reason we can still use it to calculate the fidelity of QND measurements of Fock states is that once a photon is lost, the fidelity immediately drops to zero and remains zero, such that the subsequent time evolution of the system is immaterial.
To understand how impedance matching can be attained here, consider the system being in the symmetric state with m e excitations |ψ = |N/2 − m 0 − m e , N/2 + m 0 , m e , where we have defined the imbalance m 0 = m − N/2.Acting on this generic symmetric state, Since this is true independent of m 0 and m e , it is an operator identity, but only in the subspace of symmetric states.Another identity can be obtained by exchanging s ↔ g.This allows us to rewrite the equation of motion in the simplified form where ξ in is the same input noise as in Eq. (18).In symmetric states with sufficiently large population in both ground states, we can replace the operators approximately by their expectation value in the initial state, S gg − S ee ≈ N/2 − m 0 and S ss ≈ N/2 + m 0 .Thus we conclude that ensures impedance matching, such that the decay rate of an excitation via the first and via the second channel are equal up to O(1/N ).Under this condition, the equation of motion (21) on resonance can be solved approximately S ge ≈ a in,1 / g .On the one hand, this implies that a out,1 is independent of a in,1 , since the whole signal is absorbed.On the other hand, a out,2 ≈ −S sg a in,1 2 s / g N 2 .Given this, a † out,2 a out,2 ≈ a † in,1 a in,1 , i.e., every photon incident on port 1 is transmitted to port 2. The presence of S sg in this expression indicates that for each incoming photon, an atom is transferred from g to s.This still holds for the second, third, ..., nth photon, but with an error of order O(n/N ), which is the same as for the non-QND detector, up to constants of O (1).
To make this intuition quantitative, we follow a simple argument.First of all, let us denote the overlap of the output with the input wave packet in case of the scattering of a single photon in a given state as For any finite-bandwidth wave packet, this differs from unity.For example, a single photon with a Gaussian wave function of width τ is transmitted by an impedance-matched cavity with outcoupling rate κ and internal decay rate κ int with fidelity where erf (x) = (2/ √ π ) x 0 exp(−y 2 )dy is the error function.To make contact with our simulations, we thus set κ = g (N/2 − m 0 ) and the internal cavity decay κ int = free .In this limit of few excitations, it is therefore valid to think of free-space decay as limiting the fidelity by introducing a branching ratio between decay back into the waveguide and decay into free space, which is captured explicitly by the prefactor in Eq. (23).
Given some single-photon fidelity F 1 , a linear device would transmit m p photons in the same mode with a fidelity of F m p 1 .However, for each atom that transitions from |g to |s , the single-photon collective decay rate from e → g is reduced by g , whereas the collective decay rate from e → s is increased by s , which leads to imperfect impedance matching.Using Eq. ( 9) to estimate the resulting additional reflection, we find the probability for absorbing the kth photon is reduced by 1 − 2k 2 /3N 2 .Thus the probability that m p photons are absorbed is reduced by a factor of 1 − 2(2m 3  p − 3m 2 p + m p )/3N 2 to third order in 1/N.Combining this with the single-photon fidelity F 1 , we can estimate the m p -photon QND fidelity as In the following, we compare these predictions with numerical simulation and find they agree well.

(a) (b) (c) (d)
FIG. 4. QND number-resolving detection.(a) Overlap F 1 of the reflected wave packet with the incoming wave packet when a single photon in a Gaussian mode of width g is scattered from a perfect array initially in state |ψ 0 = |N/2, N/2, 0 , obtained from a numerical simulation (dots) for Purcell factors of P = 10 (blue) and P = ∞ (green, no free-space decay).Free-space decay clearly reduces the fidelity, as is expected, but does not affect the scaling with the number of atoms, as the coherently enhanced coupling grows faster than the incoherent decay.The fidelity is captured very well by the transmission fidelity of an impedance-matched two-port cavity in which each port has decay rate κ = N g /2 [black line, see Eq. ( 23)].In this effective description, free-space decay corresponds to internal decay of the cavity as outlined following Eq.( 23).(b) The fidelity as a function of the initial imbalance m 0 = m − N/2, again comparing numerics with the analytical calculation for a two-sided cavity [Eq.(23)].Here and in the following we always take P = 10.(c) Numerical calculation of the fidelity F m of the QND measurement like in (a) , but now for several photons.The dots are numerically calculated, the solid lines is a simple estimate, given in Eq. ( 24) below.(d) A time trace of a simulation of five photons in the same Gaussian mode of width g scattering off 40 atoms.This illustrates that transmission is good even with modest atom numbers, and that the probability for an atom to be in the excited state remains small throughout, implying that free-space decay is suppressed.The numerical method we use is detailed in Appendix D.

C. Numerical simulation
In order to verify these conclusions numerically, we study the scattering of a multi-photon Fock states with a Gaussian wave packet of width g from an array of atoms in the atomic mirror configuration in the mirror geometry using a recently proposed technique [73].We present details of the simulation in Appendix D.
In Fig. 4(a), we show the fidelity for a single photon scattering off an array starting from the initial state |N/2, N/2, 0 .Clearly, Eq. ( 23) is a good approximation to the transmission fidelity of the atomic array.This is still true for any other symmetric starting state, as we illustrate through Fig. 4(b), which shows the fidelity of single-photon scattering when starting with an initial state |N/2 − m 0 , N/2 + m 0 , 0 , with the imbalance m 0 ranging from −N/2 to N/2.We assume there is a maximally achievable decay rate max g , s and consequently lower either g or s to fulfill the above condition (22).Together, these results show that single-photon transmission is captured very well by the above equations.
Turning to multiphoton scattering, we simulate several photons in the same Gaussian wave packet scattering off the atomic array and again calculate the overlap of the output wave packet with the input wave packet.The results are shown in Fig. 4(c).We find that our simple argument captures the fidelity well, and that our proposal can in principle reach very high fidelities for modest atom numbers.Finally, in Fig. 4(d), we show an example of the time evolution of the system.

D. Dicke state preparation
Following similar arguments as in the other sections, the QND detector is robust against spatial disorder.However, the suppression of free-space decay crucially relies on the preparation of a Dicke state between the two ground states.
One way to obtain such a state is to start with all atoms in the ground state |g , apply a π/2 pulse on the ground states {g, s}, and finally perform a projective measurement of S gg (or S ss or S gg − S ss ).This heralds a fully symmetric state with a binomial distribution of imbalances around zero and standard deviation √ N/4.However, such measurements are difficult.In principle, one can apply an off-resonant probe in the waveguide and recording the phase shift of the reflected light, which is proportional to the number of atoms [74].In Appendix B, we briefly analyze this kind of measurement and find it is fundamentally limited to atom numbers of the order of the Purcell factor N P. It has been proposed to produce atomic states by manipulating the dark-state manifold [75], which however is limited in fidelity by 1 − F ∝ N/(2 √ P), where N is the number of atoms and P the Purcell factor.Neither method scales well to many atoms.Thus, in the following, we find a fast preparation method that imposes much less stringent requirements on the Purcell factor.
We propose to produce Dicke states in a way that is in keeping with the core idea of this article: measuring atoms, not photons.This requires two arrays (halves of one array) that are individually addressable with external driving.We further require that the atomic transition frequency be in the band gap of the waveguide, which allows the atoms to be coupled coherently without dissipation.As shown elsewhere [32], such a setup readily gives rise to a Hamiltonian that couples all spins where the combination of coupling to individual waveguide modes g k ≈ const, bound state decay length L = √ α/ , detuning of the impurity from the band edge = ω 0 − ω b , and FIG. 5. Dicke state preparation with Eq. (25).A symmetric Dicke state can be prepared by coupling one fully excited array to one in the ground state according to Hamiltonian Eq. ( 25).(a) Time evolution of magnetization of the two arrays, N = 100.The initial state is recovered after a period π/g eff .After a time π/2g eff the state of each array is close to |N/2, N/2, 0 with small fluctuations in the imbalance.The first equilibration happens much faster, after g eff t 0 = 1/ √ N, which is beneficial to reduce free-space decay.(b) Dependence of t 0 on atom number N. (c)Imbalance distribution at time t 0 for N = 100 atoms.band curvature ω k = ω b + αk 2 yield the effective coupling strength g eff = 2π g 2 k /( L).Importantly, when using a Raman transition (s → e → g) to couple the atoms to the waveguide modes, both and g k in Eq. ( 25) are tunable.Careful analysis (cf.Appendix C) reveals that the effective Purcell factor P eff = g eff / eff,free scales as the ratio of the effective waveguide density of state (1/ √ α) and the constant free-space density of state (ρ 0 ), viz., P eff ∝ 1/(ρ 0 √ α ).The detuning from the band edge can in principle be made arbitrarily small without violating the adiabatic condition or the Markov approximation (cf.Appendix C).Note also that disorder in the positions of the atoms is not an issue here, as there is no position-dependent phase. 4Instead this scheme is likely ultimately limited by disorder in the coupling strengths or energy, as they destroy the symmetry of the effective Hamiltonian.
The protocol to prepare an approximately half-excited state of one of the arrays is as follows.One array is fully excited (by applying a π pulse) and the other is left in the ground state.The time evolution under the Hamiltonian (25), shown in Fig. 5(a), transfers excitations from one chain to the other, while leaving them in their individual symmetric subspaces.Notably, after a short time g eff t 0 ∼ 1/ √ N, corresponding to the first zero crossing in Fig. 5(a), the average number of excitations in each array is equal.However, this comes with the caveat that while on average the two arrays hold N/2 excitations each, in fact their imbalance has a very wide probability distribution, as we illustrate in Fig. 5(c) for N = 100 emitters.As we have shown above, a large imbalance does not invalidate our scheme, but it does mean that the usable atom number on average is halved.Since this is a constant penalty, it does not change the overall scaling with atom number N. Ultimately, this scheme requires P eff √ N, where P eff is the Purcell factor enhanced through the proximity to the band edge.
There is another, intrinsically faster, way to prepare Dicke states, if the system is governed by the Hamiltonian H = g eff S (1)  sg S (2)  gs + S (2)  sg S (1)  gs . ( As we detail in Appendix C 2, this Hamiltonian may be engineered through the interference of interactions via two different waveguides (or bands in the same waveguide), and constitutes a waveguide QED version of spin flip-flops recently realized in cavity QED [76].While certainly more challenging to implement experimentally, this Hamiltonian has the advantage of equilibrating the number of excitations in each array on an asymptotic timescale of g eff t 0 ∼ ln(N )/(2 √ 2N ) (cf.Appendix C 3).This is the fastest time that can be achieved for a given g eff , essentially saturating the timescale obtained by adding all average transition times N n ||S (2)  + S (1)  − |N, n ⊗ |N, N − n || −1 ∼ ln(N )/2N.As a result, using the Hamiltonian Eq. ( 26) reduces the requirement on the effective Purcell factor to P eff ln(N ).Another advantage of this Hamiltonian is illustrated by Fig. 6(c), which shows the probability distribution of Dicke states after a time t 0 .The overall probability for the state to be close to |N/2, N/2 is larger as with the other Hamiltonian [cf.Fig. 5(c)].

V. EXPERIMENTAL CONSIDERATIONS
A tacit assumption in the preceding sections has been that the decay rates from |e to |g and/or |s are tunable.In circuit quantum electrodynamics, decay rates can be tuned  87 Rb level scheme for photon detection.Detunings i and drivings i are chosen such that the excited states can be eliminated.(a) Destructive photon detection.Any of the decays into the green shaded region from | f 2 is intended and the total decay rate via the green channel forms eng .Since the decay 1,g is superradiantly enhanced, we have 1 2 .The deleterious decays shown in red are analyzed in the main text.(b) Nondestructive photon detection.We have refrained from drawing in all decay channels that contribute.In comparison with the superradiantly enhanced decays f 2 → s and f 1 → g, they scale as O(1/N ) and thus can be suppressed by increasing the number of atoms.by changing the detuning of an intermediate resonator [77].In atomic systems, this can be done with Raman transitions, which we analyze in the following for the D 2 line of 87 Rb.For concreteness, we assume here that the waveguide efficiently couples to π transitions.

A. Engineered decay: Destructive photon measurement
The level scheme we consider specifically is drawn in Fig. 7(a).In it, we make the choice |g ≡ |F = 2, m F = 2 , |e ≡ |1, 1 , and {|s i } = {|2, 1 , |2, −1 , |1, 0 , |1, −1 } (green area), which coupled via excited states (in While this is just one choice among many it has the advantage that the applied lasers do not couple to any transition of the large number of atoms in the ground state |g .Since in the destructive scheme we only need to count how many atoms have scattered photons, the final state is irrelevant, provided it is not |g .The Raman scheme allows for large tuneability of the relative decay rate, which is required to obtain e→s ( 2 ) ≈ N e→g ( 1 ).
To show that the additional decays drawn in red do not spoil the scheme, we derive the effective quantum master equation governing the double-model in their presence.
Neglecting the energy shifts due to the pumps, the dynamics are purely dissipative, given by the jump operators where in the last expression 1,g 1 = 1,g and 2,g 2 = 2,s .We denote the (sum of the) rates corresponding to these jump operators e→g ( 1 ), e→s ( 2 ), and ee,eff , respectively.We note that the deleterious decays i,e induce dephasing described by Lee,ee that is not negligible, primarily due to the i = 2 term, as 2 1 .However, the dephasing removes the coherence and thus the superradiant decay to |g .Since e→g / e→s = O(1/N ), the corresponding excitation decays to |s , with an error O(1/N ).We conclude that all potential errors analyzed here are suppressed with increasing N.

B. Engineered decay: nondestructive photon measurement
Similar considerations apply to the operation of the QND detector.A possible choice, shown in Fig. 7 Unlike the destructive scheme, now all other decays are deleterious.However, both decays e → s and e → g are enhanced by a factor of N/2.This implies that photon loss, either into free space or through decay into another hyperfine state scales as 1/N.All deleterious decay rates can be taken together to form free in the calculations of Sec.IV.

C. Dicke state preparation
In order to prepare Dicke states, we propose to start with one array of the atoms in one hyperfine ground state |s and the other array in another hyperfine ground state |g .In this protocol, the requirements for the Purcell factor (ratio of coupling strength to free-space decay) is most stringent.The most favourable level scheme is therefore a closed -system, where |g = |F = 2, m F = 2 and |s = |2, 1 are hyperfine ground states, but | f = |3, 3 is an excited state in 5 2 P 3/2 .In this case, f → g is a transition and decays outside this subspace are suppressed.Coherent driving between |s and | f can be implemented using a two-photon transition [78].This leaves free-space decay as error source, which has been discussed in Sec.IV D.

D. Other sources of disorder
We have neglected inhomogeneous broadening and atomic motion, which could be taken into account in the same way as positional disorder and do not modify our conclusions.Furthermore, the relative effect of disorder decreases as the number of atoms increases [22].On the other hand, fast atomic motion, analyzed in Appendix F, essentially only renormalizes the coupling strength of the atoms to the waveguide (thereby decreasing the Purcell factor), but is not a fundamental obstruction.
In typical experiments today, the atom number might fluctuate in unknown ways from one experiment to the next.Allowing some error in N is like increasing the detuning error in eng .As long as the relative error decreases with N, as it should, this does not affect the scaling with atom number.

E. Photon loss from coupling into the waveguide
It is challenging to couple photonic wave packets traveling in free space into fibre.If this is done with 80% fidelity, say, then the detector becomes already unreliable for three photons, as there is already a 50% chance of losing one in the first stage.
However, the detectors described here are also useful for entirely waveguide-based setups, which obviate the need for coupling free-space wave packets into a fibre.Recent proposals show that atomic arrays make excellent sources of quantum light [46][47][48] especially in view towards quantum metrology.There also exist proposals for two-photon gates in this platform [49], and it has been found that such arrays coupled to waveguides support long-lived subradiant states that can be used for storage [50,79] and also manipulation [80] of quantum states of light.

F. Atom readout
After the photons have scattered, one needs to measure the number of atoms in one of the states, or ideally both.Readout of superconducting qubits is well studied due to the advent of quantum computers.Fidelities now reach >99% in realistic devices [43,44].
For neutral atoms, a range of techniques have been developed.Among the earliest was to employ a cycling transition between a ground state and some excited state, which in recent experiments has yielded fidelities of 98%-99% [81,82].A similar idea is to push one type of atoms out of the trap with a resonant drive and use a quantum gas microscope to detect occupied sites [83].Instead, one could "boil" one type of atom away through inelastic light scattering.This can in principle achieve extremely high fidelities (99.97%), but is limited by atom loss and has the drawback that the lattice has to be refilled.We note that detection of atom numbers has already been used in experiment [33,84].Likewise, a strong magnetic field [85] or a state-dependent trap [86] can be used to separate the states and image them afterwards.In these setups, imaging has been performed with 99.94% fidelity, but background gas collisions still cause atom loss [86].One option to prevent this would be tweezer arrays [87].While such capabilities have not yet been demonstrated for atomic arrays coupled to waveguides, and scattering from the fibre might make photon collection harder, a solution could be to move the optical lattice slowly away from the fibre before performing the imaging.
After successful detection, the detector can be reset by pumping the s → e transition, such that eventually all atoms decay to |g .

VI. CONCLUSION
We have explored the use of arrays of quantum emitters coupled to waveguides for number-resolving photon detection.Paying particular heed to experimental limitations such as disorder and free-space decay, we have found that both can be overcome, leaving no fundamental limitation to the achievable detection efficiency.Moreover, we have shown that the same platform can also be used to perform QND numberresolving photon detection.To this end, we also propose a novel way to prepare Dicke states based on the interaction of two arrays and subsequent heralding by measuring the number of excitations in one of them.
In a nutshell, our proposal builds on four facts that together enable highly efficient detectors: (1) few-level systems allow for strong, projective measurements of their state due to their intrinsic nonlinearity, (2) nevertheless, sufficiently large ensembles of atoms are linear, (3) collective decay mitigates errors due to nonidealities, and (4) in linear systems, one can always engineer dissipation to obtain complete absorption.We hope that the ideas outlined here will mark a step towards high-fidelity number-resolving photon detectors, both destructive and nondestructive.
This solution can be plugged into the equation of motion for σ (n) ge .We need the following integral The equation of motion for σ (n) ge becomes We define the input field a in (t ) as the portion of the waveguide field that was at a position x = ct at time t = 0 and has since traveled all the way to the atoms.Thus a in (t ) = − √ c/2e iω 0 t a(ct, 0), where the prefactor is fixed by the commutation relations of a in , up to an arbitrary phase.This yields the Langevin equation As defined, a in (t ) is a slow variable, so if the dynamics of the system and the bandwidth of the input state around ω 0 are slow compared to the time it takes for light to travel a distance 2x n /c, we can neglect the retardation, rendering our description Markovian.The same applies to the atomic lowering operators.Finally, we arrive at a time-local equation In order to calculate the output field, we take the integrated equation of motion for the light field and apply a sine transform.This is essentially the same as the right-hand side of the equation of motion for σ (n)  ge , except evaluated at a different point in space.Choosing this point to be x R + ε, i.e., a small distance to the right of the rightmost atom, and again neglecting retardation, we find Further choosing ε such that sin[k 0 (x R + ε)] = 1 and defining a out (t ) = − √ c/2e iω 0 t a(x R + ε, t ), we have Finally, let us define the decay rate g = g 2 /2c.

Three-level systems coupled to two semi-infinite waveguide fields
In the main text, we consider (effective) three-level systems with state |g , |e , |s coupled to two waveguide fields, which respectively couple to the transition g ↔ e and s ↔ e. Starting from the Hamiltonian the derivation follows through as above, and yields and two more equations if 1 ↔ 2 and g ↔ s are exchanged simultaneously.These equations reduce to Eq. ( 18) for atoms in the atomic-mirror configuration k 0 x n = 2π (n + 1/4).

Two-level and three-level systems coupled to one and two infinite waveguide fields
The derivation for the infinite waveguide proceeds in much the same way and can be found elsewhere [13].A Hamiltonian that combines both bath couplings reads Here, k L,i are phases imparted by the laser.They have little effect on the detection of incoming light.There are two waveguide fields, α = 1, 2, distinguished either in frequency, polarization, or by being in a different waveguide.As before, ν ∈ {±1} labels right-and left-moving modes.From Eq. (A16), we can derive the bath equations of motion, integrate them up and Fourier transform them [13] a ν,1 (x, t ) = e iω 0 t a ν,1 (x − ct, 0) For the other field, we have to exchange 1 ↔ 2 and g ↔ s.As above, we next derive the atomic equations of motion and replace the photon field [Eq.(A17)] x n σ n ge , (A19d) Using the description above, the Langevin equations for two-level systems on an infinite waveguide read whereas the governing equations for three-level systems are and two more equations if 1 ↔ 2 and g ↔ s are exchanged simultaneously.

APPENDIX B: PREPARING A DICKE STATE WITH COLLECTIVE MEASUREMENTS
There is another way to prepare a Dicke state in our setup.First applying a π/2 pulse to all atoms, and then measuring S gg , i.e., the number of atoms in |g , which projects the state onto a symmetric state with a definite number of excitations.For large N, m has variance √ N around N/2, such that |N/2 − m| N. A π/2 pulse can for example be produced by driving the two transitions g ↔ e and s ↔ e with lasers, each detuned by a large amount , or by applying a microwave tone.After the pulse, the atomic state becomes A known method to measure the number of atoms is by measuring the phase shift of an off-resonant probe tone [74].In order to prevent excited atoms from decaying via emission of a free-space photon (essentially removing the photon from the symmetric state) or via emission of a photon in the e → s channel (thereby inducing a transition to another Dicke state), the probe has to be far detuned.However, the further detuned the probe is, the lower is the induced phase shift, thus requiring longer averaging times, which leads to a trade-off.The analysis below shows that this sort of measurement is limits the number of atoms to the Purcell factor.
If it is possible to turn off s , we only have to consider free-space decay.Using the input-output equations, we find the output field for a weak coherent probe with amplitude α in detuned by N g g where N g is the number of atoms in state |g .Thus the per-photon phase shift is ϕ − g /2 .The number of photons required to resolve single excitations therefore is N p ∝ 4 2 / 2 g .On the other hand, for a given probe amplitude α in the probability of an atom to be excited is S ee g N g |α in /2 | 2 .Thus during the time it takes to scatter N p photons, free N g / g = N g /P photons are lost.Ultimately, this way of preparing Dicke states is thus limited to atom numbers lower than the Purcell factor.If s cannot be set to zero, but instead is comparable to g , then this constitutes the dominant decay channel and the requirement on the Purcell factor becomes even more stringent P > N g N s .

Strong coupling
Here we give details on how using a Raman transition allows strong coupling of the atoms when brought close to the band gap of a waveguide and discuss briefly some limitations.
Following [32], two-level systems interacting via the modes in the band gap couple at a rate g eff = 2π g 2 k / √ α, where g k is the coupling to the individual waveguide modes (assumed constant), is the effective detuning of the renormalized transition frequency from the band edge and α is a parameter to characterize the band curvature at the band edge, via ω k = ω b + αk 2 .If g k is the effective coupling rate obtained from adiabatically eliminating an excited state |e in a Raman transition, it takes the form g k = g k,0 ( /δ), where is the pump Rabi frequency and δ is the detuning of the pump from the transition.Thus In order for the Markov approximation and the adiabatic elimination to be valid, we require δ and √ α g k , so g eff will be slow.Fortunately, reducing /δ also reduces free-space decay in the same way, free,eff = free ( /δ) 2 .
Since the atom-atom coupling additionally scales with the effective waveguide density of states g eff ∝ (α ) −1/2 , the Purcell factor P eff = g eff / eff,free ∝ (α ) −1/2 can be made very large, independent of the original Purcell factor, while preserving the Markov condition g.As this analysis shows, physically this relies on increasing the effective waveguide density of states.At the same time, this has the effect of making the bound state extent and therefore the decay length of the induced interaction L much larger than the extent of the atomic array, such that the all-to-all interaction on the right-hand side of Eq. ( 25) becomes a very good approximation.
In order to instead realize a Hamiltonian with the interaction S (1)  + S (2)  − + H.c. between two atomic arrays, one needs to combine two bands of the same waveguide (or two waveguides), and additionally make use of the spatial profile of the induced atom-atom interaction.The full expression for the coupling induced between the atoms through a waveguide band is given through [32] which differs from Eq. ( 25) through the addition of the spatial profile E k 0 (x) of the bound state induced by coupling an atom at position x, where k 0 is the wave vector at the band gap.A simple model of a waveguide exhibits bands with a dispersion ω k ∼ cos(ka).As such, there are two band edges at wave vectors k 0 = 0 and k 0 = π .In more complex models the band gaps might occur at different values of k 0 , but these are generically not expected to all coincide, and certainly not in between two different waveguides.With an otherwise constant mode profile, we can approximate E k 0 (x i ) ≈ e ik 0 x i .For simplicity, we will take k 0 = π , but this is by no means required.
Combining two waveguides (waveguide bands), one with positive detuning , the other with the opposite detuning − , and placing both arrays of atoms in atomic mirror configurations, but spaced by an odd half-integer-multiple of wavelengths from each other, the two waveguides induce the interaction S sg S gs − S (1)  sg − S (2)   sg S (1)  gs − S (2)   gs = 2g eff S (1)  sg S (2)  gs + H.c. , (C3) as required.Note that we have taken L to be larger than the whole configuration, and choosing g eff to be the same for both waveguides, which is not unreasonable as they can be tuned with the Raman transition.Furthermore, we would like to note that the requirement of placing the atoms in two atomic mirror configurations spaced by an odd multiple of half a wavelength is not any more stringent that producing a single atomic mirror configuration, although it certainly is experimentally more challenging than the disordered configurations we have considered in other parts of the main text.

Interaction timescales
While the full Hilbert space of two arrays of N spins each is 2 2N dimensional, the Hamiltonian Eq. (C3) connects the fully symmetric state |N, 0 = |N 1 ⊗ |0 2 , to only N other states, {|N − m, m }.As in the main text, our convention here is that |m denotes a fully symmetric spins state in which m out of N atoms are excited, whereas |m 1 , m 2 denotes two arrays of N atoms each, with m 1 and m 2 excitations, respectively.
In this space of N + 1 states, the Hamiltonian (C3) is a matrix of the form where As an aside, we note that this is similar to the model studied in Ref. [88], except with the elements of the matrix squared.
The connection arises since the model with √ a m on the diagonal maps to the N-excitations subspace of two coupled harmonic oscillators H = a † b + H.c. (leading to perfect state transfer), while here we instead study two coupled spins H 2 = 2g eff S (2)  + S (1)  − + H.c. If the index m = x is understood as a spatial coordinate, the Hamiltonian H 2 can be rewritten as ), (C6) where [x, p] = i, such that e i p generates a translation by −1.In order to get some insight into the dynamics of H 2 , we consider the classical long-wavelength limit neglecting terms pn with n 3 and replacing x → x, p → p, which yields where This Hamiltonian gives rise to periodic trajectories.A more physical Hamiltonian would be obtained by sending p 2 → −p 2 and H → −H, but the dynamics are the same.Starting at x(t = 0) = p(t = 0) = 0, the conserved energy of the particle is H (0, 0) = E 0 = 2g eff (3N − 1).This allows us to solve for the momentum as a function of position C9) which can be used to calculate the period of the orbit The upper limit of the integral is the point at which the particle turns around, which can be found from p(x) = 0.For large N, the integral may be approximated through For completeness, we mention another way to arrive at this result.Instead of Eq. (C3), consider the Hamiltonian where a, b are the annihilation operators for two harmonic oscillators.In the subspace spanned by the states |2N − 2m 1 ⊗ |2m 2 , where now {|m } are now Fock states of an oscillator, the Hamiltonian is again given by a matrix of the same form as Eq.(C4), except that now (C13) For large N, m these are essentially the same matrix elements as before, and one can check numerically that the dynamics is very similar for large enough N. In the classical limit, we replace the annihilation operators by the amplitudes of the corresponding coherent states.The classical mean-field equations of motion read The shape of the resulting periodic orbits can be found by integrating dr/dφ.
However, since the phase φ is ill-defined in the initial state, we can choose it to be φ = π/2.In this case, the system is governed only by which is readily integrated to yield Since initially, r/K 1 − 1/N and K N, we conclude that the natural timescale for the switch is the same result as before.We note that if another phase φ is chosen initially, φ first rapidly evolves to a value close to π/2, after which the system evolution is very similar.

APPENDIX D: NUMERICAL SIMULATION OF MULTI-PHOTON SCATTERING
In this Appendix, we briefly outline the theory behind our simulation of multiphoton scattering, following Ref.[73].The key idea is to the time evolution of the system and bath for a specific input field using the Langevin equation ( 18) in combination with the input-output equations (19).
In general, an input wave packet can be described through one or more modes of the waveguide in which photons are created on top of the vacuum |0 , |ψ in = j dx j ψ j (x j )a † (x j )|0 . (D1) As a function of time, the field travels through the waveguide, and eventually the photons interact with the atoms locally.Instead of this spatiotemporal description, in which a quantum system couples with constant rate to the waveguide, but the wave function of the photons changes in space, it has been shown that the dynamics can be equivalently captured by modulating the coupling between the system and the input mode according to the mode shape [73].Specifically, here we assume that all m p photons of the input field are in one mode, such that the waveguide field can be written (at t = 0, before the photons interact with the system) Since this is assumed to be the wave packet before the interaction, it has support only for negative x (we assume the wave packet is incident from the right in Fig. 1).After a time t 0 , the center of the wave packet (travelling at speed c) has arrived at the system.
Replacing the spatiotemporal evolution of wave packets in the waveguide with a time-dependent coupling constant, the input-output equations can be replaced by a quantum master equation [ describes the dissipation induced by the coupled waveguides.We determine the time evolution given by Eq. (D4) using QUTIP [89].Note a subtlety when comparing with the equation corresponding to Eq. (D5) derived by Kiilerich and Mølmer [73], which features an additional term g * out (t )g in (t )a † i a o .This term arises when calculating the scattering from one input mode to the corresponding output mode, which in our case might be a in,1 → a out,1 and can be thought of arising from the first term in the input-output equation Eq. (19).Here, we are interested in the scattering a in,1 → a out,2 , where this term is absent.
In the above expressions, the variable couplings strengths are given by g in (t ) = ψin (t ) Here, the initial time to start the simulation has been chosen to be t = 0, and since photons are assumed to travel at speed c, the temporal wave function is related to the spatial wave function through ψin (x) = ψ in (−ct ).For consistency, g in (t ) = g out (t ) are taken to be zero before t = 0, and ψ in should be square-integrable for positive times ∞ 0 dt |ψ in (t )| 2 = 1.This is true for our choice Eq. (D3) as long as −1 g t 0 .Since t 0 is an arbitrary offset it can always be chosen to fulfill this condition.
If the system and the input field are initially in states of known excitation number, which we assume to be the case throughout, then it is clear that for m p photons in the input wave packet, we need to consider at most m p excitations in the bosonic modes a o and a i as well as the system.Recall that we assume a symmetric Dicke state as starting point, |ψ sys (0) = |N/2 − m 0 , N/2, 0 .As a result, the required Hilbert space has a dimension of (m p + 1) 3 .
We note that this formalism does not specify what shape the output wave packet ψ out has.Indeed, it is generically not even true that it can be described through a single mode and neither is it true in general that all photons are emitted in the second channel rather than the first.The advantage of this formalism in our case is that the system dynamics are independent of the choice of ψ out (t ).Instead, if the set of output modes considered does not comprise all modes of the physical output field, the corresponding photons are lost.In the extreme case, when g out (t ) = 0, the QME Eq. (D4) describes the system dynamics for a given input, but does not yield any information about the output photonic state.
Here, we use this property by setting the output mode equal to the input mode ψ out (t ) = −ψ in (t ).This is equivalent to asking the question how many photons are transmitted without changing the shape of the wave packet-i.e., performing a QND measurement.If the final number of excitations in mode a o is equal to the number of input photons m p , all photons have been transmitted faithfully, and all photons have been dissipated via the second channel.This defines the fidelity F, which we plot in Fig. 4. Mathematically, it is defined as where |m p is the m p th Fock state of the output mode a o .This is the probability that the all photons have been transmitted in the specified output mode.

APPENDIX E: EFFECT OF FINITE PURCELL FACTOR ON PERFORMANCE OF QND DETECTOR
In this Appendix, we offer a side-to-side comparison of the fidelity with a finite Purcell factor, as shown in Fig. 4 in the main text, and infinite Purcell factor.This allows to discern which features come from free-space decay, and which from nonlinearities and finite bandwidth.It furthermore offers further verification of the predicted analytical approximate fidelity Eqs. ( 23) and ( 24) through simulations.The results are shown in Fig. 8 and are commented on the caption.

APPENDIX F: FAST THERMAL MOTION OF ATOMS
In the main text, we have only considered static disorder, which is valid for slowly moving atoms.If the thermal motion of atoms is fast, some of the effect of disorder will be averaged out.This can be modelled by instead averaging the atomic coupling over a distribution of atomic positions [78].FIG. 8. Comparison between presence (blue) and absence (green) of free-space decay.Overall, the effect of free-space decay on the fidelity is clearly captured by introducing the branching ration as we have done in Eq. ( 23).This brings down the single-photon fidelity considerably.Interestingly, though, as we consider multiphoton scattering, the effect of free-space decay diminishes relative to the effect of the nonlinearities introduced.Naturally, it never aids the fidelity, but in the end, the fidelity retains its favourable scaling with atom number, most clearly illustrated in Fig. 4 where F (x) is the purely real Dawson integral.It is peaked at x = 1, is odd, and obeys |F (x)| < 0.6, and F (x) → 0 ± as x → ±∞.Equation (F4) predicts that the coupling decreases exponentially in (k 0 σ ) 2 .For low to moderate k 0 σ , the effect of fast thermal motion can simply be accounted for by re-scaling the couplings, without affecting the conclusions in the main text.

FIG. 3 .
FIG.3.Photon detection with an infinite waveguide.(a) Detection probability for a disorder-free array of N atoms with different spacings.This graph illustrates the fact that an array in the atomic-mirror configuration (a = λ) cannot serve as photon counter, whereas all other generic spacings work similarly well, a fact that can be understood analytically (cf.Sec.III B).(b) Probability for an undetected photon on resonance for an atomic array coupled to an infinite waveguide as a function of atom number.The blue line denotes the limit of a perfectly ordered array with spacing a = λ/4 (or 5λ/4 etc).Purcell factor is P = 10, the average for the other cases was performed over 2500 disorder realizations (standard deviation shown as lightly coloured area).(c) Scaling of largest eigenvalue in ordered arrays with varying spacing compared to a fully random one (red).While the atomic mirror configuration (a = λ) is clearly different, it is a fine-tuned exception, with all other generic arrays (ordered or disordered) behaving remarkably similar.The robustness of our scheme relies to a large degree on this universal eigenvalue scaling.

FIG. 6 .
FIG.6.Fast Dicke state preparation with Eq. (26).(a) Time evolution of magnetization of the two arrays, N = 100.The time dependence can be approximated by Rabi oscillations for a few periods, which can be described in the continuum limit as outlined in Appendix C 3. (b) Dependence of the zero crossing time t 0 on atom number N. (c)Imbalance distribution at time t 0 for N = 100 atoms.

FIG. 7 .
FIG.7.Double-Raman87 Rb level scheme for photon detection.Detunings i and drivings i are chosen such that the excited states can be eliminated.(a) Destructive photon detection.Any of the decays into the green shaded region from | f 2 is intended and the total decay rate via the green channel forms eng .Since the decay 1,g is superradiantly enhanced, we have 1 2 .The deleterious decays shown in red are analyzed in the main text.(b) Nondestructive photon detection.We have refrained from drawing in all decay channels that contribute.In comparison with the superradiantly enhanced decays f 2 → s and f 1 → g, they scale as O(1/N ) and thus can be suppressed by increasing the number of atoms.
which could optionally be written as|X = 2 −N/2 m ( N m )|N − m, m, 0 .Measuring Ŝgg thus probabilistically returns the state |N − m, m, 0 where m follows a binomial distribution.

2 ×
Assuming a Gaussian distributions of positions around the atomic mirror configuration, captured by the random variable y [e ik 0 |y m −y n | + e ik 0 (y m +y n ) ]. (F2) If m = n, there should only be one integral, giving ḡnn = −[1 + exp(−k 2 0 σ 2 )] g /(4g).In the case m = n, we can straightforwardly evaluate the second term, which yields − g e −k 2 0 σ 2 /4 overall.For the first term, we first shift y m → y m + y n , in which case the y n integral becomes a straightforward Gaussian giving a factor of √ πσ 2 e y 2 m /4σ 2 .The leftover

πσ 2 dy m e −y 2 m /4σ 2 e ik 0 |y m | = − g 4 e −k 2 0 σ 2 [ 1 +
FIG. 8. Comparison between presence (blue) and absence (green) of free-space decay.Overall, the effect of free-space decay on the fidelity is clearly captured by introducing the branching ration as we have done in Eq. (23).This brings down the single-photon fidelity considerably.Interestingly, though, as we consider multiphoton scattering, the effect of free-space decay diminishes relative to the effect of the nonlinearities introduced.Naturally, it never aids the fidelity, but in the end, the fidelity retains its favourable scaling with atom number, most clearly illustrated in Fig.4(a).