Large-N limit of spontaneous superradiance

We investigate the thermodynamic limit of Dicke superradiance. We find an expression for the system's density matrix that we can prove is exact in the limit of large atom numbers N. This is in contrast to previously known solutions whose accuracy has only been established numerically and that are valid only for a range of times. We also introduce an asymptotically exact solution when the system is subject to additional incoherent decay of excitations as this is a common occurrence in experiment.


I. INTRODUCTION
Superradiance was introduced by Dicke in 1954 [1].In a simple model of all-to-all coupled emitters, he showed that during their decay process, quantum coherence is established spontaneously and leads to a superradiant burst of photons in which the maximum decay rate scales as the square of the number of atoms Γ max ∼ γ 0 N 2 /4, where γ 0 is the decay rate of an individual atom.Superradiance can be construed as synchronization of many emitters [2], akin to the onset of lasing.Moreover, since quantum effects are typically difficult to observe in large ensembles of atoms, Dicke's observation stirred a lot of interest in the 70s [3][4][5][6][7][8][9], as reviewed by Gross and Haroche [10].Superradiance was observed experimentally, also in the 70s [11][12][13][14].A closely related, but different, phenomenon is the Dicke phase transition [15,16], which we will not consider here.
In light of these developments, we re-visit the theory of superradiance.Approximate solutions have been derived in a number of different ways [5][6][7][8][9][10], but ultimately they all rely on a continuum limit in which the magnetization can take noninteger values.This elegant description is, however, incorrect at short and long times.Additionally, it sensitively depends on a chosen initial distribution, and thus needs a rigorous jus-tification.Since Dicke superradiance is a fundamental model of quantum optics, it is remarkable that no solution exact in the limit N → ∞ has been derived so far.
The main contribution of the present work is an explicit formula for the density matrix as a function of time, alongside a proof that it is asymptotically exact as N becomes large.In contrast to previous work, our solution works for all times.In a second contribution, motivated by experimental realities, we derive the solution in the presence of incoherent loss (for example through decay into free space rather than guided modes) and show that it is exact as N → ∞.

