Topologically protected strongly correlated states of photons

Hybrid photonic nanostructures allow the engineering of novel interesting states of light. One recent example is topological photonic crystals where a nontrivial Berry phase of the photonic band structure gives rise to topologically protected unidirectionally-propagating (chiral) edge states of photons. Here we demonstrate that by coupling an array of emitters to the chiral photonic edge state one can create strongly correlated states of photons in a highly controllable way. These are topologically protected and have a number of remarkable universal properties: the outcome of scattering does not depend on the positions of emitters and is given only by universal numbers, the zeroes of Laguerre polynomials; two-photon correlation functions manifest a well-pronounced even–odd effect with respect to the number of emitters; and the result of scattering is robust with respect to fluctuations in the emitters’ transition frequencies.


INTRODUCTION
Light-matter interaction is around us in nature and has played a tremendous role in the development of current technology.Until recent decades it was sufficient to deal with this interaction on average, with many photons and many atoms involved.However, increasing miniaturization of basic constituents towards nano scale is a common trend of modern technology.Downscaling to single-atom and/or single-photon levels promotes some traditionally classical research areas into the quantum realm [1], [2] [3], [4].Therefore control over quantum processes of light-matter interaction will eventually be a vital ingredient of emerging quantum devices and is as well important for the rapid developments in several related fields, including communication, signal processing, ultrafast optics, optomechanical cooling, imaging and spectroscopy and, of course, quantum information.Efficient manipulation and control however requires a relatively strong interactions on a level of single atom, single photon and/or a single electron levels [5], [6], [7], [8], [9], [10], [11].This represent a significant challenge since a typical interaction scale for individual particles in our everyday life is given by the QED coupling constant α ≈ 1/137.Two possible ways to overcome this natural limitation and to increase the effects of correlations is either to use artificial materials and devices or to engage many-body interacting effects which would eventually lead to desired non-linearities in systems with reduced dimensionality.
Recent tremendous experimental progress in fabricating few-photon sources coupled to 1D transmission lines, [12], [13], [14], [15], [16], [17], [18] opens a prospect for creating and manipulating strongly-correlated states of photons.This also triggered significant number of theoretical studies [19], [20], [21].Indeed, our experience in condensed matter and atomic physics teaches that often combined effects of reduced dimensionality of a quantum system and inter-particle interactions may effectively enhances the effects of correlations and eventually may lead to a new, collective states of matter with properties that are very different from those of the underlying particles (one notable example is a Luttinger liquid state of interacting fermions or bosons in 1D).This kind of philosophy is one of the driving idea of a quest for the novel correlated states of photons [22].
FIG. 1: Schematic view of our proposed system to generate strongly-correlated states of photons: the topological photonic insulator (bulk of the left picture) has topologicallynontrivial band structure (at the right picture).If the total Chern number of the bands below the gap is 1, the state in the gap is formed (red line in the right picture).This state corresponds to the chiral edge state of unidirectionally propagating photons (thick red line at the boundary of the photonics crystal).We suggest to insert emitters of any internal level structure at the edge state (here only two-level structures are shown).Pumping the edge state from aside by a few-photon source (in-state) creates strongly-correlated photonics states at the outcome of the structure (out-state).
One of these collective phenomena is a topological insulating/superconducting state.Recent progress in experimental construction of topological insulating states for electrons [23] has generated a bunch of experimental and theoretical activity in search for the novel topological states of matter [24].A central signature of the quantum state to have a topological property is an exis-tence of a chiral edge state, a unidirectionally propagating waves/particles which is insensitive to the local perturbations (such as imperfections, local potentials, etc.).The whole existence of these states is due to non-trivial topological properties of the band structure in the bulk of a material.
Recently a unidirectionally propagating chiral states of electromagnetic waves have been suggested theoretically in [25] and observed experimentally in [26] which generated a subsequent large theoretical and experimental activities [27], [28], [29].These states are an optical analogue of the quantum Hall edge states.It is known that Hall effects belongs to a broader class of topological insulators exhaustively classified recently [30].These topological states are characterized by nontrivial Chern number in the bulk and chiral topologically protected edge states.In 2D photonic topological insulators the edge states form a 1D waveguide with unidirectional propagation of photons, while in 3D topological photonics insulators the edge state is a 2D surface state with Dirac structure of a spectrum [31].Topological insulators for photons can be constructed by engineering specific magnetooptical photonic crystals producing a band structure with nontrivial Chern numbers.
Here we suggest to use the unidirectional edge states (chiral edge state in 2D topological photonic insulators) as a robust platform for controllable generation of strongly-correlated photonic states.To achieve this goal we suggest to couple an array of emitters to the edge state.Photons in the edge state populated by an external few-photon source interact with an ensemble of emitters and produce an outgoing photonic states with robust, controllable and universal properties.An origin of these robustness and universalities lies in the topological nature of the edge state.We found in particular that the outgoing photonic wave function has no information about position of emitters while its zeroes are given by the universal numbers -zeroes of Laguerre polynomials.Even on a few-photon level the wave function has "fractionalized" form and position of one zero imply the others thus making the outgoing state strongly correlated.We also observe strong even-odd effect with respect to the number of emitters corresponding to switching between bunching and anti bunching behavior for even and odd number of emitters respectively.We show also that all these properties are robust with respect to effects coming from temperature, imperfections and/or disorder.These correlated states of photons can be used in emerging quantum devices.

