N Scaling of Large-Sample Collective Decay in Inhomogeneous Ensembles

,


I. INTRODUCTION
The rate Γ at which an individual two-level system (TLS) emits radiation depends on its coupling to the available vacuum modes and can thus be altered when, e.g., placed inside a cavity [1].For an ensemble of N identical TLSs that are dipole-dipole coupled to each other via the emitted radiation, the decay rate can change significantly even without a cavity, as first noted by Dicke [2].Here, coherence can potentially build up between the emitters, leading to their spontaneous synchronization [3], and an enhanced collective decay rate Γ N > Γ.This is termed superradiance (SR) [4,5], if there is some initial coherence present, and superfluorescence (SF) [6] if the ensemble is initially completely inverted and coherence is initiated by vacuum fluctuations.Since its first observation [7], studies of enhanced collective decay were reported for various experimental platforms such as free-space [8][9][10][11] and waveguide-coupled atomic ensembles [12][13][14][15], solid-state materials [16], superconducting qubits [17] and biological systems [18].Besides being interesting due to its fundamental nature, applications of collective decay can be found, e.g., in quantum optics [19], waveguide QED [20], astrophysics [21], atomic clocks [22,23], and chemistry [18,24].SR/SF can be observed in small, dense samples [25] as well as in large, dilute samples, where propagation effects have to be considered.Both cases exhibit several N scaling signatures [4].In the small-sample case, when N emitters are confined to a volume of dimensions smaller than defined by the emitted wavelength λ, the initial decay rate is enhanced and given by Γ N = N Γ.In the large-sample case, where the N emitters are distributed in a volume > λ 3 , excitation of a homogeneously delocalized coherence is crucial, taking into account a positiondependent phase factor.In the regime of linear-optics SR this is reflected in the timed-Dicke state [26].Due to the extended ensemble volume, a geometric factor µ has to be introduced, describing the fraction of spontaneous radiation that can coherently couple to the ensemble of N emitters in the detected mode [4].The initial decay rate is then Γ N = µN Γ, which is a factor µ smaller than in the small-sample case.µ and µN (≡ N µ ) are the single-atom and collective cooperativities [27], respectively, and N µ becomes the new scaling parameter (see, e.g., [28,29]).Besides enabling collectively-enhanced decay rates, also bursts can be generated in SR/SF if the ensemble is initially inverted, with a delay t D of the burst after the inversion.The photon rate of the burst evolves nonlinearly, reaching its peak at t D with a value ∝ N 2 .This burst occurs when t D , Γ −1 N ≪ γ −1 , where γ is the decoherence rate of the TLSs, which defines a threshold condition.We note that recently a general condition for the occurrence of bursts in waveguide QED systems was derived [30].
As homogeneous excitation conditions in large samples can be challenging to realize (see e.g.[16,31]), the effect of inhomogeneities on collective scattering are still investigated (see, e.g., [32][33][34])-specifically in opticallydense ensembles [11,28,[35][36][37][38].In such inhomogeneous conditions, collective decay occurs only in sub-systems of the ensemble for which the threshold condition is surpassed, while each of them radiates independently [33].This leads to a breakdown of the textbook N scaling [11,28,[39][40][41], potentially hiding the signature of collective decay.As noted by Arecchi and Courtens [42], however, a maximum cooperation number (MCN) N mc ≤ N µ can be determined when propagation losses, delays and inhomogeneous broadening become relevant, which replaces N µ in the equations characterizing collective emission.This concept of a MCN allows for understanding the physical limits to collective decay and has been applied, e.g., to cases where propagation time through the ensemble is longer than the timescale of cooperation build-up [39,40].As we will show in the following, a MCN can also be determined in cases where spectral and spatial inhomogeneities are relevant, thereby recovering the scaling expected for homogeneous ensembles.
In this paper we experimentally study large-sample SF for a broad range of inhomogeneous excitation conditions.For N ≤ 2.2 × 10 5 atoms confined within a hollowcore photonic bandgap fiber (HCPBGF), we determine collective cooperativities up to ∼ 300.Due to intrinsic spatial and spectral inhomogeneities, the experimental data does not exhibit a textbook N scaling.However, by accounting for inhomogeneities in a simple model we de-arXiv:2307.11623v2[quant-ph] 21 Dec 2023 termine a MCN which recovers the textbook scaling over a large range of experimental parameters when dispersion is negligible.Our model allows for physical insight into the limits of collective decay in large ensembles, which is relevant to understand and optimize collective effects in, e.g., free-space [8][9][10][11]43] or waveguide-coupled ensembles [13,14,20,[44][45][46], or solids [16].
The paper is organized as follows: In Sec.II we describe the model system used to study SF.Section III presents an overview of the experimental setup followed by the experimental results in Sec.IV.In Sec.V we then describe the model used to determine a MCN and compare its results to the experimental data.Following a discussion of the MCN model in Sec.VI, we conclude with a summary and outlook in Sec.VII.

II. MODEL SYSTEM
We study SF for a disordered, cold ensemble of 87 Rb atoms loaded into a HCPBGF (NKT Photonics, HC-800-02) [47,48].The ensemble can be considered dilute as even at the peak atomic density N ∼ 10 12 cm −3 we have λ 3 N < 1.The single-mode waveguide provides long-range coupling between distant atoms, thereby facilitating the observation of collective effects [13,49] while avoiding the use of a multi-mode theory to analyze the data [14].The atomic ensemble possesses an inhomogeneous radial density profile due to the guiding potential, which is of similar width as the excitation field [see Fig. 1(c)] [48].Moreover, as the ensemble exhibits a large optical depth (OD) [47], attenuation of the excitation field has to be considered.Thus, the system shows large radial and longitudinal inhomogeneities.Radiation trapping is expected to be negligible due to an ensemble aspect ratio of order 10 4 : 1.In contrast to typical realizations of SR/SF, where the excited state has a finite lifetime of Γ −1 , we here study the decay of an effective TLS comprised of two long-lived ground states [see Fig. 1(a,b)] [50,51].An off-resonant excitation field (termed pump) scatters off the ensemble prepared in |1⟩.
Once the pump is on, spontaneous Raman scattering induced by vacuum fluctuations creates a delocalized coherence between states |1⟩ and |2⟩.This initiates the decay, followed by a buildup of radiation bursts in the Stokes field.We note that this quasi-continuous, off-resonant driving of the system is different from the strong, resonant driving in [29].The single-atom scattering rate of such an effective TLS is for large detunings given by where R B (= 0. the AC Stark shift of the |1⟩ ↔ |3⟩ transition frequency due to the strong and detuned pump field. Studying SF in an effective TLS has several advantages over real TLSs.First, Γ R can be tuned by the intensity and detuning ∆ p of the pump field.Second, Γ R can be made much smaller than the excited state decay rate Γ, which sets the width of the dispersive resonance.Thus, even for large decay rates Γ N ≫ Γ R , we still can have Γ N < Γ, and dispersive effects as discussed in [38] are more likely to be negligible.This leads to a significant simplification when modeling the MCN.Third, the modulation bandwidth of the excitation field has to be larger than Γ N .This can be technically challenging if the single-atom decay rate Γ is already quite large [28].A detailed description including the procedure to load atoms into the HCPBGF can be found in [48].Therefore, we here only give a brief summary.First, we load a magneto-optical trap with around 10 7 87 Rb atoms.Then, we shift the atom cloud towards the HCPBGF tip while compressing it radially using magnetic fields.Up to 2.5% of the atoms are then loaded into a reddetuned Gaussian-shaped far-off-resonant trap (FORT) emerging from the HCPBGF with a numerical aperture NA = 0.092 (6).The FORT guides the atoms radially inside the fiber, while they are basically free-falling in the longitudinal (vertical) direction, resulting in a radial 1/e half-width σ a ∼ 1.7 µm and a length L ∼ 3 cm of the ensemble.We prepare the atoms in state |1⟩ by optical pumping without specific population of Zeeman levels.The linear-polarized pump beam is launched into the fiber in such a way that its polarization is maintained at the exit to allow for efficient polarization filtering [53].Its temporal intensity is controlled by an acousto-optic modulator with a rise time of τ p = 130(5) ns.The radial pump intensity has a near-Gaussian profile with a 1/e 2 half-width of σ p ∼ 2.75 (20) µm.The shot-to-shot intensity fluctuations are < ∼ 5%.We detect the transmitted light using an avalanche photodiode (APD) and a digital oscilloscope.

IV. EXPERIMENTAL STUDIES
In the following we study SF in the effective TLS shown in Fig. 1.Before each measurement, we switch off the trapping potential to avoid inhomogeneous broadening by the FORT.When the pump is resonant with |1⟩ ↔ |3⟩ we can measure N by time-resolved optical pumping [48].The measured shot-to-shot fluctuations of N are < ∼ 3.5 %.We then detune the pump field from |1⟩ ↔ |3⟩ by a variable amount ∆ p .As confirmed by measurements and a numerical simulation, the polarization of the Stokes field is predominantly linear and orthogonal to the pump field.We thus observe the temporally-resolved Stokes radiation by detecting light orthogonally polarized to the pump using a polarizing beam splitter (PBS).We show in Fig. 2(a) exemplary data for the transmitted pump power through the HCPBGF without atoms (blue) and for the transmitted Stokes power with N = 83 × 10 3 atoms loaded into the HCPBGF (orange) versus time.No Stokes light is detected for several hundred nanoseconds after the pump is switched on.Then, a short burst of peak power P s is emitted, followed by another burst of smaller peak power due to a coherent ringing [4,54].This trace, showing a large first burst followed sometimes by a smaller one, represents the dominant temporal shape of the emitted Stokes light for N < 220 × 10 3 and agrees with earlier observations of SF [4,28,54].We note that the transmitted pump power shows complementary dynamics to the Stokes field, i.e., fast absorption, which was termed 'superabsorption' in [55].For our maximum attainable N = 220 × 10 3 , however, the dominant shape changes [green trace in Fig. 2(a)].Here, the ringing is more pronounced and the second burst is largest which is usually not the case [56].As to obtain good statistics for the Stokes signal, we repeat the measurement for each set of experimental parameters (N, ∆ p ) 100 times (see video [57]).The Stokes signal exhibits huge shot-to-shot fluctuations in amplitude, delay, intensity and shape.The burst dynamics occur on a timescale much shorter than expected for the calculated single-atom scattering rate of ⟨Γ R ⟩ r = 2π×45 kHz, where ⟨...⟩ r denotes a radial average over both the pump intensity and atomic density distribution (see Appendix C 2).From each single dataset we extract the delay t D and peak power P s of the first emitted Stokes burst by an automated algorithm to obtain their mean values ⟨t D ⟩ and ⟨P s ⟩ and the corresponding standard deviations.
Figure 2(b) shows exemplary data for the measured mean delays vs. reduced atom number N µ (symbols).Here, Ω p and ∆ p are the same as in Fig. 2(a) and we expect the system to be well-described by an effective TLS.The standard deviations from the mean delay are a significant ∼ 38%, which is much larger than both the ∼ 4% uncertainty of determining t D in a single shot and the expected uncertainty of ∼ 4% due to fluctuating atom number and pump intensity (see Appendix A).As the relative delay fluctuations from the mean ⟨t D ⟩ due to vacuum fluctuations can be estimated as 2.6/ ln N [4] (∼ 0.2..0.3), we attribute the dominant contribution to the observed fluctuations to vacuum fluctuations initiating the decay process [58,59].We note that the standard error of the mean is a factor 10 smaller than the standard deviation shown by the error bars.In the same plot we show the expected dependence (dashed green line) for a homogeneous ensemble of Fresnel number F ∼ 1 (see Appendix F) [60] ⟨t The measured mean delays are larger than Eq. ( 2) predicts, as expected for inhomogeneously-broadened systems [61].
To extract the initial collective decay rate Γ N from our SF data, we cannot resort to determining the decay rate as, e.g., in linear-optics SR due to the nonexponential dynamics.An order-of-magnitude estimation can be done by measuring the temporal widths τ b of the first SF bursts, as Γ N ∝ τ −1 b is a reasonable assumption [4].However, as there is, to the best of our knowledge, no precise theoretical prediction for Γ N (τ b ) in our considered case of SF bursts with ringing, we proceed as follows.From Eq. ( 2) we determine the collective decay rate as the measured ⟨t D ⟩ and N .Fig. 3(a) depicts the measured relative collective decay rate Γ N /⟨Γ R ⟩ r (symbols) vs. N µ for a large range of detunings, along with the expected theoretical dependence Γ N /⟨Γ R ⟩ r = N µ (blue line), where µ ≈ NA 2 /4 = 2.1 × 10 −3 .The collective decay rate is significantly enhanced by more than two orders of magnitude compared to the single-atom case.Only at the largest detunings it scales approximately linearly with N µ , but with a detuning-dependent slope.Moreover, the experimental data does not collapse onto Γ N /⟨Γ R ⟩ r = N µ expected for homogeneous conditions, except for the lowest N and largest ∆ p .Close to resonance, Γ N /⟨Γ R ⟩ r is approximately independent of N .Thus, N µ is not a good scaling parameter.We note that we obtain quite similar values for Γ N /⟨Γ R ⟩ r (differing by a factor of ∼ 1.7) when determining Γ N from the burst widths τ b instead of the mean delays ⟨t D ⟩ (see Appendix B), indicating that a determination of Γ N from ⟨t D ⟩ is indeed correct.Next, we plot in Fig. 3(c) the measured mean peak Stokes power (symbols) as a function of N 2 µ .
We observe a N 2 µ scaling (dashed orange line) only for small N .As N becomes larger, the scaling becomes approximately linear (dashed green line).This change in scaling indicates a significant reduction in cooperativity for increasing N .A similar behavior is seen for other detunings with a breakdown of N 2 µ scaling occurring for smaller N as ∆ p becomes smaller.Our results clearly demonstrate strong collective decay for our dilute, disordered atomic ensemble coupled to an HCPBGF.However, the in part significant deviations from theory for homogeneous conditions evoke a deeper investigation as to understand the limits to collective emission in our system and to optimize it.

V. DETERMINING THE MCN
Based on the concept of MCN [42], we determine in the following the fraction of atoms N mc = ηµN which can radiate collectively, i.e., can be viewed as homogeneously driven.We identify two major mechanisms affecting the MCN.(i) Inhomogeneous spectral broadening, and (ii) pump attenuation.Due to the radial as well as longitudinal dependence of the pump intensity determining decay rate, AC Stark shift, and Stokes gain, we factorize the problem: we first determine radially-averaged quantities and then apply them to average longitudinally (see the Appendix D for further details).
First, we consider inhomogeneous spectral broadening.As noted in [42], the number of cooperative scatterers is reduced by the ratio σ hom /σ inh of (homogeneous) excitation bandwidth σ hom ∼ τ −1 p and inhomogeneous linewidth σ inh if σ hom < σ inh .We estimate σ inh = δS + δΓ R , where δS and δΓ R are standard deviations of the AC Stark shift S and the effective scattering rate Γ R from their average values ⟨S⟩ r , ⟨Γ R ⟩ r .This yields a first correction factor to the MCN.Next, we consider pump attenuation along the propagation direction z through the optically-dense ensemble.Such attenuation is detrimental to excitation of a timed-Dicke state [35,37], and in addition directly affects the decay rate of the effective TLS.We first take (off-)resonant absorption according to the Beer-Lambert law into account by approximating the absorption coefficient by its initial value where α 0 is the resonant OD, and we performed a radial average.Noting further that the Stokes intensity eventually grows as where z ′ = z/L and the Stokes gain G s (∆ p , 0) ≈ ⟨2α 0 Γ R (∆ p , 0)/γ⟩ r , with γ being the ground state decoherence rate, we next consider pump attenuation due to conversion into the Stokes field.We thus make the ansatz p exp[−α(∆ p )z ′ ] for the pump intensity, where the total attenuation factor is In a first approximation we assumed steady-state conditions for both absorption and conversion and that α(∆ p ) is given by its initial radially-averaged value at z ′ = 0.This assumption of steady-state conditions is a good approximation for pump absorption as the duration of transients is roughly given by Γ −1 ≪ t D .However, as the transient regime during buildup of the Stokes field is roughly determined by γ −1 ∼ t D , this ansatz is probably less accurate for pump conversion.We thus introduced a scaling factor β(∆ p ) for the Stokes gain as the only free parameter in our model (see Appendix E for a more detailed discussion).This parameter accounts for neglecting transient conditions during Stokes field buildup exhibiting a time-dependent Stokes intensity evolution [50].
We next calculate the effective collective decay rate as by averaging along the r− and z−coordinates and using Γ R [Ω p (r, z ′ )] from Eq. ( 1) with Following [62,63] we can now determine a shadow factor which incorporates pump attenuation.Here, ⟨Γ R ⟩ r is the radially-averaged single-atom scattering rate at z ′ = 0 used as a reference without attenuation.Combining inhomogeneous broadening and pump attenuation, we finally obtain the MCN as To test our simple model, we determine ⟨t D (N mc )⟩ using N mc instead of N µ in Eq. ( 2) by manually adjusting β as the only free parameter to yield the best agreement with the experimental data [see orange line in Fig. 2(b)].In the range of detunings ∆ p > ∼ Ω p where our effective TLS model should be valid, we determine β(∆ p ) = 0.182 − 6.1 × 10 −3 ∆ p /Γ from a linear fit.Although we currently have no exact explanation for this dependence on ∆ p , we suspect that it is due to a larger temporal ratio of transient to steady-state conditions at larger detunings (see Appendix E for more details).For smaller ∆ p , β increases dramatically, indicating that our model is inapplicable in this regime.We use β(∆ p ) to calculate N mc (∆ p , N ) for all our data and plot the results for the relative MCN N mc (∆ p , N )/N µ in Fig. 4.Only for the smallest N and largest ∆ p the MCN is close to N µ .As N increases and ∆ p decreases, however, the MCN becomes smaller and only ∼ 10% of the usual collective cooperativity N µ contributes to collective scattering.We can now plot Γ N (N mc )/⟨Γ R ⟩ r and P s (N 2 mc ) , respectively, to test their scaling [see Fig. 3(b,d)].In contrast to Fig. 3(a) where N µ was used, Γ N (N mc )/⟨Γ R ⟩ r collapses onto the linear dependence N mc = ηµN over a large range of detunings and atom numbers.Moreover, we also observe a quadratic scaling of ⟨P s ⟩ with the MCN over almost the full range of atom numbers.Only at the largest N mc the scaling of ⟨P s ⟩ changes.In this regime [green trace in Fig. 2(a)], however, the second instead of the first burst exhibits the largest power in the dominant fraction of measurements, which is not the case in standard SR theory [4].Therefore, the breakdown of N 2 scaling is not surprising, and we excluded this data point from the fit.

VI. DISCUSSION
Our results demonstrate that the concept of MCN can recover the scaling of collective decay in homogeneous ensembles for inhomogeneous conditions.This provides insight into the limits of collective decay in our setup.To this end we plot in Fig. 4 the boundaries for which pump attenuation (black line) [11,35] and inhomogeneous broadening (white line) are becoming detrimental to collective decay.Our system is obviously limited by pump attenuation for large detunings.Only for ∆ p < ∼ 9Γ inhomogeneous spectral broadening is becoming relevant.Besides supporting the picture of linear-optics SR presented in [38,64], our results demonstrate that one can find a regime that renders the radial inhomogeneities insignificant for atomic ensembles in a hollow-core fiber.For appropriately large detunings, almost all atoms contribute to collective scattering, leading to a pronounced enhancement of the collective decay rate, at least for negligible dispersion.
Despite its simplicity, our model seems to capture the relevant physics.Only for detunings ∆ p < ∼ 7Γ, where the approximation of an effective TLS is not valid any-more for Ω (0) p = 6.4Γ, deviations from the experimental data become more significant.Furthermore, an estimation shows that dispersion is likely to become relevant here as frequency components of the SF bursts are lying within the dispersive bandwidth of the excited state [38].This is currently not captured by our model.

VII. CONCLUSION AND OUTLOOK
We presented an experimental study of superfluorescence in a disordered ensemble of cold atoms coupled to a hollow-core fiber.We demonstrated a decay rate collectively-enhanced by more than two orders of magnitude for N < ∼ 2.2 × 10 5 atoms.By implementing an effective TLS we were able to study the decay with negligible dispersion for a large range of detunings and optical depths.Thereby we could tune the amount of inhomogeneity, studying its effect on the SF decay and developing a simple theoretical model considering inhomogeneous broadening and attenuation of the excitation field.Using this model, we determined a maximum cooperation number N mc of atoms that can be treated homogeneously and thus decay collectively.We demonstrated that using N mc instead of µN , the N scaling known to collective decay in homogeneous ensembles could be recovered.This allowed for a physical understanding of the limits to enhanced collective decay in our waveguidecoupled ensemble.Our results provide a simple tool to understand and optimize collective scattering in inhomogeneous extended ensembles, as can be found, e.g., in astrophysics [21], solids [16,31], waveguide QED [20], novel atomic clocks [22,23], and quantum optics [19].
Future studies could try determining a further scaling factor for the MCN incorporating dispersion which could help finding an analytical model of collective decay for dispersive ensembles [65].Also, applying the here demonstrated approach to other types of experimental systems exhibiting inhomogeneities would be interesting as to, e.g., determine their limits to support collective decay.Furthermore, it might be interesting to study the regime where the SF bursts exhibit an unusual dominant temporal shape [see Fig. 2(a)].Finally, our demonstration of collectively-enhanced Stokes scattering rates could be used during the pair-generation process in quantum light sources based on spontaneous four-wave mixing [66] thereby reducing their required pump power.
Comparing the experimental data to the theoretical model requires knowledge of several experimental parameters.The geometric factor describing the fraction of spontaneous light emitted into the fiber core mode is approximately given by µ ≈ NA 2 /4, where NA is the numerical aperture of the HCPBGF.Taking the measured 1/e 2 mode field radius σ p = 2.75(20) µm of the near-Gaussian-shaped intensity distribution, we determine NA = λ/(πσ p ) = 0.092( 6) with λ = 795 nm.This results in µ = 2.1(3) × 10 −3 .The decay rate of the effective TLS depends on the branching ratio R B = 0.5 [52] into the ground state |2⟩, the Rabi frequency Ω p of the pump, and its detuning ∆ p , as does the AC Stark shift S. We measure the pump power (∼5 % shot-to-shot fluctuations) at a certain location inside the experimental setup and correct for transmission losses from the HCP-BGF to this location.Then, using the measured mode field diameter, we calculate the peak Rabi frequency Ω (0) p according to the parameters given in [52].The detuning ∆ p is controlled by the frequency of the acousto-optic modulator (AOM) modulating the pump with a rise time τ p = 130(5) ns.The frequency uncertainty of the AOM is negligible leaving only the uncertainty of the pump laser lock-point.We estimate this uncertainty as < ∼ Γ/10, which can therefore be neglected as well in our analysis.The total atom number N is determined as described in Sec.IV with a shot-to-shot uncertainty of 3.5 %.The radial atomic density profile was investigated in [48] from which we assume a Gaussian-shaped radial density with a 1/e half width of σ a = 1.7 µm.Using an atomic density N (r, z) = N (r) and Ω p (r) we determine the peak optical depth as where σ 13 is the resonant absorption cross-section of the transition |1⟩ ↔ |3⟩ according to [52], r c is the core radius of the HCPBGF and we determined N (r) such that 2π L dr r N (r) = N .As the final parameter to determine, there is a residual decoherence rate between the two ground states due to residual magnetic field gradients and time-of-flight broadening once the guiding potential is off for the measurements.This value was determined from light storage and retrieval experiments as γ 0 ∼ 0.057(6)Γ [66].Determination of the collective decay rate for our experimental data is not as straightforward as, e.g., in the linear-optics SR regime, where a fit to the exponentially decreasing detected power yields the decay rate [11].In our case of SF with strong bursts and ringing we thus have to apply another method.As discussed in Sec.IV, we determined Γ N from the mean delays ⟨t D ⟩ of the first SF bursts for which a theoretical prediction exists [see Eq. ( 2)].As this method is somewhat indirect, we apply in the following an alternative method.As was done, e.g., by Okaba et al. in the case of SR with bursts and ringing [14], the width of the first burst can be used as a measure for the initial collective decay rate.In the case of small-sample SR with a single burst the relation between the temporal burst width τ b (FWHM) and the initial collective decay rate is given by Γ N = 3.5/τ b [4].For largesample SR/SF, corresponding to our experiment, we are, however, not aware of an analogous derivation for Γ N (τ b ).Despite the lack of of a precise theoretical prediction, we determined Γ N from our data using the first burst widths and the small-sample relation Γ N = 3.5/τ b .The results for Γ N /⟨Γ R ⟩ r vs. N mc are plotted in Fig. 5. Comparing this plot to Fig. 3(b), we note the following: (a) The order of magnitude of Γ N is the same.(b) Plotting Γ N vs. the MCN recovers a near-linear scaling in both cases.(c) Scaling Γ N in Fig. 5 by a factor of ∼ 1.7 results in a scattering of the datapoints around Γ N /⟨Γ R ⟩ r = N mc (blue line).(d) The scattering of the data is larger when using the burst widths compared to the mean delays.
In view of these results we conclude that the initial collective decay rate of SF can be determined either from the mean burst delays or from their temporal widths.The former method produces more reliable results and provides a direct comparison to theory.Using rather the burst delays instead of their widths is also suggested by the uncertainties in the fitting process.Our data shows that the delays can be determined with a smaller relative uncertainty than the widths.This is intuitive as the width is also dependent on the fitted amplitude.Finally we note that the widths determined from the fits to the single-shot data are in general larger than the raw data would suggest.This might explain why the datapoints in Fig. 5 are below the theoretical line.Another reason for this discrepancy could, however, be that the relation Γ N = 3.5/τ b , valid for small samples, cannot be directly applied to the large-sample case.

Appendix C: Determination of the Maximum Cooperation Number
In the following we present the model used to calculate the maximum cooperation number N mc .It is mainly based on the ideas presented by Arecchi and Courtens [42], Bienaimé et al. [35], and Bachelard et al. [62].First, inhomogeneous broadening reduces the number of cooperative atoms [42].Second, attenuation of the pump field as it propagates through the optically-dense ensemble reduces the effective scattering/decay rate Γ N , which can be viewed as a reduction of the number of synchronized scatterers.We account for this by calculating a factor η s from an effective scattering rate in the spirit of the shadow effect introduced by Bachelard et al. [62].We note that in our case the decay rate Γ R itself depends on the (attenuated) excitation field, which is in contrast to decay in a true TLS where the decay rate is constant and only the excitation is inhomogeneous.
To find a simple analytical model, we make several assumptions, as the actual system exhibits complicated nonlinear spatial and temporal dynamics.First, there is a radial variation of the atomic density distribution N (r) = N 0 exp −r 2 /σ 2 a , and we assume a constant density along the propagation direction z.Second, the pump intensity I p (r, z) ∝ Ω 2 p (r, z), determining both Raman scattering rate and AC Stark shift of the ground state, is potentially attenuated along the propagation direction z and also exhibits a radial variation Ω 2 p (r, 0) = Ω 2 p (0, 0) exp −2r 2 /σ 2 p .Its 1/e 2 half width is only slightly larger than σ a .We use the fact that Ω 2 p (r, z) can be factorized as Furthermore, inhomogeneous broadening due to the radial variations is expected to dominate over longitudinal variations as it is also present for large detunings.These assumptions allow for considering the radial and longitudinal variations independently.In the following we start by calculating the radial averages weighed by the atomic density distribution at the start of the medium at z = 0 for obtaining the inhomogeneous linewidth.Then, we take the radial and longitudinal average to determine the effective scattering rate.Using these results we then determine a MCN.

Inhomogeneous Broadening
As discussed by Arecchi and Courtens, inhomogeneous broadening reduces the number of cooperative atoms by the ratio of (homogeneous) excitation bandwidth σ hom to inhomogeneous linewidth σ inh [42].The excitation bandwidth can be estimated by the rise time τ p of the pump field, σ hom = 1/τ p .The inhomogeneous linewidth of the effective TLS can be estimated as follows.The pump intensity inside the HCPBGF determines the AC Stark shift S of the ground state F =1 and the decay of the effective TLS We therefore calculate the radial averages of Ω 2 p (r, 0), S(r, 0), and Γ R (r, 0) at the beginning z = 0 of the medium, weighed by the atomic density distribution as where f (r) = Ω 2 p (r, 0), S(r, 0), Γ R (r, 0) .Using the radially-averaged values, we now estimate the inhomogeneous linewidth σ inh as follows.Due to the radial Gaussian dependence of the strong pump field, atoms in the ensemble of radial width σ a will experience an inhomogeneous AC Stark shift S(r) of the ground state F =1 as well as an inhomogeneous decay rate Γ R (r).As a measure for the corresponding linewidth we calculate the standard deviation for f (r) = [S(r, 0), Γ R (r, 0)] and set Using the measured rise time τ p of the pump pulse, we obtain the first of two scaling factors, to calculate the MCN.Note that we only apply this scaling factor if the inhomogeneous bandwidth is larger than the homogeneous bandwidth.Considering only inhomogeneous broadening in the MCN, we obtain the results shown in Fig. 6(a) which have to be compared to those shown in Fig. 3(a).Obviously, taking inhomogeneous broadening into account slightly improves the scaling but is only relevant in the range of ∆ p < ∼ 9Γ.

Effective Scattering Rate
Having determined a scaling factor which accounts for inhomogeneous broadening (in the radial direction), we next turn to the variations in the longitudinal direction z.This is relevant for our considered system, as propagation through the optically-dense atomic ensemble can attenuate the pump field.To this end, we calculate an effective collective scattering rate weighed by the radial atomic density distribution N (r) where and Ω (0) p [= Ω p (0, 0)] is the peak Rabi frequency of the pump.We assume attenuation of the pump intensity according to Beer-Lambert's law along the dimensionless distance z ′ = z/L with an attenuation factor α(∆ p ) affecting also S(r, z ′ ) in Eq. (C2).We obtain the corresponding collective scattering rate ⟨Γ N ⟩ r without pump attenuation by setting α(∆ p ) ≡ 0 in Eq. (C8).Similar to Bachelard et al., where a shadow factor was defined to determine an effective reduction in the total scattering cross-section of an optically-dense ensemble [62], we now define a shadow factor where ⟨Γ R ⟩ r is the radially-averaged single-atom scattering rate at z ′ = 0 used as reference without attenuation.η s is thus a measure for the reduction of the collective scattering rate due to pump attenuation.As Γ eff N r,z ∝ N we identify η s µN as an effective number of atoms that will scatter collectively when considering attenuation of the pump field.

Maximum Cooperation Number
Finally, combining inhomogeneous broadening, longitudinal pump attenuation, and the geometric factor µ we determine the MCN as Appendix D: Modeling Pump Attenuation In order to apply the general model described above to our experimental data, we need to identify which factors cause attenuation of the pump, i.e., we need to model the (dimensionless) attenuation factor α(∆ p ).The first obvious consideration is (off-)resonant absorption.Our experiment allows for loading N < ∼ 2.2 × 10 5 atoms into the HCF.With an OD per atom of α * = 2.75 × 10 −3 on the transition 2 S 1/2 F =1 ↔ 2 P 1/2 F ′ =2 we therefore achieve resonant optical depths of α 0 < ∼ 600.The detuning-dependent optical depth at z = 0 for large detunings is approximately given by where we have taken a radial average (see Appendix C 1).For our experimental parameters we reach α(∆ p ) > ∼ 1 for ∆ p < ∼ 12Γ.Setting α(∆ p ) ≡ α(∆ p ) in Eq. (C8), where we assumed in a first approximation that the attenuation factor does not depend on z and is given by its initial value at z = 0, yields a MCN N abs mc = η inh • η abs s • µ • N for each data point in Fig. 3(a).Plotting the initial relative collective decay rate Γ N /⟨Γ R ⟩ r vs. N abs mc we obtain the plot shown in Fig. 6(b).Note that so far there is no free fitting parameter.Compared to Fig. 6(a), where only inhomogeneous broadening was considered, N abs mc is a further improved scaling parameter compared to the collective cooperativity N µ .Yet the data points do not collapse onto the blue line showing Γ N /⟨Γ R ⟩ r = N abs mc .We therefore consider additional attenuation of the pump.To this end we note that the Stokes buildup happens on a timescale determined by the Stokes gain [50] and that the derivation of Eq. ( 2) was done in a regime of exponential growth of the SF burst [60].We thus make the ansatz where the total attenuation factor is given by α(∆ p ) = α(∆ p ) + β(∆ p ) • G s (∆ p , 0). (D3) We thus assume that the pump intensity I p decreases exponentially due to absorption and conversion into the exponentially-growing Stokes field I s (z ′ ) ∝ exp [β(∆ p )G s (∆ p , 0)z ′ ] with a scaling factor β(∆ p ).The initial radially-averaged Stokes gain can be approximately written as where γ is the ground state decoherence rate.For our all our experimental parameters we determine G s (∆ p , z ′ =0) > 1.This ansatz corresponds to steadystate conditions in [50] while neglecting transient conditions during Stokes field buildup, as during the latter one the Stokes intensity evolution is additionally timedependent which complicates the theoretical description [50].To account for this simplification we introduced the scaling factor β(∆ p ) for the Stokes gain, which is allowed to vary with ∆ p (see Appendix E).We note that this scaling factor accounting for partially transient conditions during the buildup of coherence is only applied to the Stokes gain but not for absorption.This is motivated as follows: Although absorption also exhibits a transient after switch-on of the pump, the timescale is determined by the excited state decay rate Γ, corresponding to ∼ 30 ns which is much smaller than the timescale of SF emission.On the other hand, the transient timescale for generation of the Stokes field is determined by the ground state coherence time γ −1 [50], which is comparable to the timescale of SF emission.Thus, transient conditions are likely more relevant in the gain than in the absorption term.
In order to calculate the initial Stokes gain, we have to determine the decoherence rate γ.We identify three major contributions to γ for our effective TLS.First, there is a residual decoherence rate between the two ground states determined by residual magnetic field gradients and time-of-flight broadening.This value was determined from light storage and retrieval experiments as γ 0 ∼ 0.057(6)Γ [66].Second, the inhomogeneous AC Stark shift of the the ground state by the pump field causes decoherence.We approximate this rate by the previously calculated standard deviation δS (see Appendix C 1).Third, the Raman scattering rate Γ R leads to decoherence due to decay between the two ground states.Therefore, we set the total decoherence rate as γ = γ 0 + δS + ⟨Γ R ⟩ r . (D5) Combining these results, we determine the MCN N mc for all measurement parameters (∆ p , N ).The results are shown in Fig. 4. Complementary to Fig. 6(a), where we only considered inhomogeneous broadening, we plot in Fig. 6(c) the results when considering only pump attenuation while neglecting inhomogeneous broadening.As we can see, attenuation of the pump field is relevant across the whole range of detunings, however, it fails to provide the correct scaling for ∆ p < ∼ 9Γ.Only a combination of pump attenuation and inhomogeneous broadening yields the expected scaling.

Figure 1 (
Figure 1(c) shows our simplified experimental setup.A detailed description including the procedure to load atoms into the HCPBGF can be found in[48].Therefore, we here only give a brief summary.First, we load a magneto-optical trap with around 10 7 87 Rb atoms.Then, we shift the atom cloud towards the HCPBGF tip while compressing it radially using magnetic fields.Up to 2.5% of the atoms are then loaded into a reddetuned Gaussian-shaped far-off-resonant trap (FORT) emerging from the HCPBGF with a numerical aperture NA = 0.092(6).The FORT guides the atoms radially inside the fiber, while they are basically free-falling in the longitudinal (vertical) direction, resulting in a radial 1/e half-width σ a ∼ 1.7 µm and a length L ∼ 3 cm of the ensemble.We prepare the atoms in state |1⟩ by optical pumping without specific population of Zeeman levels.The linear-polarized pump beam is launched into the fiber in such a way that its polarization is maintained at the exit to allow for efficient polarization filtering[53].Its temporal intensity is controlled by an acousto-optic

FIG. 2 .
FIG. 2. (a) Exemplary traces of the transmitted pump power without atoms (blue, 16 averages, scaled down by factor 0.5) and the Stokes power with atoms loaded into the HCPBGF (single-shot data).The peak pump Rabi frequency is Ω (0) p = 6.4Γ, ∆p = 18.4Γ with N = 83 × 10 3 (orange) and N = 220 × 10 3 (green).t = 0 corresponds to the moment the pump power reaches 50% of its mean maximum value.The extracted delay tD of the orange signal is shown by the grey shaded area.(b) Measured mean delays of the first emitted SF burst vs. Nµ for ∆p = 18.4Γ, with error bars representing the standard deviation.The expected theoretical dependence [Eq.(2)] is shown considering Nµ (dashed green line) as well as the MCN Nmc (solid orange line) for β = 0.07.

FIG. 5 .
FIG.5.Same as Fig.3(b) but using the widths of the first bursts instead of the mean delays to determine ΓN .