II. SETUP
We consider the quantum master equation of N two-level systems ("atoms" or "emitters") with states |0 , |1 that are subject to both collective and incoherent decay, where σ − n = |0 n 1|, S − = m σ − m , and D[a]ρ = aρa † − (1/2)(a † aρ + ρa † a).Note that we have omitted the transition frequency of the atoms from the description, as it does not affect they decay dynamics.In Eq. ( 1), the dimensionless parameter γ controls the ratio of incoherent to collective decay.For γ = 0, the model originally studied by Dicke is recovered (to be distinguished from the "Dicke model" introduced by Hepp and Lieb [15]).
Any permutation of the atoms is a symmetry of Eq. ( 1).Thus, a permutation-invariant initial state remains so, i.e., ρ(t) = Π s [ρ(t)], where Π is the permutation operator and s is any permutation of the atoms.Defining the z-component of the collective angular momentum S z = n σ z n and S 2 = S + S − + (S z ) 2 , we can label collective states of N spin-1/2 particles according to their total angular momentum j and its z-component m.The Hilbert space of N two-level systems is then spanned by states {| j, m, α }, which obey S 2 | j, m, α = j( j + 1)| j, m, α and S z | j, m, α = m| j, m, α , where 0 ≤ α < d j is an integer accounting for their multiplicity d j = N!(2 j + 1)/((N/2 − j)!(N/2 + j + 1)!).
We would like to solve Eq. (1) when the system is initialized in the (permutation-invariant) state corresponding to fully excited atoms ρ(0 The original Dicke model is obtained by setting γ = 0 [1], and we call this situation "pure superradiance".In this case, the system explores only the Dicke states with maximum total spin |n = |N/2, −N/2 + n , where n ∈ {0, 1, • • • , N}.The density matrix of the system is diagonal at all times, ρ(τ) = n P n (τ)|n n|, and thus Eq. ( 1) reduces to a rate equation Here, for convenience we rescaled time τ = Nt, and defined The most salient feature of this model is the superradiant burst that occurs at τ ≈ ln N during which most excitations are emitted at a rate that scales with N 2 .

B. Dicke model with incoherent decay
Free-space decay competes with superradiance, as it reduces the coherence of the ensemble and allows the system to explore states with lower total spin.For convenience, we label these states by their number of "dark" excitations r = N/2 − j, which measures by how much the total spin has been reduced, and the number of collective excitations n = N/2 − 2r + m that can be removed from the state before it reaches the bottom of the Dicke ladder.In terms of n and r, the relevant states are projectors into the space of defined n and r, where in terms of r the multiplicity reads d r = N!(N − 2r + 1)/(r!(N − r + 1)!).As before, Eq. ( 1) reduces to a rate equation transitioning between these states, and the solution is a mixture ρ(τ) = n,r P n,r (τ)P n,r [62], where the probabilities P n,r obey ∂ τ P n,r = −Γ (1)  n,r P n,r +Γ (2)  n+1,r P n+1,r +Γ (3)  n+2,r−1 P n+2,r−1 +Γ (4)  n,r+1 P n,r+1 , (5) and the rates are given through [63] Γ (1)  n Here and in the following, we distinguish distributions in n, r by using calligraphic font.Note that collective decay (terms without γ) cannot change the total spin j (and neither r), and its rate is reduced by the presence of dark excitations r.All other terms, proportional to γ, are due to incoherent loss, which can either reduce r (Γ (4) ), increase r (Γ (3) ) or leave it unchanged (Γ (2) ).

C. Previous work
Several approaches have been used to describe pure Dicke superradiance.In principle, one can write down an exact iterative solution [64], but it cannot be summed up in closed form to yield a formula for the magnetization for general N. Approximate solutions for P n (t) were derived in different ways by Degiorgio [5], Degiorgio and Ghielmetti [6], and Haake and Glauber [7].All these solutions are equivalent and read with λ = 1.Gross and Haroche later provided another derivation, but quote P with λ = 0 [10].The choice of λ has little effect and in the following we take λ = 0.All approaches rely on the continuum limit of Eq. ( 2) (dropping the 1 in γ n to be consistent with the choice λ = 0), where n is interpreted as a continuous variable n ∈ [0, N].
Other authors arrive at the same conclusion by matching it to the corresponding Wigner distribution [7], or to predictions from the master equation [6,27].However, the approximation through the continuum equation ( 8) is not valid at short times and its results have not rigorously been shown to be correct [66].As a result, the distribution Eq. ( 7) is incorrect at short times.While the magnetization Eq. (9a) agrees well with numerics at finite N, to our knowledge it has not rigorously been shown that it is valid in the limit N → ∞, not even at sufficiently long times.The case of γ 0 has received considerably less attention.The competition between collective and incoherent decay has been studied in some detail [67], and a stochastic unravelling was used to efficiently simulate the system [68], but a complete theory is missing.We note that in the context of the Dicke phase transition in cavity QED, a number of works have found a non-trivial interplay between dephasing, incoherent loss, and collective coupling [69][70][71].

D. Our contribution
Our first contribution is to derive a solution for pure superradiance that we rigorously prove to be correct in the large-N limit for all times (Theorem 1).This is in contrast to the previously established solution, Eq. (7), which, as we show, is asymptotically exact only for intermediate times.Specifically, we show that the one-norm of the difference between the exact probability vector P * and the literature solution P obeys lim but lim For α > 2, || P|| 1 < N α−2 π 2 /6 → 0 as N → ∞, which implies || P − P * || → 1. Equation (10) follows from the fact that for constant times, the normalization of P tends to 1 (Lemma 2).Since the exact solution P * is normalized, the reverse triangle identity gives Eq. (10).
Our second contribution is to extend our analysis to γ 0, where we again provide a solution that is valid for all times (Theorem 2).

III. LARGE-N LIMIT OF PURE SUPERRADIANCE
Our goal is to find an asymptotically exact solution R n in the sense that the one-norm of the difference to the exact solution P * vanishes as N → ∞, like in Eq. ( 11), but we require that the solution works for all times.We can prove that this is fulfilled by R given through where τ 1 = (1 + δ 1 ) ln N with δ 1 = 2/5, and There is a slight subtlety, since the time of the superradiant burst occurs depends on N, and thus it is not sufficient to prove that the solution converges at a constant time.Instead we allow for general sequences of times that may depend on N.
At short times, we replace γ n → γ 0 n = N − n + 1 in Eq. ( 2), which is valid since the probability distribution is expected to have support only for N − n N. The resulting equation is solved by Using this solution, we show that the error due to the substitution γ n → γ 0 n vanishes as N → ∞ as long as we only consider times up to τ 0 = δ 0 ln N with 0 < δ 0 < 1.We then use Q(τ 0 ) as the initial condition for Eq. ( 8) and thus obtain a solution up to τ 1 .Since Q n (τ 0 ) is sufficiently smooth in n, we can show that the error incurred by using the continuum approximation vanishes as N → ∞.A further approximation of the result yields R < [Eq.(13a)].For late times, we again linearize the equations of motion, γ n → γ 1 n = n, which is valid as the probability distribution is expected to concentrate around n ≈ 0. This yields the solution R > in Eq. (13b).We can bound the error to show that ||R > − P * || 1 → 0.
A few comments are in order about the solution we provide.First, for any finite N, R neither exactly fulfils the equation of motion Eq. ( 2), nor the continuum equation Eq. ( 8), and its derivative is discontinuous at τ 1 , although R itself is continuous.Nevertheless, R describes the correct distribution at all times.Second, as part of our proof, we show that ||P * (τ) − P(τ)|| 1 → 0 as N → ∞, which proves that the solution found in the literature also converges to the right solution, but only for intermediate times.Third, a remaining question is whether the magnetization and radiance (9a) are accurate, since they are calculated by integrating the continuum distribution rather than summing the (true) discrete distribution.With the following corollary, we show that both formulae are nevertheless asymptotically exact for all times.
Corollary 1 (Magnetization and radiance).The magnetization and radiance (9a) predicted in the continuum limit are asymptotically correct for any sequence of times One might wonder whether a universal form for the magnetization and radiance can be found that is independent of N and thus holds in the large-N limit.We can obtain such a form by using the continuum limit of Eq. (7).Performing a change of variables to T = e τ /N, we obtain Apart from the normalization 1/N, this expression is independent of N. The corresponding magnetization µ = 1/2+ S z /N and radiance ρ = S + S − /N 2 read (plotted in Fig. 1a) Note that Eqs. ( 17) and ( 18) are identical to Eqs. ( 7) and (9a) up to the change of variables, but we find it more convenient to analyze the behaviour at N → ∞ in this form as they are independent of atom number N, having absorbed all dependence on N into T .Note that the time variable T (t) = exp(Nt)/N starts at T (0) = 1/N and exponentially quickly moves to infinity.The superradiant pulse occurs close to T = 1 and, as is evident from Eq. ( 18), its height scales with N 2 .We can determine the maximum radiance numerically and find it is ρ ≈ 0.196 at T ≈ 1.39, at which point µ ≈ 0.532.In normal units, the time of maximum radiance is t pulse = ln(N)/N + 0.330/N.At t = 0 we have T = 1/N, at which point µ(1/N) ≈ 1 − 1/N.Thus, the starting point of the dynamics corresponds to the time at which on average one photon has decayed.Another result from the above analysis is that in the original time t, the magnetization tends to a step function µ(t) → Θ(ln(N)/N − t pulse ).

IV. SUPERRADIANCE IN THE PRESENCE OF INCOHERENT LOSS
Since incoherent loss competes with superradiance and is virtually unavoidable in realistic settings, an important question is whether it may preclude superradiance in certain regimes.We answer this question in two parts.First, we derive the exact solution as N → ∞, which shows that in the large-N limit, superradiance always persists (Sec.IV A).Second, the question then arises whether there is a critical emitter number that is needed to observe superradiance.This is of experimental relevance particularly in systems in which γ, the ratio of free-space decay to waveguide decay, is large.We answer this by calculating the threshold emitter number N above which signatures of superradiance emerge as a function of γ (Sec.IV B).Finally we also compute how many photons decay into free space on average (Sec.IV C).

A. Asymptotically exact solution
The asymptotically exact solution R(τ) that converges to the exact solution P * of Eq. ( 5) at all times is given by where τ 1 = (1 + δ 1 ) ln N with δ 1 = 2/5, and where F is a Poissonian distribution and We are interested in the one-norm of the difference of R and the exact solution P * and again allow for arbitrary sequences of times.
Theorem 2 (Superradiance with incoherent loss).For any sequence of times {τ n > 0 : n ∈ N}, The proof closely follows the structure of Theorem 1 and can be found in Appendix B 2. To unpack this result, we highlight a few important properties of superradiance in the presence of incoherent loss below.

B. Qualitative behaviour at large N and threshold
At large atom numbers N, the average magnetization and radiance are given by Eq. (9a) up to subleading errors.While correct, this description misses the fact some excitations are transferred into states that do not decay collectively.The population of such dark (or subradiant) excitations grows linearly at a rate γ until the superradiant burst, such that on average, about γ ln N such dark excitations are produced.Since they can only decay incoherently, they dominate the late time behaviour with a slow decay at a rate γ/N.
The result that magnetization and radiance behave as in pure superradiance does not necessarily apply to finite N and in fact there is a threshold N below which superradiance disappears.One necessary requirement for superradiance is that the maximum radiance occurs at t > 0. We thus can calculate the threshold above which this is the case directly from the quantum master equation ( 1) by evaluating

C. Number of incoherently decayed photons
Since it is clear that superradiance persists in the limit N → ∞, we compute the number of photons lost incoherently during the whole process, N loss , based on the assumption that the system undergoes pure superradiant decay, which becomes exact as N → ∞.Since the incoherent loss occurs at a rate equal to γ times the number of excitations in the system, we integrate the solution Eq. ( 18) to obtain where γ E 0.577 is Euler's constant, and the approximation becomes very good for N > 10.Note that Corollary 1 also applies to Eq. ( 24) as it is just the integral over µ(τ).The leading-order behaviour N loss = γ ln N, valid when ln N 1 N loss /N, can be obtained from the simple observation that the superradiant burst occurs at t ln N/N and that the emitters are mostly excited before and mostly in the ground state after.The fact that N loss ≈ γ ln N implies that as N → ∞, the fraction of atoms that decay incoherently is vanishingly small.This is ultimately the reason why in the large-N limit, the predictions from pure superradiance for magnetization and radiance also hold in the presence of incoherent loss.

V. OUTLOOK
Our work improves previous solutions to Dicke superradiance and puts them on a rigorous footing.There are a number of other experimentally relevant effects that should be considered in future work, some of which may fundamentally change the behaviour at large N. Typical effects include disorder in the decay rates of atoms into the collective mode, inhomogeneous broadening, finite temperature, and, specifically in the case of waveguide QED, spatial disorder.All of these distinguish the atoms and therefore make a straightforward extension of the theory presented here difficult.While inhomogeneous broadening becomes negligible at large N, it is not clear at present how the other contributions would affect our solution.

ACKNOWLEDGMENTS
We would like to thank Nina Fröhling and Tao Shi for insightful discussions.We acknowledge funding from ERC Advanced Grant QUENOCOBA under the EU Horizon 2020 program (Grant Agreement No. 742102), and within the D-A-CH Lead-Agency Agreement through Project No. 414325145 (BEYOND C).RT acknowledges Max Planck Harvard research center for quantum optics (MPHQ) postdoctoral fellowship.
Appendix A: Proof of Theorem 1 We prove Theorem 1 in a series of lemmas following the outline given in the main text (see Appendix A 3). Before, we establish some definitions and auxiliary results.

Definitions
Definition 1 (Exact solution P * ).We denote the exact solution by P * (τ) = (P * 0 (τ), • • • , P * N (τ)).It obeys the equation with γ n = n(N − n + 1)/N, which can equivalently be written as the matrix equation Definition 2 (Early solution Q).We denote the solution at early times by with γ 0 n = N − n + 1, which can equivalently be written as the matrix equation Definition 3 (Continuum solution P).We denote the continuum solution by P(τ) with components which, for continuous n ∈ (0, N) obeys We take P to denote the vector formed by taking integer n, Note that Definition 3 is the solution found in the literature.
Definition 4 (Our solution R).We use A choice of parameters consistent with all constraints is (A4c)

Additional lemmas
To bound the differences between the various probability vectors, we need a few additional results.In the following, || • || always denotes the 1-norm.
First we show that in the limit N → ∞, the probability mass of the distributions Q(τ) and R(τ) in the interval from n = 0 to N − N µ vanishes faster than any polynomial of N. We will use this on many occasions to restrict the range of sums over n.
Lemma 1 (Vanishing probability mass).If τ < τ 0 = δ 0 ln N and µ > δ 0 , To get to the second line, we use that the maximum is obtained at n = N − N µ and we bound e −τ ≤ 1 and 1 − e −τ ≤ 1 − e −τ 0 .Similarly, Lemma 2 (Normalization of P).For constant times τ, the normalization of P obeys Proof.We would like to evaluate || P|| 1 = N n=0 Pn for constant τ.First note that if n ≤ N − N µ for any µ > 0, we have Pn < (N 2 /n 2 ) exp[−e −τ N µ ] → 0 and thus we can restrict attention to n > N − N µ .We define s = N − n and evaluate where a = exp(−e −τ ) is a constant with 0 < a < 1.In the limit N → ∞ the geometric series gives e −τ /(1 − a), which establishes the result.
We frequently need to bound sums over P for times in the range τ ∈ [τ 0 , τ 1 ], τ 1 = (1 + δ 1 ) ln N, for which we use the following definition.
Proof.First, note that we can write (A10) Let us bound τ ≤ ln N and τ > ln N separately.For this we define and For τ 1 > τ > ln N, the maximum with respect to n is reached at (A14) When τ ≤ ln N, the maximum is reached at n max > N. Since n can be at most equal to N, we instead obtain Combining these results, we arrive at Eq. (A9).

Main proof
To prove Theorem 1, we prove the statement separately for short, intermediate, and late times.

(A19)
For a given τ, the term inside the absolute value signs in Eq. (A19) changes sign at n = N + 1 − 2e τ .Thus we can split the sum up into two parts, one from n = 0 to n = n = N + 1 − 2e τ , and the other from n = n + 1 to n = N − 1. ∆ 1 obeys the inequality where with x = a, b, and N a = {0, . . ., n} and N b = {n + 1, . . ., N}.
In each part, all terms have the same sign (+1 in the first, −1 in second), so we can get rid of the magnitude sign such that only the boundary terms survive.
Considering first ∆ b , we have (A22) The right-hand side of Eq. (A22) is an increasing function of τ, so we can replace τ by τ 0 , which yields Clearly this vanishes for N → ∞ for δ 0 < 1 (see Eq. (A4)).
Lemma 5 (R converges to Q at short times).For any τ < ln N, R(τ) converges to the exact solution Q(τ).
Proof.Using Eq. (A25) and Definitions 2 and 4, we have 1) .(A26) Taking a specific time τ = τ ln N, we can use Lemma 1 to restrict the sum over n to the range n > N − N τ.Let s = N − n.Then (A27) where in the last line we used that the sum over e −τ (1 − e −τ ) s is bounded 1.For shorter times, we can arbitrarily restrict the sum to s < N 0.1 and the bound still works.Thus, ∆ 2 → 0 for all τ < (1 − ε) ln N for any constant ε > 0. we split the sum into two parts, one up to N − N µ and the other from N − N µ to N (for some 0 < µ < 1).For the first part, we have where in the first line we introduced α n = e −τ N(N/n − 1) and to go to the second line we first used that the maxmium with α of αe −α is at α = 1, but that α ≤ N µ−δ 0 to replace α = N µ−δ 0 , where we also used that e −τ ≤ N −δ 0 .Equation (A30) vanishes for µ > δ 0 .For the second part, we have which can be made to vanish, too, by taking µ < 3δ 0 /2 (consistent with Eq. (A4)).
Lemma 7 ( P is asymptotically correct.).For times τ 0 < τ < τ 1 , P converges to P * ∆ 4 (τ) = ||P * (τ) − P(τ)|| → 0. (A32) Proof.We establish ∆ 4 → 0 by bounding the difference between the continuum solution (Definition 3) and the discrete exact solution (Definition 1).The exact discrete solution fulfils Eq. ( 2), whereas the continuum solution obeys We replace the differential by a first-order finite difference with error U By the remainder theorem applied to the variable n, the residue is bounded by where We now bound the error E = P * − P. It is governed by (A36) Note that in the last line we defined with E in 0 (τ) = U 0 (τ) = γ 1 P1 (τ).This yields a bound on We bound the terms in E in individually.The contribution due to the Taylor residue is max Considering first the U 0 = P1 (τ) term, we have so we can take n ≥ 1 in the following.
To bound H 2 , we first note Taking the maximum over n and τ lets us bound the corresponding contribution to E in in terms of G (see Definition 5) Using Lemma 3, we find that these terms all vanish as long as δ 1 < 1/2 and δ 0 > 2/3 (consistent with Eq. (A4)).
We also have to consider the second part of E in that stems from the error introduced in the rate equation.In terms of G, we have ) → 0 (A43) again using Lemma 3. Thus, all terms in E in vanish, which implies that || E|| → 0 as N → ∞.
(A46) Using (A45) and the definition, we have since exp(Γτ) is a stochastic matrix, and thus it does not change the norm.
To bound the integral, we first evaluate the argument where First, we notice that the term in square brackets in Eq. ( A48) is always positive, as we can bound it from below by (A50) Second, we find Thus we can drop the absolute value signs from Eq. (A48) (since the difference in square brackets is positive) while incurring an error of at most 2K.The sum over n then reduces to the boundary term N(N − 1)R > N (τ).Thus, we have The integral vanishes as N → ∞ for any τ > τ 1 .
Finally, we prove the corollary that states that the formulae for magnetization and radiance are correct.
Corollary 1 (Magnetization and radiance (restated)).The predicted magnetization and radiance (9a) are asymptotically correct for any sequence of times Proof.At early times (τ < τ 0 ), we can evaluate the magnetization explicitly using the distribution Q (Lemma 4) Thus, the error in the magnetization is O(1/N, e 2τ /N 2 ).Clearly, the same is true for the radiance.
For intermediate times, we bound the difference between the integral and the sum, Using Taylor's theorem to bound the residue of the zeroth order expansion, we obtain This vanishes by Lemma 3. The same steps for the radiance yields the same error bound, but with an additional contribution n (N/n)e −τ Pn , which also vanishes by Lemma 3. Finally, to cover late times, we prove that at time τ 1 , both the exact magnetization and µ (9a) vanish as N → ∞.Specifically, we show that the probability mass of the solution for n > N µ for µ < 1 − δ 1 vanishes, even after multiplying it by Theorem 1 (Superradiant decay from all-inverted state (restated)).For any sequence of times {τ n > 0 : n ∈ N}, Proof.Together with the triangle inequality, Lemmas 4 and 5 establish Eq. ( 14) for times τ N < ln N. Similarly, Lemmas 6 and 7 together show that Eq. ( 14) holds between τ 0 and τ 1 , and Lemma 8 proves that Eq. ( 14) is true for τ N ≥ τ 1 .
Lemma 9 (Vanishing probability mass 2).In the limit N → ∞, the probability mass of the distribution F(τ) vanishes faster than any polynomial of N in the interval from r = N ε to N for any ε > 0 and τ Proof.Using Eq. ( 21), we have lim At early times (before τ 0 ) we linearize the rates around s = N − 2r − n ≈ 0, which yields Definition 9 (Early solution with loss Q).We define Q n,r by (B5) with Q n and F r defined in Definitions 2 and 8. Q n,r obeys Eq. (5) with rates in Eq. (B4) and can thus equivalently be written as where Γ 0 is the same matrix as Γ in Definition 6, but with linearized rates Eq.(B4).
We solve for the dynamics after the initial phase by moving to a continuum limit in x = n/N, while retaining r as a discrete variable.Specifically, we solve Lemma 10 (Continuum solution with loss).Equation (B7) with initial data at time τ = τ 0 being Pn+2r (τ 0 )F r (τ 0 ) is solved by where Proof.Proof is by substitution.
In Theorem 2, we take τ 0 → 0 in Eq. (B9) to obtain Eq. (22).The main reason to do so is for simplicity.As we prove below, it does not affect the solution.Intuitively, the continuum solution Eq. (B8) can be understood by going back to Eq. (B7) and realizing that the left-hand side is the same as Eq. ( 8).We can think of it as the homogeneous part of the differential equation, which is solved by P. The righthand side only acts on the distribution in r and therefore does not change the distribution in n.Since the right-hand side increases r at a rate depending on γn 2 /N 2 , the distribution in r is a Poisson distribution (Definition 8) but with a time that depends on the dynamics of n, which can in turn be related to the characteristics of the left-hand side of Eq. (B7).
At late times, we linearize around n = 0 to obtain the rates To extend the solution to all times, we need to keep small contributions of order r/N in Γ (4) , as they are responsible for the decay of dark excitations at late times.The rates predict that the remaining collective and dark excitations decay independently, but their rates have different orders in N. In the large-N limit, we thus expect first a rapid decay to n = 0, and subsequently a slow decay of the dark excitations, which can only decay incoherently, as they are dark with respect to the collective decay.
Definition 10 (Late solution with loss).We define the late distribution