1D EDGE STATES WITH EMITTERS
We consider unidirectionally propagating topologically protected chiral edge electromagnetic states in photonic crystal.We assume that an array of two-or multi-level systems can be inserted at that chiral edge channel.We assume that transition frequencies of that systems (which we call emitters) are close to the frequencies of the propagating chiral mode.Moreover, since the origin of the chiral edge comes from the Dirac-like cone structure in the k-space, it is reasonable to assume a linear dispersion for the propagating chiral model.Therefore, our model has a form of electromagnetic waves with a linear dispersion, propagating in one-dimensional edge geometry and interacting with a system of two-or multi-level emitters.The Hamiltonian can be written therefore as Here, the electromagnetic field operators a † (x), a(x) satisfy a usual commutation relations [a(x), a † (y)] = δ(x−y) and have a linear dispersion spectrum, whereas the emitter-related variables S a b are placed at coordinates r a at the 1D edge.Emitters are determined by the internal level structure and by the allowed transitions between them.Thus, if we have a two-level system then a = x, y, z and S ± = S x ± iS y ) satisfy the usual spin algebra of atomic operators [S z , S ± ] = ±S ± , [S + , S − ] = 2S z and we have a natural representation of the transitions in the two-level system with transition frequency given by ∆ b for the emitter at position r b .For the array of emitters we can allow the frequencies ∆ b to be distributed around some average frequency ∆.On the other hand, if we have a three-level structure of the emitter, the definition of operators S a depends on the type of a scheme, which can be of either types, Λ, Σ, or V .Thus for the Λ-type scheme we have S + = g 31 |3 1| + g 32 |3 2| while for the V -scheme S + = g 31 |3 1|+g 21 |2 1| and for the Σ-scheme S + = g 32 |3 2| + g 21 |2 1|.In all cases S − = S † + .In this paper we mainly focus on the two-level case.
The problem is further specified by the initial state.Here we assume a natural initial state inspired by the current experimental realizations of waveguides with emitters coupled to it [16], [17]: we assume that the edge state is pumped in by the external few-photon source while the emitters are initially in the ground state.Injected photons unidirectionally propagate as a wave packet at the edge state and interact with the ensemble of emitters.The resulting state is read-off by the detector.In the basis of a Fock states |in = dk i f ({k i }) i a † ki |0 where i = 1 . . .N and f ({k i }) is an envelope function of momenta {k i }.The evolution of this initial wavefunction is given by the scattering matrix S(N ; M ; {l i }), |out = S(N ; M ; {l})|in where the scattering matrix S(N ; M ; {l}) describes a scattering of N -photon state |in on the array of M emitters each of which may have a different number of levels characterized by a set of num-bers {l b }, l b = 2, 3, 4, . . .for b = 1 . . .M .In the case of two-level emitters the model (1) can be solved exactly by the Bethe Ansatz and the scattering states for any initial conditions can be determined exactly [33].Basic ideas of this approach are summarized in the Supplement.To solve a more general problem with combination of emitters of different level structures we are using a combination of Bethe ansatz [32], [33], [34], [35] and a scattering matrix formalism [36] (see also Supplement).A combination of these methods provides a complimentary pieces of insight into how the combination of several basic ingredients of a model (position of emitters and their level structure, strength of light-emitter interaction, etc.) serves as a powerful tool for engineering correlated multiphoton states.
Various results of our computations suggests the following fundamental property of the model: the scattering matrix of the N -photon state on the array of M atoms is given by the convolution of the N -photon scattering matrices on every atom, where the convolution ( * ) can be taken in any convenient representation (e.g.coordinate or momentum spaces).Sufficient conditions for this property to be correct are the following: (i) unidirectional nature of a spectrum, and (ii) a constant group velocity of the incoming wave packet.Indeed, in the case of a chiral spectrum the transfer matrix is identical with the scattering matrix provided that the group velocity is a constant for all particles.As a first straightforward consequence we conclude that the outgoing state of photons does not depend on the positions of emitters.This is the origin of many universal properties of the outgoing photonics states.[42] Few-photon initial conditions We discuss scattering of a few-photon wave packets on arrays of M emitters.If we measure all frequencies in units of κ (thus e.g.∆ phys = κ∆) and the distances in units of κ −1 ( e.g.y phys = y/κ), the coupling constant enters the problem exclusively through this renormalization of frequencies and distances.

