Storage and release of subradiant excitations in a dense atomic cloud

We report the observation of subradiance in dense ensembles of cold $^{87}$Rb atoms operating near Dicke's regime of a large number of atoms in a volume with dimensions smaller than the transition wavelength. We validate that the atom number is the only cooperativity parameter governing subradiance. We probe the dynamics in the many-body regime and support the picture that multiply-excited subradiant states are built as a superposition of singly-excited states that decay independently. Moreover, we implement an experimental procedure to release the excitation stored in the long-lived modes in a pulse of light. This technique is a first step towards the realization of tailored light storing based on subradiance.

The interaction between a single two-level atom and radiation is well understood: the atomic response is described by a resonance frequency and a decay rate. When considering more than one emitter in a volume with dimensions smaller than the transition wavelength, this response may be altered, as was proposed in the pioneering work of Dicke [1,2]. Indeed, light-induced interactions modify dramatically the behavior of the ensemble and its response becomes collective. In particular, the decay rate of excitations hosted in the ensemble can be starkly modified. Superradiance, i.e. a decay of the excitation at a rate faster than the single-atom one, has been verified experimentally in atomic systems from ions to dilute clouds of atoms [2][3][4][5][6]. The study of its counterpart, namely subradiance with a decay rate smaller than the atomic one, has been restricted to a handful of works [7]: direct observations were reported in a pair of ions at variable distance [3] and in molecular systems [8][9][10]. Recently, it was also observed in a cold, dilute atomic cloud [11][12][13] and as a line-narrowing in an ordered 2D layer of atoms [14].
Engineering subradiant states has drawn an increasing attention since it might pave the way to several applications. For instance, the possibility to store an excitation in subradiant modes and to address it in real time while the excitation is stored has inspired proposals to use it as a storage medium [15][16][17][18][19]. Secondly, the narrowing of the line associated to subradiant modes and their subsequent enhanced sensitivity to external fields could be a promising application for metrology [15,20,21]. Recent proposals have also suggested to use subradiance as a tool for quantum information processing and quantum optics [22,23].
All these proposals have been formulated in ordered systems with small inter-particle distances,r λ where λ is the wavelength of the atomic transition. Motivated by this, we take a first step in this direction by exploring the regimē r < λ , but in the disordered case using dense clouds of 87 Rb atoms, characterized by a peak density ρ 0 satisfying ρ 0 λ 3 1 (r = ρ −1/3 0 ). Furthermore, the ensembles we produce have a prolate shape with typical radial size of ∼ 0.5 λ and an axial one of ∼ 5 λ . We thus closely approach Dicke's regime where many emitters are trapped in a volume comparable to * igor.ferrier-barbut@institutoptique.fr the wavelength of their transition. This regime introduces several important differences with respect to the case of a dilute extended cloud studied in refs. [11][12][13]. Firstly, here the parameter governing the collective properties is no longer the optical depth, but rather the number of atoms in the cloud. Indeed the cooperativity C is the ratio between the atom number N and the number of modes efficiently coupled to the system M, i.e., C = N/M [24]. In a cloud with a volume much larger than λ 3 , this parameter is the optical depth on resonance, which was experimentally shown to govern collective effects [4,11]. Approaching Dicke's regime, as we do here, the ensemble becomes efficiently coupled only to a single mode and thus the cooperativity should be the atom number N [25]. Secondly, since in our clouds kr ∼ 1 (k = 2π/λ ), all the terms of the dipole-dipole interaction [26] play a role, as opposed to the dilute regime where only a radiative 1/r term is considered [24].
Here we observe subradiance in the time domain in a cloud operating near Dicke's regime. First we validate the characteristic dependence on the atom number. Second we explore the storage of light in long-lived multiply-excited states. Varying the intensity of the excitation laser, we characterize these subradiant states containing few to many excitations. Our finding supports the idea that multiplyexcited subradiant states are built as a superposition of singly-excited states in random ensembles, similarly to what was recently predicted in ordered 1D systems [18,[27][28][29]. Finally, we demonstrate dynamical control of subradiance while excitations are stored by releasing it on demand via the application of a laser. This real-time control of the coupling of an ensemble to the electromagnetic modes while it hosts an excitation offers new possibilities for light storage.