Main proof
The proof of Theorem 2 has the same structure as the one of Theorem 1.One small difference is that for times shorter than τ 1 , we bound the difference to the approximate solution P instead of the exact solution P * , as it is simpler.Lemma 15 shows that this is justified, because the difference between the approximate and exact solution also vanishes for times that grow only logarithmically.
Lemma 11 (Q is correct).For any sequence of times {τ N < ln N : N ∈ N}, Q(τ N ) converges to the approximate solution P(τ).Proof.As in Lemma 4, we define the difference between the matrices that generate the evolution Γ 1 = Γ 0 − Γ, with Γ 0 defined in Definition 9 and Γ defined in Definition 7.This allows us to bound ∆ 1 as Eq.(A18) where We split the bound into a part independent of γ and one part proportional to γ, which we will treat separately.The independent part is Due to Lemma 9, we can neglect the contribution from r, such that Eq. (B15) reduces to Eq. (A19) and thus vanishes.The part proportional to γ reads (B16) Note that the rates inside the absolute value sign are all equal or smaller than 1.Using Lemma 1, we can restrict the sum over to values between N − 2r and N − 2r − N µ , which allows us to bound N − n by N µ (recall that r is at most ln N by Lemma 9).Since Q is normalized, we find The difference in the first line vanishes by Lemmas 6 and 9.
The second line is only nonzero for τ < τ 0 = δ 0 ln N, where δ 0 < 1 as per Eq.(A4).For those times, we can restrict the sum over n to values between N−N µ and N (Lemma 1).Defining s = N − n and using s < N µ , we have This bound does not work for a sequence of times {τ n > 0 : n ∈ N}, in which the τ N decay as N µ−1 or faster.However, in this case, we have F r≥1 (τ N ) < N r(µ−1) , which means we can replace F r by δ r,0 , in which case we can use Lemma 6.
Proof.We use Lemma 5 to replace Q n by R n , and Lemmas 1 and 9 to restrict the range of the sum.Using the same reasoning as in Eq. (B20) Lemma 14 (P cont is asymptotically correct.).For times τ 0 < τ < τ 1 , P cont [Eq.(20a)] converges to P (Definition 7) Proof.We start by defining the error where P cont solves the continuum equation and has been defined in Lemma 10.We know that the norm ||E|| vanishes at τ 0 as N → ∞, as per Lemmas 12 and 13.To bound it for the entire time until τ 1 , we consider its time evolution ) Equation (B25) again allows us to bound the error by the integral over the one-norm of E in , as in Eq. (A38).To establish this, we need to bound terms of the form First note that due to the prefactor we can neglect the contribution from the sum from n = 0 to n = N 1−ε .For the sum from n = N 1−ε to N we use n ≥ N 1−ε and Lemma 9 to write The first line in Eq. (B29) is the same as for pure superradiance Eq. (A36) and thus vanishes.In the second line, the first term is of the form Eq. (B26), but with an r out front.Since by Lemma 9, r grows only logarithmically with N, this contribution still vanishes.The second term in the second line vanishes also by Lemma 9.The terms proportional to γ are all of the form Eq. (B26) and vanish, too.Proof.From the previous sections, we know that R n,r (τ 1 ) is a good approximation with the error vanishing.As before, we can bound ∆ 5 by where Γ 1 is the difference between the linearized evolution matrix (Definition 10) and the full evolution (Definition 6).
To bound the integral, we first evaluate the argument (B35) First we note that subleading terms can be neglected, as even bounding them individually leads to a vanishing contribution.
(B39) In the regime m < N µ for some µ < δ 1 /2 the exponential makes the negative term vanish.In contrast, if m > N µ , the fraction T m (τ 1 )/T m+1 (τ 1 ) is smaller than 1, such that we can bound it by 1 and use the result Eq. (A50).
Dropping the absolute value signs in Eq. (B36) means that only the boundary terms survive (up to extra terms that vanish even faster) (B40) Integrating this from time τ 1 to ∞ yields an upper bound to the error that goes to zero as N → ∞.
Note that for this to work it is crucial that we include the decay of r in Eq. (B10).If we had not included this term in the linearized rates, the distribution in r would remain stationary rather than decay and lead to an error growing in time.