Single-photon scattering
Using the formula derived in [34], we investigate the scattering of a chiral photon on array of M atoms of transition frequencies ∆ b , b = 1 . . .M .The formula for the out-going wavefunction produced by the δ(x)-function instate reads where z = x + t, t is time.Here First we consider the limit ∆ j → ∆, which can be performed analytically (see Supplement).The calculation yields − κΘ(z − y)e i∆(y−z) L (1) where L M −1 (x) is an associated Laguerre polynomial.Assuming Gaussian form of the incoming wavefunction 2 , where δ is the detuning w.r.t. the atomic transition frequency ∆, the outgoing wavefunction can be written as In These results manifest clearly a strongly-correlated and universal character of the scattering in our system.The zeroes of the Laguerre polynomials determine the positions of minima of the outgoing wave packet.Thus if the position of the first peak is known the subsequent peaks and minima are universally determined by the zeroes of L (1) M −1 (x).The result does not depend on positions of M emitters and thus universal.Time delay on every emitter formed by the emitter-photon bound state in (4) (exponentially decaying part of the wave function) leads eventually to "fractionalization" of the wave packet into M − 1 peaks.

Two-photon scattering
The two-photon scattering can be shown to have three parts: (i) the one corresponding to no interaction with emitters, (ii) the one corresponding to one photon interacting with emitters while a second one goes through the ensemble without interaction and (iii) the irreducible part of scattering corresponding to two photon-bound state formed in the presence of emitter.This picture of different scattering channels will be important below to understand the outcome of the scattering.
Let us look into the Gaussian incoming wavefunction of the form   where δ is again a detuning from the atomic resonance frequency.We are mainly interested in the limit of wide pulse µ → ∞, however the relative position of the two photons is confined into the interval of the order of σ.
When the limit ∆ b → ∆ is taken (see Supplement) we observe again a polynomial structure related to Laguerre polynomials.Moreover, the scattering in the case of δ = 0 does not depend on the number of atoms, only of its parity.These limiting forms of the outgoing wavefunction reads for even M even, and for odd M .Here d = y 1 − y 2 is a photon's relative coordinate.One can also calculate the asymptotic expansion of the outgoing wavefunction for large detuning δ.The result up the first nontrivial order δ −2 reads FIG. 3: Two photon scattering on M -atoms.Note that the detuning δ and the width of the pulse σ are measured in units of κ, i.e. δ phys = κδ and σ phys = κ −1 σ.Note, that while for M even the photons do not scatter at δ = 0, for M odd they do.Note that for large detuning δ the effect of scattering asymptotically vanishes.
There are several effects present in the scattered state off M emitters: one is a renormalization of the incoming Gaussian waveform, which is likely not to be detectable, the other is the exponential tail, which scales as δ −2 regardless of the number of atoms.This is another universal feature of the scattering in our system.The exponential tail is apparently a remnant of single-particle-only scattering, which can be seen either by direct comparison with 1-particle scattering formulas, or by employing the scattering matrix in k-space.There one can easily see that the irreducible two-particle part of the scattering scales as a higher power of δ than the one-particle scattering.Further, the outcome of the two-particle scattering is sensitive to the parity of the number of emitters: for the odd number of emitters a tendency towards anti bunching is pronounced while for the even number of emitters outgoing photons are bunched, see Figs (3).The polynomial nature of the scattered wave function persists for the two-particle scattering as well.