I. EXPERIMENTAL SETUP
A detailed description of the experimental setup can be found elsewhere [31,32]. Briefly, as sketched in Fig. 1(a), it is composed of four aspherical lenses with large numerical aperture (0.5) in a maltese-cross configuration [33]. We use the x high-resolution optical axis to create an optical tweezer at a wavelength λ trap = 940 nm, with a tunable waist (range 1.8-2.5 µm). Exploiting a gray-molasses loading on the D1 arXiv:2012.10222v1 [physics.atom-ph] 18 Dec 2020 line, we trap 5000 87 Rb atoms in the largest tweezer at a temperature of about 650 µK in a 4.2 mK trap. In this configuration, the trapping frequencies are ω r 2π × 81 kHz and ω z 2π ×7 kHz, where ω r and ω z represent the radial and the axial directions. The central density of the cloud in these conditions is ρ 0 /k 3 = 0.3 ± 0.1 (r 0.2 λ ), where k = 2π/λ . The trap can then be compressed either by increasing the power of the trapping beam or reducing its waist [32].
We use the F = 2 → F = 3 transition on the D2 line with wavelength λ 780 nm, linewidth Γ 0 2π × 6 MHz and saturation intensity I sat = 1.6 mW/cm 2 . In order to excite the cloud, we switch off the trap and shine a 150 ns-long pulse of resonant light along the y − z direction of Fig. 1(a). This duration is long enough to reach the steady state during the excitation. After 1 µs, the atoms are re-captured in the tweezer. This time is short enough (< 1/ω r ) that the density remains constant during this trap-free time. We repeat the same sequence up to 20 times on the same cloud (depending on the trap geometry), checking that the atom number is reduced by less than 10% at the end of the pulse sequence. We then repeat this sequence on 3000 to 10000 clouds at a 2 Hz rate to obtain one trace (photon count per 0.5 ns bin versus time). For the largest atom number, N 5000, we collect in the steady state typically 0.01 photon per pulse in a 0.5 ns time bin. The excitation pulse is controlled by means of two acoustooptical modulators in series ensuring a high extinction ratio. The characteristic rise/fall time is 10 ns, the temporal shape of the pulse is shown in Fig. 1(b) (gray line). The 1 mm excitation beam waist is much larger than the cloud size so that all the atoms experience the same Rabi frequency. We can vary the saturation parameters s = I/I sat up to s 250 with I sat the saturation intensity (we have calibrated this intensity independently using dilute clouds). Exploiting the two highresolution optical axes, we collect the fluorescence light into two fiber-coupled avalanche photodiodes (APDs) operating in single photon counting mode, one aligned along the axial direction of the cloud (APD //), the other perpendicularly to it (APD⊥), as sketched in Fig. 1(a).

II. SUBRADIANCE NEAR DICKE'S REGIME
In Fig. 1(b) we report two time-resolved fluorescence traces recorded along x, obtained for two clouds with respectively 300 atoms (ρ 0 /k 3 0.02, purple line) and 5000 atoms (ρ 0 /k 3 0.3, blue line) following the switch-off of the excitation laser for s 27. We also show the solution of the optical Bloch equations (OBEs) for a single atom, solved for the measured pulse shape. In the low atom number case, the decay of the excitation is very well described by the OBEs, indicating that the atoms act as independent atoms. On the contrary for large atom number, the fluorescence decays nonexponentially: we observe first a decay at a rate larger than the single atom decay (superradiance), followed by a slower one (subradiance). Moreover, as shown in Fig. 1(c), we observe that superradiance occurs mainly in the axial direction of the cloud while the emission of subradiant excitation is observed in both directions. The finite switch-off time of the driving pulse limits the superradiant decay that can be observed. For this reason, in the rest of the paper we will focus our attention on the subradiant tail, leaving a detailed study of superradiance for future works.
All the measurements reported here have been performed with resonant and linearly polarized light, in the absence of Zeeman optical pumping, thus exciting a multi-level system. We did not observe that the polarization direction impacts the observed subradiance. We have further observed that subradiance is essentially unchanged within our dynamical range in the presence of a magnetic field and a circularly polarized pulse with prior optical pumping (see appendix E). This suggests that the internal structure does not play a major role for subradiance in our regime. Finally, we have checked that the subradiance is unchanged when we vary the detuning of the excitation laser around the atomic resonance (appendix B). This indicates that radiation trapping of light in the cloud seen as a random walk of the photons before escaping can not explain the observed slow decay [34]. For our small dense cloud, and contrarily to the case of dilute, optically thick clouds [11], this is expected, as the photon mean-free path l sc = 1/(ρσ sc ) (with ρ the atomic density and σ sc = 3λ 2 /2π the resonant cross-section) is smaller than the mean inter-particle distance.
To further support our observations of subradiance, we perform numerical simulations in a simplified setting of two-level atoms. For the large number of atoms involved in the experiment, the ab-initio simulation by use of a master equation is beyond reach and we therefore resort to approximations. We use a non-linear coupled-dipoles model [30,31] consisting of a coupled system of OBEs, given in appendix A. It formally amounts to a mean-field approximation, assuming that the density matrix of the system can be factorized [30,35,36]. It allows one to take into account saturation effects of individual atoms. We numerically solved the equations for N = 200 atoms at a density ρ 0 /k 3 = 0.3 [ Fig. 1(d)]. The results do not feature superradiance, but do yield subradiance. The origin of the superradiance observed in the experiment and not present in the mean-field simulation is left for future investigations. The prediction of subradiance in our simulations could suggest that the mean-field model is enough to account for our observations. However as we will show below, it fails to reproduce our results in the strongly saturated regime, even qualitatively. In the remaining of this section, we characterize the observed subradiance as a function of atom number.
To reveal the characteristic scaling of subradiance in our regime, we investigate the crossover between the low atom number regime, where the system behaves as an ensemble of non-interacting emitters, to the large atom number one. To adjust the atom number N, we release the atoms from the trap and recapture them after a variable time prior to sanding the burst of excitation pulses. This technique allows us to reduce the atom number by a factor more than 10 with negligible heating, thus leaving the cloud sizes unchanged. To analyze the experimental data, we fit the decay with a phenomenological model using either a single exponential decay or the sum of two decays with different characteristic times. The fitting function is decided based on a χ 2 criterion (see appendix C for more details). The decay time is extracted from data averaged over tens of thousands of realizations and thus corresponds to an average subradiant behaviour.
We report in Fig. 2(a) the results of this analysis. We observe that as N grows, a clear subradiant tail appears, and that the characteristic decay time is an increasing function of N. To determine the parameter that governs the cooperativity, we acquire three different data sets shown as different symbols in Fig. 2. They are taken in three different trapping geometries leading to different cloud sizes (appendix D). We plot the same data in Fig. 2(b) as a function of the optical depth along x, b 0 = ρ 0 l x √ 2πσ sc , and in Fig. 2(c) as a function of the peak density ρ 0 = N/[(2π) 3/2 l x l 2 r ] (l x,r are Gaussian sizes). We observe that the data nearly collapse on a single curve when plotted as a function of atom number rather than as a function of optical depth. Furthermore, the decay times cannot be described either by the cloud density only, which is a local quantity. This is expected given the long-range character of the dipole-dipole interaction. We verify that the amplitude of the subradiant decay measured as the relative area in the tail (see appendix F) is also governed by the atom number only. Therefore the two parameters describing subradiance namely lifetime and amplitude are solely governed by N, which is the cooperativity parameter for this regime as explained in the introduction. This behaviour distinguishes our results from the ones in dilute systems where the cooperativity parameter is b 0 [24], and indicates that we are approaching the Dicke limit. The imperfect collapse of the experimental data in Fig. 2(a) might be due to the fact that the system size along x is still larger than the excitation wavelength.

III. STUDY OF MULTIPLY-EXCITED SUBRADIANT STATES
Subradiance corresponding to the presence of long-lived excitations, a natural application would be to store light in an atomic medium. Storing multi-photon states would require long-lived multiply-excited states. Therefore understanding the nature of these excitations is a prerequisite for the application of multiple excitation storage [12]. Here we investigate this question experimentally by varying the intensity of the excitation laser.
Considering first the strong driving limit, the ensemble is prepared in a product state where each atom is in a mixture of the ground |g and excited state |e , with density matrix ρ = 1 2 N (|e e| + |g g|) ⊗N . In general, this mixture comtains Decay time subradiant components, as was discussed in ref. [12] for the case of two atoms. Reaching such a steady state during the excitation therefore leads to the initial excitation of long-lived subradiant states. However it does not preclude the further population of these states during the early decay following the switch-off of the laser excitation [12,28,37]. For large N, it was shown [18,[27][28][29] in the case of ordered 1D arrays, that subradiant states containing n exc > 1 excitations are built from a superposition of subradiant states of the single excitation manifold, which decay independently with their respective lifetime Γ (1) n , as exemplified for n exc = 2 in the caption of Fig. 3. The signal resulting from the decay of an n exc > 1 state is ∝ ∑ n exc n=1 exp(−Γ (1) n t). The interest of this ansatz stems from the fact that the singly-excited states can be calculated via a model of classical coupled dipoles [38]. It is however an open question whether this simple picture holds in the disordered case studied here. As we show below, our experimental findings support the picture of multiply-excited subradiant states constructed as a superposition of independent singlyexcited subradiant states.
To control the number of excitations in the system, we vary the saturation parameter s = I/I sat between s 0.01 and s 30. For these measurements we work without compression of the trap and the atom number is fixed 4500. All the time traces acquired in this section are reported in appendix H. We extract the subradiant lifetime using the procedure introduced earlier. In addition, we calculate the tail fluorescence by summing the signal for times larger than t 0 + 4 Γ 0 , where t 0 marks the switch-off of the excitation pulse. We observe that the subradiant decay time is constant over three orders of magnitude of the excitation intensity [ Fig. 3(a)] and that the tail fluorescence increases with excitation intensity before saturating at an intensity smaller than I sat [ Fig. 3(b)].
The constant lifetime suggests a first simple description of the data as the excitation of a single mode. We thus use a single-mode approximation to describe the subradiant tail population with the following expression assuming a saturation behaviour similar to that of a single atom: ∝ αs/(1 + αs), with α the fit parameter. We report the result of this fit as the solid line in Fig. 3(b). From the extracted α = 3.4 (5) and assuming that the lifetime of a given mode τ dictates its saturation intensity I sat ∝ 1/τ as for a single atom [39], we obtain τ = 3.4/Γ 0 , represented as a solid line in Fig. 3(a). This result agrees remarkably well with the direct measurement of the decay rate. Thus the single mode approximation describes very well our data, seeming to confirm the validity of this approximation. However in the saturated regime the long-lived modes leading to the subradiant decay host up to 10 % of the total excitations (see appendix F Fig. A4), which for a fully-saturated (i.e. 1/2 excitation per atom) cloud of 5000 atoms means several hundreds of excitations. Despite this large number of excitations, the decay rate remains the same demonstrating that, in the subradiant tail, the rate at which excitations decay is independent of the density of excitations in the system. This finding is consistent with multiply-excited states constructed from a large population of singly-excited subradiant states which decay independently. In the opposite case of strong interactions between singly-excited subradiant states, we would have observed an excitation density-dependent decay rate due to additional decay processes induced by interactions between excitations. Experimentally we observe average decay times of about 3/Γ 0 , in agreement with the result of classical coupled dipole calculations showing that the singlyexcitation manifold contains a large population of modes in the range around 3/Γ 0 [see Sec. IV, Fig. 4(b)].
Finally, we come back to the description of the dynamics in terms of the mean-field model introduced in sec. II. This model assumes a factorizable density matrix throughout the decay, and the coupling between atoms is induced by their dipole moment proportional to the coherence between ground and excited states ρ eg [39]. However for high saturation the atoms are prepared in an incoherent mixture of the ground and excited states, hence the coherence and thus the dipole vanish. Therefore the mean-field model predicts a decoupling of the atoms with one-another, which then decay independently with the single atom lifetime 1/Γ 0 . Our observations up to s = 250, see appendix G, contradict this prediction showing that the density matrix of the system cannot be factorized throughout the decay, although it can be factorized initially.

IV. RELEASE OF THE SUBRADIANT EXCITATION
In this final section, inspired by theoretical proposals [16][17][18], we perform a proof-of-principle demonstration of the ondemand release of the light stored in subradiant excitations. To do so, we apply a position-dependent detuning. The resulting inhomogeneous broadening makes the interaction between atoms no longer resonant: the ensemble now consists of independent atoms efficiently radiating, thus releasing the subradiant excitation. To develop an intuition of how this protocol works, we first consider a toy model consisting of two coupled linear classical dipoles d 1 and d 2 with decay rate Γ 0 and separated by r 12 . The dynamics of the system is governed by the following equations (hereh = 1): where δ is the difference between the transition frequencies of the two atoms and V (r 12 ) is the (complex) dipole-dipole interaction potential [26]. This system has superand a subradiant eigenmodes with decay rates Γ ± = Γ 0 ± 2 Im[V (r 12 )]. For δ = 0 these modes are the symmetric and anti-symmetric combinations v ± ∝ (1, ±1). Hence, if at time t = 0 the dipoles are prepared in v − , the turning on of the inhomogeneous broadening (δ = 0) projects v − onto the new basis provided by the eigenvectors of the matrix in Eq. (1). In the limit δ |V (r 12 )| these are the uncoupled dipoles v 1 = (1, 0) and v 2 = (0, 1), with identical decay rate Γ 0 . The evolution of the two-atom dipole is then given by v(t) ∝ e −Γ 0 t/2 (v 1 − e i δ t v 2 ), recovering the single atom decay rate of the radiated power v 2 (t) ∝ e −Γ 0 t . This toy model thus shows that placing the atoms out of resonance with one another allows the system to radiate again as an assembly of independent emitters.
We report in Fig. 4(a) the results of the experiment obtained applying an inhomogeneous broadening. To do so we turn on the trap light at different times (indicated by the vertical arrows) after the extinction of the excitation laser. The black curve is our reference for which the trap is turned on at long time ( 1 µs). The far off-resonant light induces a positiondependent detuning described by: where δ 0 is the light-shift induced by the trap. Experimentally δ 0 32 Γ 0 , w trap = 2.5 µm and x R = πw 2 trap /λ trap . Using the atomic density distribution in Eq. (2) for the x i , y i , z i , the standard deviation of the detuning induced by the trap is about 4 Γ 0 . At the turn on of this inhomogeneous broadening, we observe the emission of a pulse of light [ Fig. 4(a)]. The presence of this pulse can be qualitatively understood using the toy model: when the inhomogeneous broadening is applied the atoms start to radiate at a rate faster than the subradiant one, thus the intensity of the emitted light is initially enhanced before decaying. Moreover we observe that the pulse is followed by a faster decay similar to the singleatom regime (fitting the experimental data after this pulse we obtain a decay rate of 1.3(1)Γ 0 for all data sets). The measurements have been performed with 5000 atoms and with s 27. We verified the collective origin of this effect by performing the same measurements for low atom number, observing no pulse. Furthermore the strong suppression of subradiance obtained with a light-shift varying slowly in space demonstrates that the subradiant excitations are not stored in pairs of closeby atoms but rather delocalized over all the atoms of the cloud [40]. This is expected for a near-resonant excitation of the cloud, as delocalized excitations corresponds to states with small interaction frequency shifts.
To go deeper in the understanding of the observed behavior, we extend the toy model previously introduced to large ensembles. We obtain the decay rates by evaluating the eigenvalues λ i of the interaction matrix for N = 5000 classical coupled dipoles sampled using the experimental positiondependent detuning. The associated modes are singleexcitation modes, but as observed earlier they should give a qualitative description of the behaviour of the subradiant tail. The results are shown in Fig. 4(b). They indicate that in the presence of the inhomogeneous broadening, the distribution of decay times is much narrower than in free space [note the log scale in Fig. 4(b)]. In particular a significant fraction of modes with decay rates close to the observed subradiant one (between Γ 0 /3 and Γ 0 /10) is suppressed by the inhomogeneous broadening. We further performed numerical simulations of the dynamics for smaller atom numbers using the mean-field model already introduced in Section II. The results are shown in Fig. 4(c) for a cloud with N = 200 atoms at a density ρ 0 /k 3 = 0.3. The trap is turned on during the decay. We took into account the finite rise time of the trap beam of about 25 ns. The simulation provides the evolution of the total population of the excited state ee (t) (dashed line). However, experimentally we measure the intensity of the light emitted by the cloud, which in the absence of a drive is proportional to −d p(t)/dt (in a 4π solid angle) represented by the continuous lines of Fig. 4(c). When the inhomogeneous broadening is applied, the population curve presents a change of slope, corresponding to a peak in the intensity, followed by a decay at a rate now close to the single atom one. These simulations confirm the interpretation of our experimental findings Finally, as shown in Fig. 5, we observe that the enhancement of the emission is stronger along the cloud axis. We have further observed by changing polarization and the internal structure of the atoms that no significant difference arises in the released pulse shape, Fig. 5. We leave for future work the investigation of how to control the directionality of the pulse which could allow for the tailoring of the angular emission pattern of the excitation stored in subradiance.

V. CONCLUSION
In this work, we investigated the subradiant decay of excitations stored in a dense cloud of atoms trapped in an optical tweezer, approaching the Dicke limit of a large number of atoms in a volume smaller than λ 3 . We confirm that the cooperativity parameter is the atom number rather that the optical depth or the density. Moreover, by tuning the intensity of the excitation laser, we studied the nature of multiply-excited subradiant states and experimentally find that all the information is contained in the response of the system at low intensity where only singly-excited modes are populated. A quantitative theoretical description in our regime is however extremely challenging and new models describing this interacting dissipative many-body system need to be developed. Finally, by applying an inhomogeneous broadening, we were able to release the subradiant excitation stored in the cloud in the form of a pulse of light. Our experiment thus represents a proof-of-principle of the temporal control of subradiance in an atomic medium, a prerequisite for its use as a light storage medium. It was made possible by the small size of the cloud, which permitted the use of a relatively low power of the control light to apply a significant detuning between atoms. In the future it could be applied to systems featuring tailored subradiant modes, such as structured atomic systems [18,[41][42][43], in which these modes could be targeted by a local addressing. use of a single exponential decay. We fit all the data in this way and for every measurements we evaluate the reduced χ 2 : where K bin is the number of time bins in the data set (500 ps bins), N exp i is the recorded number of counts in the bin centered on time t i and N the (t i ) is the value of the fitting function. This definition of χ 2 assumes poissonian statistics for the recorded counts in each bin.
In Fig. A2(c) we report the values of χ 2 as a function of the atom number, for the shallowest trapping geometry described in the main text (blue points in Fig. 2). As the atom number in the cloud becomes larger, the system exhibits super-and subradiance and the fit with a single exponential decay is no more able to describe the entire observed dynamics as revealed by an increase in χ 2 . When it becomes larger than 1 we use a different phenomenological model to describe the decay, fitting with the sum of two exponential decays [see Fig. A2(b)]. The agreement between this second model and the experimental data is much better in the large atom number regime, and consequently χ 2 is smaller as highlighted by the filled points in Fig. A2(c). The very small values of χ 2 comes from the large uncertainties resulting from the small number of counts in the tail.

Appendix D: The three different trapping geometries
In section II, we employ three different trapping geometries that allow us to obtain three different cloud dimension. In the first one, the trap is the one used to load the atoms from the MOT (Gaussian sizes l r 0.7 λ , l x 7.7 λ blue circles). In the second one we increase the power of the trapping beam by 60% gaining a factor 2 in density without losing atoms (l r 0.5 λ , l x 6.0 λ , purple diamonds). In the third case, we reduce the beam waist down to 1.8 µm. The atom number after the compression is lower (N 1500) but we reach a high density ρ 0 /k 3 1.5 owing to the reduced trapping volume (l r 0.4 λ , l x 2.9 λ black squares), [32]. A complete set of time traces can be found below. For every point we measure the atom number and the temperature.

Appendix E: Effect of polarization and optical pumping
Unless otherwise specified, the measurements reported in this work are performed with a linearly polarized excitation light without any direct optical pumping (OP). We did not observe that the polarization direction impacts in any way the observed subradiance. However, there could be an effect due to the fact that the linear polarization configuration does not correspond to a two-level configuration.
In order to check if this strongly modifies the collective behavior hosted in the tail of the fluorescence signal, we performed one experiment in a situation much closer to a twolevel configuration: a circularly-polarized (σ − ) excitation, together with a 20 G magnetic field along the excitation propagation axis. We perform hyperfine and Zeeman optical pumping with the same polarization as the excitation light to place ourselves in a closed two-level system (between F = 2, m F = −2, F = 3, m F = −3). The light scattered by an atom on a nearby one might contain other polarization components and drive other transitions out of the 2-level system. However the detuning between the σ − and π polarization transitions (closest in detuning) is about 5 Γ 0 preventing the rescattering of this polarization by nearby atoms.
In Fig. A3 we report in blue an acquisition done with linear polarization and without optical pumping, and in purple that with optical pumping, circular polarization and magnetic field. We observe that the tails of the two traces are similar, within our dynamic range of observation. We conclude from this that the internal structure does not seem to affect strongly the results on subradiance reported in this work and our conclusions. Time (ns) Normalized counts FIG. A3. Photon count decay with a linear polarization, no magnetic field and no optical pumping i.e. multi-level situation (blue). The same in the two-level case: 20 G magnetic field, σ − polarization and prior optical pumping (purple). The two traces have been normalized to the steady state, measurements for N 4500 in the first trapping geometry (see appendix D).

Appendix F: Relative population of the long-lived states
Another parameter that we use to quantify subradiance is the relative fluorescence observed in the long-lived tail with respect to the total fluorescence recorded after the switch off of the drive. We refer to this quantity as the tail ratio. It is defined as: where I(t) represents the time-resolved fluorescence emitted by the atomic cloud where t = 0 is the switch off time of the excitation laser. This parameter estimates the fraction of excitation still hosted in the system 4/Γ 0 100 ns after switching off the excitation, compared to the whole decay. In the single-atom limit, this parameter is equal to e −4 2% by definition.
We analyze the behaviour of this tail ratio as a function of atom number on the same data as the one used for Fig. 2. In Fig. A4 we indeed see the collapse as a function of atom number similarely to what found for the decay time (see discussion in Sec. II).

Appendix G: Measurements at large saturation intensity
In this section we report the data acquired at very high values of the drive intensity, reaching s = I/I sat 250. The measurements have been performed with a larger atom number, N 6000. For this reason this dataset and the one used for Fig. 3 of main text cannot be directly compared. In Fig. A5(a) we report the temporal traces acquired in this way, normalized to the steady state value of the trace taken with the largest excitation intensity. In Fig. A5(b) we instead report the results of numerical simulations performed using the nonlinear coupled dipoles model using 100 atoms distributed in a Gaussian cloud with a peak density ρ 0 = 0.3/k 3 . The reported lines are the results of 10 realizations of the same numerical experiment. As one can see, according to the meanfield model, the subradiance is expected to disappear as the excitation strength is increased. The fact that experimentally the system still hosts a subradiant excitation even at very large intensity shows that the density matrix of the system cannot be factorized throughout the decay.