Robustness of correlated states
Next we consider a more physical situation when atomic frequencies are random variables distributed according to the Gaussian distribution.The source of the randomness can be manyfold.In particular it can be caused by the temperature since for existing topological photonic crystals operating in the GHz regime the temperature effects can be important.This could also be imperfections or external noise coupled to emitters.We will model the broadening of the levels by the random Gaussian ensemble with a variance Σ, so that the distribution of levels ∼ exp(− ) is taken around some mean detuning δ.Based on exact expressions we can perform the ensemble averaging of the observable quantities.In Figs.(4) we plot typical results for the averaged two-photon probability density corresponding to g 2 correlation function.We observe robustness of the results we described above for the big range of variations of parameters.We conclude that our results are robust with respect to level's broadening.

DISCUSSION
We suggested a novel approach of producing stronglycorrelated states of photons in a highly controllable way.Our setup uses the edge states of photonics topological insulators with multi-level emitters coupled to it.We observed a number of universal features of scattering in our setup: multi-particle scattering does not depend on position of emitters.The outcome of the scattered wave packet has a polynomial structure and the minima of the outgoing pulse are given by zeroes of this polynomial.For a single-particle (and reducible part of multiparticle scattering in general) case these polynomials are Laguerre polynomials.In a single-particle case the scattering thus looks like if the photon would be "fractionnalized" in between the minima of that polynomial.
Physical realization.To observe quantum many-body effects of photons described here one needs a photonic topological insulators introduced and studied in [26], [27].The role of emitters can be played by the quantum dots or superconducting qubits, like e.g. in [10] and [15].Also, available single-photon emitters made of quantum dots coupled to photonic crystal waveguides makes our system realizable at the today's level of development [11], [18], [8].Typical frequency ranges of existing photonic topological insulators are GHz with matches the photonics single emitters frequencies and couplings.Temperature effects relevant for this frequency range show up in the inhomogeneous broadening of levels.We take this into account by assuming Gaussian fluctuations of detunings for every emitter.Photon's initial state injection can be made using the same single-photon emitters.
Extension for a general level's scheme.Here we focused on a two-level structures, however using the convolution property 2 one can immediately apply several recent results on scattering off three-or four-level schemes [20,21], which contain a number of interesting experimentally observable properties [7,37] including singleemitter-EIT, transistor-like behavior, implementation of quantum gates (e.g.CNOT [38]), etc.By combining emitters of different level structure this flexibility makes it possible to build an optical schemes with desired correlated properties of outgoing photonic wave function.
For explicit evaluation of outgoing state one should use a scattering formalism [36] (see also Supplement) together with the convolution property (2).
Extension for different initial conditions.Here we explicitly evaluated the properties of the scattered light for photonic initial conditions, when emitters are not excited initially.Other initial conditions can be treated similarly.In particular one can think of emitters initially excited from aside.In this case the evolution can be studied using similar method while the outgoing wave function will depend on the emitter's positions.Also coherent state initial condition can be analyzed using our experience with one-and two-photon scattering as well as a knowledge of scattering on a single emitter.If we focus on the case of a single-mode k 0 coherent state |α k0 than it becomes clear that a multi particle irreducible parts make a minor contribution to the outcome |out k0 in momentum state k 0 .One can draw this conclusion by expanding the coherent state into a Fock basis and using the observation that multi-particle irreducible part has smaller available phase space than the reducible parts.Observed universalities of a single-particle scattering will be therefore present in the coherent state initial state.More detailed arguments are given in the Supplement.
Effects of dissipation of photons into outside world is equivalent to opening of additional scattering channels in our approach.Although detailed studies of decoherence of topologically-protected states are not available to our knowledge, based on experimental observation of longlived edge states [26], [27] one can generically believe that topologically protected states will be robust with respect to dissipation.
Currently available topological photonic crystals operate at GHz frequency range.This is related to difficulties of creating large magneto-optical response (which is needed for nontrivial topology of the band structure) with existing materials.However some recent developments [39] in material science may lead to significant magneto-optical response in the optical domain which then may be used to engineer topological photonic crystals in the optical range of frequencies.
Finally we note that in 2D interacting photons can lead to anyonic statistics and fractional Hall states [40].Recent suggestion of realizing Dirac cone structure [31] may pave the way to this exotic physics of interacting 2D photonics fluids.Moreover recent classification of interacting bosonic topological insulators [41] suggest a possible route to engineer even more exotic states of light.
We expect that the novel correlated states of photons or plasmons will find a number of applications, some of them we can only envision now: quantum information technology, optomechanics and precision measurements. where Here z j ≡ x j + t and the integration contour Γ is determined as follows: in the complex plane of parameters λ the contour is a composition Γ = γ 1 ⊕ γ 2 . . .⊕ γ N such that: (i) Imγ j+1 -Imγ j > κ for j = 1, . . ., N and (ii) at Reγ j = ∆ a one has Imγ j > −κ/2 for x j < r a .Here we assume that initially photons are injected in the most left region of space with respect to all emitters and that all emitters are initially in their ground states.One can also consider more general initial conditions, when some (or all) emitters are initially in the excited states.In this case the condition (ii) for the choice of the contours γ j should be modified by adding two more: Imγ j < κ/2 for x j > r a and −κ/2 < Imγ j < κ/2 for x j = r a .Note that here we corrected the definition of contour Γ with respect to Ref. [34].
Let us consider several important examples.

Single-photon scattering
First we consider a scattering of a single chiral photon state off an array of two-level systems.The evolution of the plane wave state |k = dxe ikx E + (x)|0 ≡ E + (k)|0 can be computed using the general scheme of the evolution of the basis state.In this case the amplitude φ 1 (y) can be computed [34] as (this result is also quoted in the main text) where The out-state of the wave is therefore given by where (repeating the steps of section III.A of [35]) which reproduces known results for the chiral part of scattering state [19].It is a convolution (product) of scattering phases on individual emitters [36].

Two-photon scattering
By definition, the two-particle scattering matrix S (2) M scatters the in-state (2) For the case of scattering off an array of M emitters the scattering matrix reads

S
(2) M (y 1 , z 1 )S M (y 2 , z 2 ) + S (1) M (y 2 , z 1 ) + iT (2) where S (1) M (x, y) are the single-particle scattering matrices.The first two terms in the last line represent the reducible part of the scattering (independent single-particle scattering) while T (2) M (y 1 , y 2 , z 1 , z 2 )-matrix represent correlated irreducible part.The reducible part of the two-photon scattering off M emitters can be written as In case of resonance δ = 0 we obtain S M ;red (0; ∆y, ∆z) = It is also convenient to define
The irreducible two-photon part of the scattering is given by Note that for the calculation of this quantity the mutual arrangement of the contours γ 1 and γ 2 does matter.Following the above prescription for a choice of integration contours one can show that (cf.Eq. ( 21) from Ref. [34]) A special case of M equal detunings is interesting.Let us consider the limit ∆ a,b → ∆ of ( 26) for z 1 > z 2 > y 2 > y 1 .At first, we introduce an integral representation which helps to decouple the sums over a and b.At second, we employ the identity [cf.Eq. ( 18)] With its help Eq.( 26) is transformed in this case into Integrating over τ we obtain This expression is remarkable, since it resembles Eq. ( 26).However, instead of the summation over distinct discrete detunings ∆ a and ∆ b , we have now the double integral over continuous variables s and r.The second term exactly compensates the contribution coming from the pole s = r − iκ in the first term, so one can rewrite (30) as (31) where the contours γ 1 and γ 2 are small circles clockwise embracing the points s = ∆ − iκ/2 and r = ∆ − iκ/2, respectively (cf.Eq. ( 44) and (45) of Ref. [33]).
Restoring a complete form of Eq. ( 31) by adding the terms (y 1 ↔ y 2 ) • (z 1 ↔ z 2 ), we transform it to a mixed representation of the full momentum where the component T (2)++ M has ∆y > 0 and ∆z > 0, and the other ones are obtained from it by flipping the signs of ∆y and ∆z.Quantatively T (2)++ M (δ; ∆y, ∆z) is where both contours γ 1 and γ 2 embrace the origin clockwise.Noting that we cast (33) to From the identity it follows that where P M −2 (r) is a polynomial of the degree M − 2. One can then show that Inserting it into (35), we obtain where X = ∆y + ∆z.The formula (39) is capable to account for the parity effect in the case δ = 0. We have Combined with the result for the irreducible part this leads to explicit formulas for the parity effect described in the main text.
At finite δ we represent (39) as The function F (δ/κ, κX) is defined by where the function can be replaced by the polynomial because of condition s = 0 after differentiation in (42).One then finds that where C M −1 l are binomial coefficients.We therefore conclude that the outgoing wave function for the two-photon case develops a polynomial form as well.
Plots for the probability density of the two-particle scattering are shown in the main text.Here we supplement them for different set of parameters.
FIG. 6: Two photon scattering on M -atoms.Note that the detuning δ and the width of the pulse σ are measured in units of κ, i.e. δ phys = κδ and σ phys = κ −1 δ.For zero detuning δ = 0 the contributions from the one-particle and the two-particle scattering cancel each other and as a result the wave goes through unscattered.Note that for large detuning δ the effect of scattering asymptotically vanishes.

SCATTERING FORMALISM
Here we describe an application of the scattering formalism developed in [36] for computing the nontrivial elements of the scattering S and T matrices for N -photonic states on arrays of two-and tree-level emitters in various positions.The advantage of this approach is that it has no limitations of the Bethe ansatz.In particular, one can consider an array of emitters of different level structure (two-level, three-and four-levels) for different emitters and for arbitrary emitter-dependent coupling strengths κ a .The approach agrees with the Bethe ansatz in the limit κ a → κ for the two-level scheme.
Single-photon scattering.A generalization of (17) for the scattering of a single photon off M emitters for different coupling constant is obtained by multiplying corresponding phase shifts on individual emitter, An agreement in the limit κ a → κ is obvious.Single photon scattering matrices off 3-level emitters are listed in [36] and agree with known results [20].
FIG. 7: Two photon scattering on M -atoms.Note that the detuning δ and the width of the pulse σ are measured in units of κ, i.e. δ phys = κδ and σ phys = κ −1 σ.For zero detuning δ = 0 the contributions from the one-particle and the two-particle scattering cancel each other and as a result the wave goes through unscattered.Note that for large detuning δ the effect of scattering asymptotically vanishes.
Two-photon scattering.The scattering matrix (21) off M emitters is a convolution of M scattering matrices off individual emitters, as suggested in [36].Let us illustrate this on example of the two emitters M = 2.
For different κ 1 and κ 2 the solution is presented in [36].Here we focus on κ 1 = κ 2 = κ.First we take the Fourier transform of the Bethe ansatz solution ( 26) and evaluate it for the case of M = 2.We have Omitting the factor 1 π δ(p 1 + p 2 − k 1 − k 2 ) we evaluate the sum in (46) where Note that the expression (48) is invariant under the exchange of ∆ 1 and ∆ 2 (which also implies α 1 ↔ α 2 and M 1 ↔ M 2 ).
Using the identities and we cast (48) to Furthermore, we notice that and conclude that T which agrees with the result of convolution [36] evaluated at κ 1 = κ 2 .This illustrate the convolution property (2) from the main text.

SCATTERING OF COHERENT LIGHT
Recently we studied the scattering of the initial coherent state off a single emitter [36].Our results suggest the following picture of many-body effects in the coherent state.The probability of finding n k particles in the mode k before scattering is given by the coherent state Poisson distribution.
where coefficients s are hereby defined in such a way that the state |m; l ∈ k with l-particles in mode k are normalized to 1.
In fact, the particles leave the mode k only if they scatter in an i-particle irreducible way with i ≥ 2, since only then the photons can redistribute the total momentum between themselves and the chance to stay in k is measure zero.A sole 1-particle scattering conserves momentum of a photon and hence the photon cannot leave the mode k.The coefficients s(m; l) therefore measure amplitude of relative weight of 1-particle scattering to more particle irreducible scattering.E.g. s(m; 0) is generated by purely 1-p scattering (1-PIS), s(m; 1) = 0 because of conservation of momenta, s(m; 2) is contributed only by 2-PIS, while s(m; 4) is contributed by both 4-PIS and combination of two 2-PIS.
Hence we can qualitatively understand the behavior of the probability distribution P (n k ) after scattering.Indeed For large γ the more-PIS dominate over 1-PIS and s(m; l) is small for small l.Morover, the states with low n k acquire contribution from all higher lying states.Therefore, as a result, the weight shifts to low n k which implies that few-particle effects are mostly important.

FIG. 2 :
FIG. 2: One photon scattering off M emitters for detuning δ = 0 (upper panel) and scattering off M = 10 emitters.Initial state has a form of Gaussian wave packet with dispersion σ.For other values of σ see the Supplement.All frequencies are expressed in the units of κ and δ phys = κδ and σ phys = κ −1 σ.

FIG. 4 :
FIG.4: Averaged two-photon scattering off M emitters.The transition frequencies ∆a are independent random variables distributed from the Gaussian distribution of variance Σ.The detuning δ is the deviation of the incoming wave frequency from the sample mean frequency.The solid lines correspond to the ensemble average of the density (second-order correlation function).The filled regions denote the median (darker) and mean (lighter) absolute deviation.The detuning δ, the width of the pulse σ and the variance Σ are measured in units of κ, i.e.δ phys = κδ and σ phys = κ −1 σ.

FIG. 5 :
FIG.5:One photon scattering on M -atoms for detuning δ = 0 and the initial Gaussian wave packet (see the main thext) characterized by σ.