Ana Cipris

Dr William Guerin

Professeur Philippe Courteille

Patrizia Vignolo

Sergey Skipetrov

Maître Vincent Josse

Dr Robin Kaiser

Beatriz Olmos

Nathalie Hamel

Francois-Regis Huc

Isabelle Larochette

Sanaë Mahir

Modes collectifs de longue durée de vie de la lumière dans des atomes froids: études expérimentales et numériques sur la sous-radiance et la localisation d'Anderson

Keywords: Diffusion collectif, atomes froids, sous-radiance, localisation d'Anderson vi Collective scattering, cold atoms, subradiance, Anderson localization viii Teti

Modes collectifs de longue durée de vie de la lumière dans des atomes froids: études expérimentales et numériques sur la sous-radiance et la localisation d'Anderson L'objectif de cette thèse est d'élargir notre compréhension des phénomènes collectifs de diffusion de la lumière qui se manifestent par la forte suppression du rayonnement du milieu diffusant. En particulier, nous nous intéressons à la sous-radiance, un effet collectif de l'émission spontanée, qui, du fait du couplage dissipatif des atomes via des modes électromagnétiques communs, crée des modes d'excitation longue durée de vie dans un ensemble d'atomes froids, avec des durées de vie beaucoup plus longues que celles d'atomes individuels. Auparavant, la sous-radiance avait été bien caractérisée pour des ensembles d'atomes froids dilués et dans la limite de faible intensité du laser d'excitation (régime de l'optique linéaire).

Dans cette thèse, nous étudions ce phénomène collectif au-delà de l'optique linéaire et du régime dilué.

Pour étudier expérimentalement la sous-radiance au-delà du régime de l'optique linéaire, nous avons fait varier le paramètre de saturation du laser d'excitation jusqu'à des valeurs relativement importantes et nous avons observé l'augmentation de la population sous-radiante avec l'augmentation du paramètre de saturation. Plus important encore, nous avons identifié une loi d'échelle super-linéaire de la population sous-radiante au-delà du régime de l'optique linéaire, en raison d'un processus similaire au pompage optique via des états super-radiants multi-excités. Cette étude a montré que l'augmentation de l'intensité du laser de pilotage au-delà du régime de l'optique linéaire permet de surmonter le faible couplage des modes sous-radiants avec l'environnement.

De plus, nous avons étudié numériquement la sous-radiance en utilisant le modèle vectoriel des dipôles couplés, en mettant l'accent sur la compréhension de l'effet de l'interaction dipôle-dipôle en champ proche. Pour cela, nous avons comparé en détail les modèles vectoriel et scalaire (qui ne contient pas de termes d'interaction en champ proche) et étudié la durée de vie sous-radiante en fonction de la densité et de l'épaisseur optique du nuage d'atomes froids.

Nous avons observé que les termes de champ proche sont préjudiciables à la sous-radiance, ce qui se manifeste par la réduction des durées de vie sous-radiantes avec l'augmentation de la densité. Ce caractère néfaste des termes de champ proche peut être interprété comme étant dû à un élargissement inhomogène effectif induit par l'interaction en champ proche, qui réduit les durées de vie sous-radiantes lorsque la densité augmente.

Le caractère néfaste de l'interaction en champ proche a également été précédemment identifié numériquement pour la localisation d'Anderson des ondes lumineuses (la localisation v d'Anderson est un phénomène ondulatoire en diffusion multiple cohérente dans un milieu fortement diffusant, qui entraîne l'arrêt du transport diffusif ; il donne également des modes à longue durée de vie) : le modèle scalaire des dipôles couplés prédit la localisation des ondes lumineuses, mais pas le modèle vectoriel. Cependant, l'ajout d'un champ magnétique suffisamment fort permet de récupérer les modes localisés pour la lumière vectorielle. Dans ce contexte, nous avons étudié numériquement une possible signature de la localisation d'Anderson de la lumière dans les atomes froids, basée sur la mesure des fluctuations d'intensité de la lumière diffusée. Nous avons montré qu'à la transition vers le régime localisé, les fluctuations sont augmentées. En particulier, pour la lumière vectorielle, nous avons démontré l'absence, et l'apparition, de fluctuations d'intensité accrues sans, et avec, le champ magnétique appliqué, ce qui corrobore que les fluctuations d'intensité pourraient servir de moyen expérimental pour sonder les états localisés.

Table of Contents

Introduction

Wave propagation and scattering are natural phenomena that are ubiquitous in everyday life;

the sounds that we hear, the ocean waves, ripples on a pond, earthquakes, the blue color of the sky, and colorful sunsets are only a few examples of the phenomena resulting from propagation and scattering of many types of waves in various media, which we experience everywhere around us. The processes such as scattering, diffraction, or interference are fundamental for all waves (electromagnetic, sound, seismic, matter waves, etc.), however, the properties of these processes and the phenomena that they might give rise to as a result of the wave-matter interaction, can greatly depend on the medium itself.

The research on propagation and scattering of waves in various media have been of fundamental, as well as practical interest for many years (to be able to understand and exploit the wave-related phenomena). In particular, the wave propagation in a random, i.e. disordered, medium has been extensively studied in the last decades, especially due to the discovery of coherent multiple-scattering phenomena, such as weak and strong localization, and the realization that these phenomena may be common to all types of waves [START_REF] Akkermans | [END_REF]].

The main research activity of the cold-atom group at Institut de Physique de Nice has been the study of phenomena related to the scattering of light waves using a very specific random medium of scatterers: laser-cooled and trapped atoms. In general, the cooling, trapping, and manipulation of atoms with laser light have revolutionized atomic physics and opened up an impressive field of research, with the goals of fundamentally understanding atomic and optical processes, as well as using cold-atom systems as platforms in many different research areas. The advances in laser cooling and trapping techniques [Metcalf 1999[START_REF] Adams | [END_REF]] and higher as in the case of Bose-Einstein condensates), have provided scientific breakthroughs in many areas such as precision measurements (atomic clocks, high-resolution spectroscopy, atom interferometry), quantum degeneracy (Bose-Einstein condensates) and quantum simulations.

Furthermore, cold atoms are exceptional systems to study light scattering and they allowed for the observation of many phenomena related to light scattering. There are many properties of cold atoms that make them very often more advantageous for studying phenomena related to light scattering, in comparison to classical scattering objects. The control and manipulation of light scattering in a cold-atom medium can be easily achieved by modification of properties of the cold-atom sample, such as its on-resonance optical depth, density, and temperature, via changing the parameters for laser cooling and trapping. Due to the well-known quantum internal structure of atoms and the well-understood light scattering properties of individual atoms, light scattering in cold atoms can be further controlled in a very predictable way by, for example, modifying the detuning of the incoming light from the resonant frequency of the transition. Also, the atoms' response to the light can be modified by applying magnetic or electric external fields. Moreover, atoms can be considered as point scatterers (for visible light, the wavelength λ is much larger than the atom size) and hence, they are very strong resonant scatterers: the resonant cross-section of each individual atom is very large, on the order of λ 2 . In fact, with cold atoms there is no significant effect on the cross-section due to atomic motion, since the Doppler-broadening, which can reduce the resonant cross-section, is nearly absent in cold-atom media. As opposed to classical scatterers, this large resonant scattering cross-section of cold atoms enables obtaining a strongly scattering medium even at low densities. Additionally, cold atoms have very narrow resonances (linewidth), which allows us to neglect the other atomic levels not corresponding to the transition of interest, and hence the cold atoms can often be considered as degenerate two-level systems. Another features of cold-atom media are monodispersity, i.e., atoms are all identical scatterers having the same resonant frequency and linewidth, and the absence of absorption (each absorbed photon is re-emitted, which basically makes this process a scattering event), which is important since with classical scattering media, the absorption of light turned out to be a big obstacle in observing coherent multiple-scattering phenomena. All these properties make the study of multiple scattering in cold atoms, and more generally of collective effects in light-atom interaction, particularly rich [Kaiser 2000].

itself as a breakdown of a diffusive wave transport. In other words, for a strong enough disorder, the collective interference in path amplitudes of multiply scattered waves brings the extended (diffusive) wave propagation to a complete halt and leads to the exponential spatial localization of waves. In fact, in low-dimensional disordered samples, waves are always localized no matter what is the disorder strength (as long as the size of the medium is larger than the localization length). However, there are two regimes, i.e. phases, of wave transport in a 3D disordered medium: diffusive propagation of waves in a weakly scattering medium and localized regime for waves in a sufficiently strongly scattering medium. Hence, for 3D disordered media, a phase transition from extended-wave to localized-wave regimes exists, i.e.

a critical amount of disorder is required for AL to occur. The Anderson localization of waves in 3D random media, i.e. strong localization, has been observed for several types of waves.

However, despite many experimental efforts (mainly with classical scattering media, such as semiconductor powders or dielectric particles), strong localization has not been observed

for light waves [Skipetrov 2016b]. A cold-atom disordered medium seems to be a promising platform for observing strong localization of light due to its specific features, such as the monodispersity of scatterers, the large scattering cross-section of cold atoms, the relatively easy manipulation as well as control of the cold-atom medium and its scattering properties.

A more detailed discussion on obstacles with classical scatterers, as well as advantages (and disadvantages) of cold atoms for light localization will follow in Chapter III.

On the way towards Anderson localization of light in 3D, different aspects of multiple scattering of light in cold atoms have been extensively studied by the group, like the slow diffusion of light due to classical multiple-scattering (which can be described by the photon random-walk model) [Labeyrie 2003] and coherent backscattering, a coherent multiplescattering phenomenon which is very often considered to be a precursor of strong localization [Labeyrie 1999, Bidel 2002]. For a review on these early works, see [Labeyrie 2008].

All these experiments were done with very dilute cold-atom samples, where the average distance between atoms is much larger than the wavelength. However, it is expected that Anderson localization occurs in the regime of very strong scattering, such that the mean-free path between two scattering events is smaller than the wavelength. Thus, to go towards such small photon mean-free paths, a straightforward way is to simply decrease the distance between scatterers by increasing the density of the atomic sample. However, it is known in quantum optics that two very close atoms (whose interatomic distance is smaller than the wavelength) will form superradiant and subradiant collective states [Dicke 1954,Stephen 1964].

For example, the influence of those superradiant pairs on the light transport properties has been discussed in Ref. [Gero 2006].

This led us to consider also the phenomenon of subradiance, which creates long-lived collective states due to the anti-synchronizations of dipoles (the two-level atoms are considered as dipoles, i.e. tiny antennas that emit electromagnetic waves): the electromagnetic fields radiated by such dipoles will collectively undergo a destructive wave interference which leads to the inhibition of radiation. Since Anderson localization can also create long-lived modes in the strongly disordered medium, it appears that both effects might exhibit similar signatures, such as a slow fluorescence decay. Therefore it seems necessary to understand subradiance in order to be able to distinguish it from Anderson localization of light waves. Let us emphasize that both subradiance and Anderson localization are wave interference phenomena. In fact, the subradiance (superradiance) is a global interference phenomenon, i.e. it is a collective effect involving anti-synchronization (synchronization) of dipoles in the whole medium (superand subradiant states are extended over the entire sample). On the other hand, Anderson localization is a local interference phenomenon: in a strongly enough disordered medium, the light is confined inside the medium by interferences in multiple scattering from only a few scatterers (the AL modes are localized on those few scatterers).

While strong localization requires a very strongly scattering medium, this is not mandatory for collective effects, such as super-and subradiance, to arise. As we will see in the next paragraph, super-and subradiance are not limited only to dense atomic samples, but they can also appear in the dilute limit. Our group has thus started a research program on superand subradiance, first in the limit of dilute cold-atomic samples.

Brief introduction on super-and subradiance

Collective effects in spontaneous emission: Dicke's superradiance and early experiments

When it comes to the scattering from an ensemble of atoms, it is generally incorrect to assume that the emission of each atom occurs independently [such that the fluorescence decay from an ensemble of atoms is ∝ exp(-Γt), characterized by the emission rate of individual atoms Γ = 1/τ (where τ is the natural lifetime of the excited state)]. Since the scattering is a wave phenomenon, the collective interferences in spontaneous emission might drastically modify the emission from a collection of scatterers. The collective response of atomic ensemble to the light may give rise to superradiant (subradiant) phenomena, which, due to the constructive (destructive) interferences in emission, exhibit enhanced and accelerated (suppressed and slowed down) collective emission rate compared to that of a single atom.

The collective phenomena in light emission, mainly superradiance, have been predicted by R. Dicke in 1954 in his pioneering work "Coherence in spontaneous radiation processes" [Dicke 1954]. By considering an ensemble of N motionless, identical, two-level atoms located in a sub-wavelength volume (R λ, where R is the atomic sample size and λ is the transition wavelength; hereafter, such sample is referred to as small sample) and addresses the question of the collective eigenstates of the system and their radiation rates (see Fig. 1). In particular, Dicke showed that if all atoms are initially in their excited state (fully inverted symmetric state), the de-excitation of the ensemble would follow the cascade of all symmetric, i.e.

superradiant states (superradiant cascade), with an acceleration of the emission rate (for example, in the case of the symmetric single-excitation state the emission is predicted to be N times faster compared to the single-atom emission rate Γ) and a coherent emission.

Anti-symmetric states are also introduced, which can be completely decoupled from all electromagnetic modes and hence, exhibit extremely long lifetimes (subradiance). The phenomena of super-and subradiance can be explained as follows. The interaction between electric dipoles of individual atoms due to their coupling to common vacuum field modes (collective coupling of dipoles), triggers the synchronization of dipoles and classical interference processes, which can be envisioned as a build-up of a collective dipole.

Superradiance stems from in-phase oscillations of dipoles (constructive interference) which lead to the build-up of a huge collective dipole (larger than the dipole of individual atoms).

The emission of such a huge collective dipole thus exhibits an accelerated decay. In other words, superradiance yields short-lived collective states of the atomic ensemble, with a maximum emission rate Γ sup larger than that of an individual atom Γ (Γ sup > Γ). On the other hand, subradiance occurs when dipoles are anti-synchronized, i.e. they oscillate out of phase (destructive interference), which yields a tiny collective dipole (much smaller than the individual dipole). Radiation from such tiny dipole is suppressed, i.e. the fluorescence decay from an ensemble of atoms is immensely slowed down. Hence, subradiance results in creation of long-lived collective states, whose lifetimes τ sub = 1/Γ sub can be much longer than those of individual atoms τ (τ sub > τ ).

The pioneering work by Dicke triggered an extensive theoretical research on collective spontaneous emission addressing various aspects of super-and subradiant phenomena (see, for example, [Stephen 1964, Freedhoff 1967, Ernst 1968, Ernst 1969, Lehmberg 1970a, Lehmberg 1970b, Arecchi 1970, Stroud 1972, Friedberg 1971, Friedberg 1972, Friedberg 1973, Saunders 1973a, Saunders 1973b[START_REF] Milonni | [END_REF][START_REF] Bonifacio | [END_REF], Ressayre 1976, Ressayre 1977, Gross 1982, Crubellier 1985, Crubellier 1986, Crubellier 1987, Men shikov 1999]).

The experimental research on these phenomena started in the 1970s. There have been many observations of superradiance and some of the first are [Skribanowitz 1973[START_REF] Gross | [END_REF], Gibbs 1977] (for a review see [START_REF] Feld | Coherent Nonlinear Optics[END_REF]). In most of the research done in that period, the situation of the Dicke limit, i.e. small sample, has not been reached: the samples have finite size, they are often pencil-shaped, and propagation/geometric effects are important (see, e.g., [MacGillivray 1976]). In that case, the parameter of the sample which governs the superradiant decay rate is related to the optical thickness (see, for example, [Friedberg 1976]).

When it comes to subradiance during that period (before the introduction of the concept of single-photon superradiance, which will be discussed below), there have been only two reports on observation of subradiance: [Pavolini 1985] (indirect signature with a large number of atoms) and [DeVoe 1996] (two-ion system).

Single-photon superradiance

A revival of the topic of super-and subradiance began in 2006 with a publication of the study by M. Scully and collaborators [Scully 2006] that addresses the question of collective spontaneous emission from an extended (R > λ) ensemble of N two-level atoms prepared by a single-photon excitation, i.e. by absorbing one photon of wave vector k 0 : there is only one among N atoms which is excited but we do not know which one. The atomic sample can then be described by the collective state (which was coined as Timed-Dicke (TD) state):

|Ψ = 1 √
N j e ik 0 •r j |g 1 , g 2 , ..., e j , ..., g N .

The r j are the atomic positions and |g j , |e j are the ground and excited states, respectively. than the single atom decay (black dashed line) and then slower (subradiance). The inset shows the emission diagrams at short times (blue) and long times (red) after the laser switch-off. Taken from Ref. [START_REF] Bienaimé | [END_REF].

By studying the spontaneous emission from this collective state, they have shown that the emission of a photon occurs in the same k 0 direction, which might seem counterintuitive for spontaneous emission. Moreover, there is no condition on the density of the atomic sample, the interatomic distance does not have to be small compared to the wavelength for this collective effect to occur. This study triggered a lot of discussions [Eberly 2006[START_REF] Mazets | Multiatom cooperative emission following single-photon absorption: Dicke-state dynamics[END_REF], Scully 2007, Das 2008]. Note that the whole subject has been given the term 'single-photon superradiance' [Scully 2009a, Scully 2009b] (while for its counterpart, which was addressed later on, 'single-photon subradiance' [Scully 2015]).

Following that study, the superradiant emission rate of the TD state was determined to be [START_REF] Mazets | Multiatom cooperative emission following single-photon absorption: Dicke-state dynamics[END_REF],Svidzinsky 2008,Svidzinsky 2010,Courteille 2010,Friedberg 2010,Prasad 2010]

Γ sup C N (kR) 2 Γ ,
where C is a numerical factor depending on the geometry. This collective emission rate is larger by a factor of ∝ N/(kR) 2 from the single-atom decay rate Γ. The amplification factor is therefore neither the number of atoms (as in the Dicke limit of a small sample), nor the density of a sample. In fact, for a spherical sample, N/(kR) 2 is proportional to the on-resonance optical thickness of the sample b 0 .

Following those predictions on the superradiant decay rate of the TD state, it has been numerically demonstrated by our group that, for a large and dilute atomic sample, a small fraction of the excitation would populate subradiant states, giving rise to a decelerated emission rate, whose lifetime is also proportional to the on-resonance optical thickness [START_REF] Bienaimé | [END_REF] (see Fig. 2).

Note that instead of using a single-photon source, the study of the collective spontaneous emission from large atomic samples can be carried out with a continuous laser field in the weak-excitation limit, i.e. for very low saturation parameter of the laser, which would allow for the assumption that only one out of N atoms is excited [Courteille 2010, Prasad 2010].

For that case, the better terminology is 'super-and subradiance in the linear-optics regime' than 'single-photon super-and subradiance'.

On the theoretical side, the calculations and numerical simulations addressing collective spontaneous emission phenomena in the low-excitation limit (linear-optics regime) can be obtained with the so-called coupled-dipole model (which will be introduced in Chapters II and III), which treats atoms as classical dipoles and it considers the dipole-dipole interactions mediated by the radiated field (each atom is coupled to other atoms through the electric field radiated by all the other atoms). Besides super-and subradiance [Bienaimé 2011a[START_REF] Bienaimé | [END_REF], Bienaimé 2013[START_REF] Roof | Observation of Single-Photon Superradiance and the Cooperative Lamb Shift in an Extended Sample of Cold Atoms[END_REF], Guerin 2016, Araújo 2016, Araújo 2018b], such calculations involving collective dipole-dipole coupling have been widely used to study various phenomena related to light-scattering from cold-atom media (samples of motionless atoms), such as the collective radiation pressure force [Courteille 2010, Bienaimé 2010, Bux 2010, Bienaimé 2014, Bachelard 2016], the collective energy shifts [Scully 2009a, Bienaimé 2011a, Svidzinsky 2010, Keaveney 2012, Meir 2014[START_REF] Javanainen | [END_REF], Jenkins 2016[START_REF] Roof | Observation of Single-Photon Superradiance and the Cooperative Lamb Shift in an Extended Sample of Cold Atoms[END_REF]], the collective line broadening [Zhu 2016,Sutherland 2016,Bromley 2016,Jennewein 2016], collective Mie scattering [Bender 2010, Bachelard 2012], and coherent backscattering [Rouabah 2014, Chabé 2014], as well as photon localization [Bellando 2014[START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF], Máximo 2015].

Previous experimental results of the group: super-and subradiance in the linear-optics regime in dilute clouds

Prompted by the theoretical research on single-excitation super-and subradiance, i.e. superand subradiance in the linear-optics regime, the experimental research of the cold-atom group in Nice resulted in 2016 in the first experimental observation of subradiance in extended (R λ) and dilute cold atomic clouds containing N 1 atoms, as well as in the weakexcitation limit (with the saturation parameter of the driving laser s ∼ 10 -2 ) [Guerin 2016].

From the temporal decay dynamics of cold-atom samples at late times (light scattered off-axis , with a detuning δ = -3.15 (in units of natural linewidth) of the driving laser. All decay curves are normalized to the fluorescence level right at the switch off of the laser. For comparison, the theoretical single atom decay exp(-t/τ at ) is also shown (dash-dotted line). The longest subradiant lifetime τ sub is obtained from an exponential fit (dashed lines) at late time. (b) Measured subradiance lifetimes τ sub /τ at as a function of b 0 . All measured points collapse on a single line, independently of the detuning. The linear scaling of τ sub with b 0 is stressed by the linear fit (solid line). Taken from Ref. [Weiss 2018].

after the switch-off of the driving laser), subradiant lifetimes as long as τ sub ∼ 100τ (τ = 1/Γ the decay time of an individual atom) have been measured. More importantly, it has been shown that the subradiant lifetime scales linearly with the on-resonance optical thickness of the cloud b 0 ∝ N/(kR) 2 and that the subradiant lifetimes are independent of the laser drive detuning (see Fig. 3). That demonstrated that the relevant cooperativity parameter is b 0 , i.e. that the subradiance, in the limits of dilute and extended atomic sample, is governed by the resonant optical thickness of the scattering medium.

Then followed the experimental study on off-axis superradiance in the dilute, extended and weak-excitation regime [Araújo 2016]. The observation of superradiance has been reported for large detuning and it has been demonstrated that the superradiant decay rate increases with the on-resonance optical thickness, but it is suppressed in the case of near-resonant excitation (it depends on the detuning of the laser drive) due to attenuation and multiple-scattering.

Later on, the difference between subradiance and radiation trapping, in the dilute and extended regimes of cold-atom cloud, has been discussed and experimentally characterized in [Weiss 2018]. Since the radiation trapping also gives rise to a slow fluorescence decay (compared to the decay of an individual atom) due to (classical) multiple scattering (which can be described by the random walk model), it was of interest to understand the interplay between the two phenomena exhibiting a slow decay. By illuminating the cold-atom cloud by a laser beam whose waist is smaller than the size of the cloud (to eliminate the dominating single-scattering contribution from the edges of the sample; this is in contrast to previously mentioned experiments where the beam is larger than the cold atomic cloud), the two phenomena were observed simultaneously but with different properties: the early decay is dominated by the radiation trapping and the scaling of the corresponding lifetimes with detuning-dependent optical thickness b(δ) has been observed, while at later times subradiance dominates with subradiant lifetimes independent of detuning and showing the trend of scaling with the on-resonance optical thickness b 0 .

The features of collective effects, such as super-and subradiance, as well as multiple scattering, can be explained in the framework of the coupled-dipole model describing the atoms as dipoles which share one excitation and which are coupled to each other through the dipoledipole interaction generated by the radiated field. In particular, the scalar approximation of the coupled-dipole model, which neglects the polarization of the driving field as well the near-field dipole-dipole interaction terms (it only contains long-range interaction terms ∝ 1/r) has captured the experimentally observed characteristics of subradiance and superradiance in the dilute, extended and weak-excitation limit (see Refs. [Guerin 2017a, Guerin 2017b, Araújo 2016,Araújo 2018b] and Supplemental Material of Ref. [Guerin 2016]). This suggested that the scalar coupled-dipole model may be suitable for simulating light scattering from atomic samples in the dilute limit where the near-field interaction is negligible. This validity for dilute limit will be demonstrated in Chapter II)

The fact that the super-and subradiance are governed by the on-resonance optical thickness b 0 of the dilute and extended scattering medium has been discussed, for example, in Ref. [Guerin 2017a]. The argument is following. Although there is an infinite amount of vacuum electromagnetic modes in free space, for a scattering sample of a finite size R, the number of electromagnetic modes efficiently coupled to the sample is M ∝ (kR) 2 . Then, the ratio N/M (where N is the number of atoms of the sample) represents the number of atoms that emit in the same mode and if this number is N/M > 1, interference effects in the emitted light will occur, giving rise to super-and subradiance phenomena. Moreover, this ratio corresponds to the optical thickness (up to a prefactor) N/M ∝ N/(kR) 2 ∝ b 0 , which provides an explanation for the scaling with b 0 . On the other hand, for the small samples with R λ, there is only one outgoing mode M = 1 (spherical wave). Hence, the relevant parameter for collective phenomena in the limit of small scattering medium is the number of atoms N . A recent experiment demonstrates such behavior [START_REF] Ferioli | [END_REF]].

Content of this thesis

During my doctoral research, the main objective was to study subradiance beyond the limits of zero temperature of the atomic sample, low intensity of the driving laser, and low sample density, in order to extend our understanding of these collective long-lived modes. For that I contributed to several experimental and numerical studies on subradiance. I have also contributed to a study on identifying a signature of Anderson localization of light in 3D disordered atomic media. Moreover, I have participated as well in the study of superradiance and characteristics of collective phenomena at the laser swith-on. These studies led to the several published articles (see the list in Appendix A). Note that in this manuscript I will not present all of the studies, but only those of which I had a greater contribution.

On the experimental side, we first studied the effect of the temperature of the atomic samples on subradiance. We have observed that, when increasing the temperature of the cold-atom cloud in the range from several tens of µK to several mK, the subradiant lifetimes undergo only a slight decrease, and we measured subradiant lifetimes longer than the characteristic time corresponding to the Doppler broadening. This demonstrated the robustness of subradiance against thermal motion. This work has been done with Patrizia Weiss (post-doc on the experiment at the time) and it led to the publication of the article [Weiss 2019].

Then followed the experimental study of subradiance beyond the linear-optics regime. We have demonstrated that the population of subradiant modes is enhanced when the intensity of the driving laser is increased. The superlinear scaling of the subradiant population beyond linear-optics regime is due to a process similar to optical pumping via multi-excitation superradiant states. This work has been done in collaboration with the theoretical group led by Romain Bachelard, which provided the numerical results supporting our experimental findings. The results have been published in [Cipris 2021a]. This study will be described in Chapter I.

During that study, we have also analysed the behavior of the collective Rabi oscillations at the switch-on of the driving field, demonstrating that collective effects, previously studied after the laser switch-off (temporal decay dynamics), exhibit characteristic behaviors at the laser switch-on. The full study was published split in two papers, one for the linear-optics regime, where the experimental observation of the collective multimode vacuum Rabi splitting in free space that scales with resonant optical thickness has been reported [Guerin 2019] (from the data acquired for [Araújo 2016]), and the other for the nonlinear optics regime [do Espirito Santo 2020]. Since I had a minor contribution to this work (I experimentally acquired the new data), I do not include that study in the manuscript.

On the numerical side, I have developed a code based on the vectorial coupled-dipole (CD)

model to be able to investigate how the near field dipole-dipole interaction (which is absent in the scalar coupled-dipole approximation) affects the observable features of the long-lived modes beyond the dilute limit of the atomic sample, since the near-field interaction is nonnegligible in dense atomic samples and it becomes more significant as we increase the density of the atomic sample and hence, decrease the interatomic spacing. One numerical study with vectorial coupled-dipole model was focused on the subradiance, while the other on Anderson localization of light.

I have also done numerical simulations with the scalar coupled-dipole code for the study of superradiance, in particularly to benchmark the superradiant decay rate obtained from CD simulations to that obtained from another approach introduced by I. M. Sokolov, which provided us an optical picture for superradiance in the linear-optics regime, based on a single scattering event embedded in a dispersive effective medium composed by whole atomic sample. This work has been published as [Weiss 2021].

We have numerically investigated the effect of the near-field dipole-dipole interaction on subradiance by using the vectorial CD model. We observed that the near-field effects have negligible role on subradiant lifetimes in the low-density regime, however for higher densities, we observed that subradiant lifetimes decrease with increasing density for a given on-resonance optical thickness of the atomic sample. These results have been benchmarked with the results obtained from the scalar coupled-dipole model to confirm that the identified decrease of subradiant lifetimes is absent in the scalar-approximation, implying that this decrease is indeed due to the near-field interaction. This work was published in [Cipris 2021b] and it is presented in detail in Chapter II. Chapter I

Subradiance beyond the linear optics regime

The collective effects in spontaneous emission from an ensemble of scatterers may yield subradiant states, i.e. collective atomic states whose lifetimes are orders of magnitude larger than the lifetime of a single-atom excited state. For example, subradiance has gained a lot of interest in quantum memory and quantum information, and there have been many proposals for storage of light through preparation of subradiant states [Kalachev 2007, Scully 2015[START_REF] Plankensteiner | Selective protected state preparation of coupled dissipative quantum emitters[END_REF], Facchinetti 2016, Jen 2016, Asenjo-Garcia 2017, Jen 2017, Guimond 2019, Needham 2019].

Subradiant modes are essentially weakly coupled to the external electromagnetic fields, which makes them hard to excite. In this chapter we experimentally study subradiance beyond the linear-optics regime by increasing the saturation parameter of the driving laser up to relatively large values. We demonstrate that the weak coupling of subradiance in large and dilute cold-atom samples can be overcome in the beyond linear-optics regime by driving the cold-atom sample with a strong excitation laser.

This chapter is organized as follows. In Sec. I.1, we introduce our cold-atom sample prepared in a magneto-optical trap, as well as experimental setup and protocols for measuring the slow fluorescence decay and properties of cold-atom samples. Then, in Sec. I.2, we present the experimental results on subradiance obtained by driving the sample with a probe beam of varying saturation parameter, which demonstrated the super-linear scaling of subradiant population with increasing saturation parameter, while the final subradiant lifetimes remain unaffected by the strength of the laser. In Sec. I.3, we briefly introduce theoretical models and we compare the numerical results to those obtained experimentally. This theoretical part was obtained by the group of Romain Bachelard. Moreover, we discuss the interpretation of the results. Finally, in Sec. I.4 we summarize and comment on our results.

13 (dilute cloud; ρ 0 λ 3 ∼ 10 -2 ), while its size, i.e. RMS radius is R ≈ 1 mm. Moreover, the number of atoms of the generated cold cloud is N ≈ 6 × 10 9 , and the obtained on-resonant thickness b 0 = σ 0 ρ(0, 0, z)dz (for the laser beam whose propagation direction is taken to be along ẑ axis), where σ 0 is resonant rubidium atomic cross-section, is up to b 0 ∼ 180.

To study the behaviour of the slow fluorescence decay, we excite the atoms by a probe laser on the hyperfine transition F = 2 → F = 3. As our focus is on what is happening after the laser switch-off, there are a couple of important requirements for observing the slow decay, i.e. the decay of the long-lived collective (subradiant) modes. First, we need a fast switch-off: subradiant modes have lifetimes longer than the lifetime of an individual atom, hence a switch-off duration shorter than τ = 26.24 ns is required so that the subradiant decay would not be masked by the laser switch-off and we would undoubtedly be able to observe the slow decay associated to the cold atomic ensemble. Moreover, subradiant modes are hard to populate due to their weak coupling to the environment and consequently the fluorescence level corresponding to subradiant decay is expected to be very low. Hence, we need a very good extinction ratio of the probe beam to be able to observe the slow decay of subradiant modes. These requirements are met by using two acousto-optical modulators (AOMs) in series, which allow us to create probe pulses (in duration of 5 µs) with a switch-off duration of ≈ 15 ns (90% -10% fall time) and an extinction ratio better than 10 -4 . For a detailed description of the setup see Ref. [Araújo 2018a].

For generating the MOT trapping beams, as well as the probe beam, we use an external cavity diode laser (Toptica DL pro with FWHM linewidth of 500kHz). Additionally, in the experimental setup we have implemented a tapered laser amplifier (Toptical BoosTA pro;

with an input and output fiber coupling) which increases the output power of the DL pro master laser. This optical amplifier allows us to have enough power for the probe beam to explore the behaviour of subradiance for a large range of probe saturation parameter. It should also be noted that we have implemented a series of repumping laser pulses in the experimental cycle: we added a repumper pulse just before each probe pulse is applied (except for the first probe pulse), which ensures that all the atoms are optically pumped to F = 2 hyperfine ground state after applying each probe pulse. These repumper pulses have been added specifically for this experiment where the probe saturation parameter is varied up to relatively high values.

I.1.2. Experimental cycles and measurement protocol

Absorption imaging cycles

In previous experiments of the group on subradiance [Guerin 2016, Weiss 2018] (obtained in the weak excitation limit), the measurements of the properties of the cold 87 Rb cloud corresponding to each applied pulse were done after the data acquisition of the decay of the scattered light (see [Araújo 2018a] for more details on such calibration of cloud properties).

However, in this experiment we probe the atomic sample with a varying saturation parameter, up to relatively large values. Therefore, each applied pulse could significantly affect the atomic sample, for example due to heating. It is therefore important to characterize the atomic sample taking into account the effect that the probe pulses have on the cloud. To do so, we have implemented a method that allows us to measure the properties of the atomic sample simultaneously to the data acquisition of fluorescence decay. Note that here simultaneously does not mean that measurements of the sample properties are done exactly at the same moment as the data acquisition of fluorescence decay, but that the measurements of atomic sample are interlaced with the data acquisition by changing one out of 250 cycles, as it is going to be explained below. Such measurement of the atomic sample also allows us to identify any possible drifts of the sample properties during long data acquisition.

Before describing the method for the sample characterization simultaneously to the data acquisition, let us mention that the cloud properties, such as the on-resonance optical thickness and temperature of the cloud, are obtained from absorption imaging measurements. From the absorption imaging we obtain the transmittance of the probe beam T which is related to the optical thickness by T (x, y) = exp [-b(x, y)], where

b(∆) = b 0 1 + 4(∆/Γ) 2 , (I.1)
is the detuning-dependent optical depth 1 . By fitting a Gaussian to b(x, y), we extract the size of the cloud R (RMS radius), the on-resonance optical thickness b 0 (from the amplitude of the Gaussian fit) and the number of atoms N (from the integral of the Gaussian). Moreover, the temperature of the cloud is determined by the well-known time-of-flight (TOF) technique.

Note that the same laser beam is used for the excitation of the cloud by the 12 pulses and for the absorption imaging, however the absorption imaging is always performed at low saturation parameter and large detuning ∆ = -4Γ. In the same way, the absorption imaging is also done for different times of flight, without and with a few applied probe pulses prior to it, in order to measure the initial temperature of the cloud, as well as to extract the heating induced by the pulses (see Fig. I.4).

Now

Note that the dynamical control of experimental cycles, including the dynamical exchange between fluorescence measurement and absorption imaging cycles, is accomplished with the MATLAB Graphical User Interface (GUI) via a National Instrument (NI) card, which allows for the communication between the software and the experimental devices such as AOMs, camera, etc.

I.1.3. Calibration of the saturation parameter of the probe laser

Since the objective of this experiment is to characterize subradiance as we increase the saturation parameter of the probe beam s(∆) beyond the linear optics regime, it is important 1 Note that in previous experiments [Guerin 2016, Weiss 2018, Weiss 2019] the detuning-dependent optical depth was defined as b(∆) = g b 0 1 + 4(∆/Γ) 2 , with g = 7/15 the average Clebsch-Gordan coefficient of the transition for the statistical mixture of the Zeeman substates. Here we assume that this g factor is included in the on-resonance optical depth b 0 [Eq. (I.1)], so the difference between the optical depths of this experiment and the previous ones is b 0 /b 0 ≈ 2. With taking into account this, if we compare the obtained subradiant lifetimes in this experiment (for low saturation parameter and most-delayed fit window; see Sec. I.2) with those in the previous experiments for a given b 0 , we can see that the measured subradiant lifetimes are in a very good agreement. This confirms the reproducibility of results of the subradiant experiment over the course of several years.

to have a precise calibration of s(∆). The saturation parameter is defined as

s(∆) = g s 0 1 + 4(∆/Γ) 2 , (I.2)
where g = 7/15 2 is the degeneracy factor of the |F = 2 → |F = 3 transition of 87 Rb (D 2 line) for a statistical mixture of equally populated Zeeman sublevels and s 0 = I/I sat is the on-resonance saturation parameter, with I sat = 1.6 mW/cm 2 the saturation intensity. An initial assessment of the saturation parameter at the center of the probe beam was obtained from the measurement of the total power of the probe beam, P tot before entering the vacuum chamber and the measurement of the probe beam waist (1/e 2 radius w = 5.3 mm). These measurements allowed us to determine the peak intensity (i.e. on-axis intensity) of the probe beam (with a Gaussian profile), which is related to the total power by

I peak = 2P tot /πw 2 .
Note that the cloud size is substantially smaller than the beam radius (R ≈ 1 mm) and the probe beam is aligned such that the cold atomic cloud is at the center of the beam.

Consequently, the probe power is approximately uniformly distributed over the cloud and the intensity seen by the cloud corresponds to the peak intensity to a great degree. Hence, we can determine the saturation parameter with s 0 = I peak /I 0 . However, this is usually not a precise measurement of s(∆) because of a number of effects: losses along the beam path, beam not perfectly Gaussian, atomic cloud not perfectly at the center of the beam, etc. Hence, a calibration method based on the interaction with the atoms is preferable. Hereafter we use the label s (∆) for the saturation parameter that is determined by measuring the total power of the probe beam, while for the properly calibrated saturation parameter we use s(∆). Note that the power of the probe beam is measured simultaneously to the data acquisition by a dedicated photodetector. The other two calibration methods rely on hyperfine depumping into the F = 1 ground state.

Although the transition of interest in this experiment is F = 2 → F = 3, when the probe

2 g = 1 2F + 1 m F |C F,m F →F ,m F | 2 , where C F,m F →F ,m F are Clebsch-Gordan coefficients of the transition of interest.
beam is largely detuned to the red from that transition, there is a non-negligible probability of exciting the F = 2 state from which atoms can decay into the dark F = 1 state. The hyperfine depumping rate into F = 1 state is given by

Γ depump = p 21 Γ 2
s 22 s 22 + 1 , (I.3a) s 22 = g 22 s 0 1 + 4(∆ 22 /Γ) 2 , (I.
(∆) = b 0 /[1 + 4(∆/Γ) 2 ],
which allowed us to precisely determine the resonance frequency.

I.1.4. Spurious effects due to the high saturation parameter: Heating and pushing of the atomic cloud

As mentioned, in this experiment we probe the cold 87 Rb cloud with a beam with a saturation parameter up to relatively large values. The interaction with the light can thus have significant effects on the atomic cloud, in particular heating and pushing.

As described in Sec. I.1.2, the optical depth of the cloud and its temperature are measured simultaneously to the data acquisition using interlaced cycles of subradiance measurements and absorption imaging. From the TOF measurements of the atomic cloud without and with a few applied pulses, we were able to determine the initial temperature of the cloud, as well as extrapolate the temperature of the cloud after each of the 12 applied pulses. While the minimal temperature of the cloud that we measured is T ≈ 100 µK, the maximum temperature for the highest s(∆), considering the heating, is T ≈ 700 µK. As shown in [Weiss 2019], in this range of temperature of laser-cooled atoms, subradiance is not significantly affected: only a very slight decrease of the subradiant lifetime can be expected for the maximum considered temperature. Hence, we do not expect to observe significant effect of the heating on subradiance.

Another effect that could be relevant is the radiation-pressure force exerted by the probe beam on the cloud, inducing a velocity along the beam direction and correspondingly a Doppler shift, which changes the detuning seen by the atoms and thus the effective saturation parameter. This pushing effect can easily be computed and we have checked that, even for the highest saturation parameter, it only induces a very small reduction of the saturation parameter, which does not affect our results. This is confirmed by the fact that different Finally, the expansion of the cloud during the pulse series is also responsible for a slight decrease of the effective intensity interacting with the cloud (due to the finite beam waist), of at most ∼ 20%. It does not affect significantly any of the presented results. As previously discussed, our goal is to characterize the long-lived, i.e. subradiant decay dynamics as we increase the saturation parameter of the probe beam. We increased the saturation parameter by increasing the power of the probe beam from s(∆) ≈ 3 × 10 -3 up to s(∆) ≈ 0.2, 0.8, 1.8 for ∆ = -2Γ, -4Γ, -8Γ, respectively. Note that the larger the detuning, the smaller the maximum obtainable s(∆), as a consequence of the limitation in the available power of the probe. For each value of s(∆) and a given probe detuning, the data acquisition and simultaneous characterization of the cold atomic cloud were done as discussed in the previous section.

In The left panels of these figures show the decay curves normalized with the steady-state intensity of emitted light, i.e. steady-state photon count rate (the steady-state fluorescence level is reached during the 5µs of probe pulse duration). On the right panels are shown the same but non-normalized decay curves corresponding to the number of photon counts per second and per experimental cycle as a function of time after the probe laser switch-off.

Note that t = 0 corresponds to the instant when the probe laser is turned off, i.e. the end of the pulse duration. Here, our focus is on the long-lived tail of the decay curves, i.e. the subradiant part of the fluorescence decay which exhibits slower decay than the single-atom (represented by the black dashed line). When it comes to the normalized decay curves, we observe that the level of the slow (subradiant) decay first increases and then decreases as the saturation parameter becomes larger. On the other hand, the subradiant level of the non-normalized decay curves first increases as the saturation parameter increases and then for s(∆) 0.1, the subradiant level remains constant. These observations hold independently of b 0 and ∆, and they are, with additional quantitative data, further discussed in the following subsection.

I.2.2. Subradiant lifetime and relative amplitude

To characterize the long-lived decay dynamics, we use the following analysis method. We fit the decay curves in a given fit interval by a single exponential function I(t)/I 0 = Ãsub exp(-t/τ sub ) (normalized decay curves), where Ãsub is the relative subradiant amplitude and τ sub is the subradiant lifetime. Similarly, non-normalized curves are fitted by I(t) = A sub exp(-t/τ sub ).

To avoid confusion, we are going to refer to A sub of non-normalized decay curves as the subradiant amplitude, while to Ãsub of normalized fluorescence decay as the relative subradiant amplitude. This decay analysis, i.e. fit by an exponential, was done over several different time intervals Γt Fit at late times where the decay is subradiant. Note that the choice of the latest, i.e. most-delayed fit window Γt Fit is limited by the noise level. Furthermore, it should be noted that it was taken into account that the upper limit of the fit interval is conditioned by the noise level, i.e. it was made sure that the fluorescence corresponding to the upper limit of the fit window is above noise level. Moreover, the fit results for which the goodness of the fit is R 2 < 0.80 are discarded.

In this subsection, our focus is on the normalized fluorescence decay. Subradiant lifetimes and corresponding relative amplitudes as a function of the saturation parameter of the probe s(δ) for several different b 0 and detunings ∆ = -2Γ, -4Γ, -8Γ, as well as different fit intervals Γt Fit are shown in Figures I.8, I.9 and I.10, respectively. As we have investigated the behaviour of the subradiant modes with different lifetimes by monitoring the emission decay over different time windows Γt Fit , several observations can be made. First of all, independently of b 0 and ∆, it can be seen that for the longest-lived observable subradiant modes, with the lifetimes obtained from the latest, i.e. most-delayed fit window, τ sub is not significanlty affected by the increase of s(∆): for all s(∆), subradiant lifetimes are approximately the same. However, for earlier time fit windows, the behaviour of τ sub is notably affected by increasing s(∆): for the several lowest values of s(∆), τ sub is approximately the same, then as s(∆) increases furthermore, τ sub decreases and finally subradiant lifetimes again become independent of s(∆). Moreover, it can be seen that the relative subradiant amplitude Ãsub first increases with increasing s(∆), and then decreases as s(∆) becomes larger (no matter of fit interval Γt Fit ). Note that these observations hold for all b 0 and ∆ of the experiment. It should also be noted that for a given b 0 and a time fit window Γt Fit , the obtained subradiant lifetimes are the same for all ∆ = -2Γ, -4Γ, -8Γ. For example, this can be clearly seen on the inset of Fig. I.11(b), where we directly compare subradiant lifetimes for different ∆, and on Fig. I.11(a), which demonstrates that the temporal decay dynamics is approximately the same for different detunings.

Furthermore, we have checked that the measured steady-state fluorescence level indeed obeys Γ st ∝ s/(s + 1) as expected, and hence in the limit of small s 1, the steady-state level is Γ st ∝ s. An example of such scaling of the measured steady-state fluorescence level is shown on the inset of Fig. I.11(a)). Now let us recall that the Ãsub obtained in this subsection is the subradiant amplitude relative to the steady-state fluorescence level, i.e. it is essentially the subradiant amplitude A sub normalized with Γ st ∝ s/(s + 1). Therefore, for the linear-optics regime, one expects that the relative subradiant amplitude remains constant while s(∆) is varied. However, for the range of s(∆) explored in this experiment, we do not observe such linear-optics behaviour for the subradiant modes. As mentioned, even for the lowest s(∆) of the experiment, we still observe an increase of Ãsub with s(∆), which indicates that for the long-lived modes we are in the regime beyond linear-optics, even for a saturation parameter as low as s(∆) ∼ 10 -3 . This observation will be furthermore discussed below. Now let us turn to the population of subradiant modes P sub . We obtain the population of subradiant modes (normalized with the number of atoms of our cold atomic samples N ) from:

I.2.3. Population of subradiant states

P sub = 1 N ∞ 0 dtA sub,tot exp(-t/τ sub ) = A sub,tot τ sub N , (I.4)
where A sub,tot = ζA sub is the subradiant amplitude associated with the photon scattering rate over all spatial directions i.e., over the full solid angle Ω = 4π. To obtain the correction factor ζ, we determined our detection solid angle Ω det ∼ 4 × 10 -2 sr and under the assumption of isotropic emission we were able to obtain the ratio of detected photons to the number of photons emitted over the full solid angle Ω det /Ω. Moreover, to get even a better estimation of the actual number of emitted photons, we also have to take into account the collection efficiency of the detector η ∼ 3%, which had been independently calibrated in the past.

Hence, the correction factor for the subradiant amplitude that accounts for the scattering over the full solid angle and the detector efficiency is ζ = 1 η Ω Ω det . Therefore, the population of subradiant modes P sub given by Eq. (I.4) represents the number of excitations in these subradiant modes normalized with the number of atoms.

On (see the inset of Fig. I.12(a)). Firstly, we observe that the measured population undergoes a 200-fold increase, from 3 × 10 -7 to 7 × 10 -5 , as the saturation parameter is increased from s ≈ 3 × 10 -3 to 0.3. This corresponds to a maximum number of ≈ 4 × 10 5 excitations in these long-lived modes (for N ≈ 6 × 10 9 atoms). Note that these numbers are only orders of magnitude since the detection efficiency (solid angle, quantum efficiency of the detector, various losses on the optical path) is not precisely calibrated. More importantly, we observe a super-linear increase of the subradiant population: we fit the increasing part of P sub (for s 0.04) by a power-law function P sub ∝ s β which yields β = 1.49 ± 0.16.

We have also checked that the population of long-lived modes with shorter lifetimes (obtained The black dotted line corresponds to P sub ∝ s 1.5 , which illustrates the super-linear growth of the population with the power-law exponential ∼ 1.5. The inset shows the associated subradiant lifetimes for ∆ = -2. Note that for ∆ = -4Γ similar lifetimes are obtained.

for

I.3. Interpretation and comparison with theory

The numerical computations presented in this section have been performed by the group of Romain Bachelard.

I.3.1. Theory

The atomic cloud is modelled as an ensemble of N two-level motionless emitters with positions r j and a transition frequency ω a = kc = 2πc/λ between their ground |g and excited |e states (with σ - m = |g m e m | and σ + m = |e m g m | the lowering and rising operators, respectively), while the natural linewidth of the transition is Γ. The cloud is driven by a laser beam, i.e. near-resonance monochromatic electric field, with Rabi frequency Ω(r), detuned from the transition by ∆. Within the Markov and rotating-wave approximations, the coupled dynamics of the evolution of the density matrix ρ describing the atomic dipoles is obtained from the master equation:

ρ = -i[H, ρ] + L(ρ), (I.5)
where the coherent Hamiltonian Ĥ and dissipative dynamics L are given, in the pump frame, by [Stephen 1964, Lehmberg 1970a, Friedberg 1973]:

Ĥ = -∆ m σ+ m σ- m + 1 2 m Ω(r m )σ + m + h.c. + 1 2 m,n =m ∆ mn σ+ m σ- n , (I.6a) L(ρ) = 1 2 m,n Γ nm 2σ - m ρσ + n -σ+ n σ- m ρ -ρσ + n σ- m . (I.6b)
The diagonal term corresponds to the dynamics of an individual atom, Γ nn = Γ = 1 and ∆ nn = 0, while the coupling terms are given by ∆ mn = -Γ cos(kr mn )/(kr mn ) and Γ mn = Γ sin(kr mn )/(kr mn ), with r mn = |r mr n |. This model corresponds to a 'scalar light' approximation: in the experiment the density of the obtained cloud is in dilute limit (atomic density ρ 0 ≈ 0.06λ -3 ), where a scalar description of the light is a good approximation (see Chapter 2), hence a full vectorial model is not necessary to study the decay dynamics.

The temporal dynamics of the system can be numerically studied by solving the full master equation [Eq. (I.5)], however, since the Hilbert space grows exponentially with the number of atoms as 2 N , only up to a dozen of atoms can be simulated. Note that we are going to refer to simulations from the full master equation [Eq. (I.5)] as exact simulations.

To simulate the dynamics of a larger system (for more than only a dozen saturated two-level atoms), additional approximations in describing the system state are required. Here we use a truncated scheme based on the BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon) approach,

where the density matrix is recast as a sum of reduced density matrices of order m = 1..N , thus establishing a hierarchy of quantum correlations [Bonitz 2016]. The truncation of the hierarchy to two-particle (pair) quantum correlations has proven to be an efficient technique to simulate the dynamics of strongly driven atomic clouds (see [Krämer 2015, Pucci 2017] for more detailed description of the truncated model). The simulations based on this truncated approach, to which we are going to refer as QPC (quantum pair correlations) simulations, allow us to simulate the dynamics of a systems containing up to N ∼ 100 atoms.

Let us also mention the semi-classical, i.e. mean-field (MF) approach3 where all the correlations are neglected (a truncation at the first order) [Krämer 2015]. Note that, while the

MF model assumes σ ±,z m σ ±,z n ≈ σ ±,z m σ ±,z n
, in the QPC approach the pair correlations do

not vanish σ ±,z m σ ±,z n -σ ±,z m σ ±,z n = 0.
Moreover, it should be noted that QPC and MF simulations have been benchmarked with exact simulations: while the mean-field is in a good agreement with the exact simulations only in the case of low saturation parameter (linear optics), the QPC is in a good agreement with the exact simulations for both the low and high saturation parameter regimes.

I.3.2. Experimental vs. numerical results

Numerically, first the steady-state of the system has been computed by driving the atomic system by a probe laser beam (plane-wave) with a saturation parameter s(∆). Then, the laser was switched-off and the decay of the radiated intensity in the far-field limit was obtained from I(n, t) ∝ N j,m e ik n•(r j -rm) σ+ j σm , where n is a unit vector in the direction of observation. Note that the atomic sample is considered to have a Gaussian density distribution, as in the experiment.

As in the analysis of experimental results, the lifetimes and corresponding amplitudes of long-lived modes have been obtained from the exponential fit of decay curves at late times. The red dash-dotted line stands for a fit P sub ∝ s β , with β = 1.007 ± 0.010 (fit on numerical data), indicating the linear-optics (LO) regime for subradiant states. The blue dotted line corresponds to a power-law fit P sub ∝ s β with β = 1.49 ± 0.23 (fit on experimental data), which indicates the regime where the population of subradiant states is being enhanced by a process similar to optical pumping (OP) via superradiant states. Insets: Corresponding lifetimes of long-lived states as a function of the saturation parameter obtained from the experiment (top panel) and from the simulations (bottom panel). (b) The decay dynamics of the population P n of the n-excitation states for a cloud of N = 9 atoms obtained with exact simulations. The atomic cloud is initially driven to steady-state by a strong resonant electric field with Ω 0 = 5Γ. We observe that the lower the excitation number n, the slower the decay. The numerical results presented here were obtained by N. A. Moreira and T.S. do Espirito Santo.

the vastly different parameters for the atom number and size of the sample: N = 100, kR ∼ 6 in the simulations and N ≈ 6 × 10 9 , kR ∼ 8000 in the experiment), it can be seen that the numerical results are in very good agreement with the experimental results: the experimental observation of the super-linear increase of subradiant population as P sub ∝ s 1.49 is also obtained with QPC simulations in the intermediate s regime 2 × 10 -3 s(∆) 4 × 10 -2 .

Additionally, the lifetime of long-lived modes is only marginally affected by the strength of the drive in both, experimental and numerical, cases. While the simulations of the truncated dynamics present an increase of ∼ 15% in lifetime as the saturation parameter is increased, the experimental error bars do not permit us to identify this increase. Moreover, as it is the case for experimental results, the numerically obtained subradiant populations nor corresponding lifetimes do not depend on the detuning of the probe laser. Such observation of independence of subradiant lifetimes on detuning of the probe has also been made in Ref. [Guerin 2016] and it is discussed in Ref. [Guerin 2017b] for the linear-optics regime.

Furthermore, in the simulations it was possible to study the decay dynamics for much lower values of s(∆) than it was achievable experimentally. Hence, in the numerical results the linear optics regime for subradiant modes has also been observed for s(∆) < 2 × 10 -3 , where the subradiant population scales linearly with s: P sub ∝ s. The simulations at b 0 = 5 present a non-linearity threshold at s LO ≈ 2 × 10 -3 . It was recently suggested that such a threshold may scale as Γ 2.5 n , with Γ n = 1/τ n the n-th mode linewidth, i.e. decay rate [Williamson 2020]: assuming that the subradiant states present linewidths scaling as Γ n ∼ Γ/b 0 [Guerin 2016], saturation parameters orders of magnitude smaller may be necessary to experimentally reach the linear-optics regime for long-lived states (in the experiment we probed the clouds with larger on-resonance optical depths and consequently the long-lived (single-photon) modes have longer lifetimes than what is the case for b 0 = 5 as in the simulations). It should be noted that numerically it was not possible to obtain a scaling of s LO from the low-b 0 simulations.

The results obtained here show that in the previous experiments [Guerin 2016, Weiss 2018, Weiss 2019], the data was actually not taken in the linear-optics regime for collective long-lived modes with s(∆) ∼ 10 -2 . This is in contradiction with what was claimed in the Supplemental Material of Ref. [Guerin 2016]: no significant dependence of the relative subradiant amplitude on the saturation parameter in the range of s(∆) explored there had been observed, indicating that the experimental data for s(∆) ∼ 10 -2 , had been taken in the linear-optics regime for subradiant modes (the relative subradiant amplitude is expected to be unaffected by the probe intensity in the linear-optics regime). However, as the range of s(∆) explored there was relatively short, the increase of relative subradiant amplitude with s(∆) was not obvious and the wrong conclusion of being in the linear optics regime for s(∆) ∼ 10 -2 had been made.

Note that the exact simulations with N = 7 atoms show a very good agreement with the QPC simulations: for the very low s, the population scaling P ∝ s is observed indicating the linear-optics regime for long-lived modes, while for intermediate s, the superlinear-scaling P sub ∝ s β , with β ∼ 1.5 is identified. Moreover, the semiclassical simulations, i.e. the mean-field approach which neglects correlations, have been benchmarked with the exact (and QPC) simulations and it has been observed that the mean-field model fails to properly describe the decay of subradiant states in the high-s regime, i.e. the mean-field model does not capture the long-lived decay dynamics beyond linear-optics regime. This strongly suggests that the subradiant states here studied are due to quantum correlations [Tana 2004].

I.3.3. Interpretation N=2-atom case

To interpret our results, in particularly the observed super-linear scaling of subradiant population, we consider a "toy-model" of N = 2 two-level atoms with transition frequency The steady-state population of the long-lived mode then presents three typical regimes, depending on the pump strength. First, for the lowest intensities (linear-optics regime), the population of |ee is negligible and the single-excitation modes |± are driven only directly from the pump, so one obtains the following scaling for their population:

ω a = kc = 2πc/λ
P + -≈ s + -/2 ∝ s, with s + -= 2Ω 2 + -/(Γ 2 + -+ 4(∆ ∓ ∆ + -/2) 2
) the effective saturation parameter for each mode, and for s ± 1, i.e., Ω + ∆.

s = 2Ω 2 0 /(Γ 2 + 4∆ 2 )
As the drive strength is increased, the doubly-excited state |ee is substantially populated thanks to the strong coupling of the drive to the superradiant state: P + ≈ s + /2 ∝ s and

P ee ∝ s 2 .
Then, the |-state gets an additional population by decay from |ee , at rate Γ -, leading to a long-lived population that grows quadratically with the saturation parameter:

P -∝ s 2 [see Fig. I.16(c)].
Finally, for the largest values of the saturation parameter, i.e., with a Rabi frequency such that the dynamics of each atom is dominated by the drive (Ω 0 ∆ Γ, |∆ ± |), the system is cast into a separable state described by the density matrix ρ = ⊗ j=1,2 (|g

j g j | + |e j e j |)/2.
This mixed state projects equally on the states |gg , |+ , |-and |ee , resulting in P -≈ 1/4.

Hence, the strong pump overcomes the weak coupling of subradiant modes which, in the linear-optics regime, prevents one to populate them efficiently. The present mechanism is analogous to optical pumping, where an excited state (here |ee ) is directly driven by the laser, and induces a population in the long-lived state (here |-) by incoherent decay.

Our N 2 atomic sample

The discussed interpretation of enhanced subradiant population for N = 2 atoms in terms of a process similar to optical pumping can be extended to our results obtained with a cold Thus, it is expected that the late time decay dynamics is going to be dominated by the single-excitation long-lived modes as they exhibit the slowest decay. This explains why we experimentally observe that the most delayed fit window yields the longest observable lifetimes that are unaffected by the increasing saturation parameter (see for example Fig. I.12):those longest lifetimes correspond to the single-excitation subradiant modes (whose population, as previously discussed, is enhanced by the optical-pumping-like process). Moreover, the fact that the subradiant lifetimes obtained from the fit over the earlier time windows decrease with s for higher saturation parameters possibly indicates that this shorter lifetimes correspond to the decay from the multi-excitation long-lived modes.

atomic

I.4. Conclusion

By studying the slow decay dynamics of cold atomic clouds as the saturation parameter s of the laser drive is increased, we have experimentally observed that the population of long-lived, i.e. subradiant modes exhibits the superlinear growth with s as P sub ∝ s β , with the exponent β ∼ 1.5. We attribute this beyond linear-optics behaviour for long-lived modes to a process similar to optical pumping through decay processes from multi-excitation superradiant states. Moreover, the subradiant lifetimes observed at the latest times, which correspond to a single-excitation subspace, remain unaffected by the strength of the laser drive.

The numerical simulations (QPC and exact) show a very good agreement with experimental results. However, the semiclassical simulations (which neglect the quantum correlations) have failed to reproduce the results, which is an argument to support the idea that the subradiant states might be appropriate to store entanglement or quantum correlations.

The fact that we have experimentally demonstrated that the population of long-lived collective modes can be enhanced by increasing the intensity of the laser drive beyond the linear-optics regime, while the corresponding lifetimes remain unaffected, is an encouraging step towards harnessing the subradiant modes for information storage in quantum memory devices.

Chapter II

Subradiance beyond the dilute regime

In previous experimental studies on subradiance, it had been shown that for dilute (such that the typical distance between the atoms is larger than the wavelength) and large (where the size of the sample is much larger than the wavelength) cold-atom samples, as well as in the linear-optics (weak-excitation) regime, the collective spontaneous emission is governed by the resonant optical thickness b 0 of the sample, i.e., subradiant lifetimes scale linearly with b 0 [Guerin 2016, Weiss 2018]. From the numerical side, the scalar-coupled dipole model, which involves only the long-range dipole-dipole interaction (∝ 1/r), captured very well those experimental observations [Guerin 2016, Araújo 2018b].

The question is how the near-field dipole-dipole interaction affects the long-lived subradiant modes, since it is supposed to become non-negligible beyond the dilute limit, as the distance between atoms becomes comparable to the wavelength. The near-field effects have been reported for several light-scattering phenomena. For example, in the case of the Dicke limit of a small sample (whose size is much smaller than the wavelength) for initially fully inverted systems (each atom is excited), it was reported that the near-field dipole-dipole interaction breaks the symmetry of the collective states yielding a reduced superradiance [Friedberg 1972, Friedberg 1973, Friedberg 1974, Gross 1982]. Moreover, it has been demonstrated that in dense samples, the near-field interaction is at the origin of a suppression of light scattering Hence, the object of the study presented in this chapter is to explore what effect the near-field interaction has on subradiant modes as we increase the density of the sample. The study was done numerically with the vectorial coupled-dipole model (and the scalar model to compare the vectorial results to the scalar ones). Note that we considered large samples (whose size is larger than the wavelength) and the linear-optics regime.

The chapter is organized as follows. In Sec. II.1, we introduce the vectorial and scalar coupled-dipole models, as well as the atomic sample. In Sec. II.2, we present the results obtained with the vectorial and scalar models, in particular how subradiant lifetimes depend on the density of atomic samples. We show that the near-field dipole-dipole interaction is detrimental for subradiance, which we see in the reduction of subradiant lifetimes with increasing density. In addition, we extensively studied the influence of subradiant pairs and positional correlations on the slow decay dynamics to out rule those effects as being at the origin of the reduction of subradiance. We also discuss the eigenvalue distribution, which allows us to interpret the detrimental character of near-field effects on subradiance. Finally, in Sec. II.3, we summarize and discuss the results.

II.1. Microscopic model

II.1.1. Coupled-dipole equations

We consider a system of N identical 4-level atoms with ground state |J g = 0, m g = 0 coupled to triple-degenerate excited states |J e = 1, m e = 0, ±1 by the electric field, with a coupling provided by the dipole transition moment. Each atom is treated as a discrete point-like radiating dipole located at fixed position r j , where j = 1, ..N , and with transition frequency ω 0 = ck 0 , where k 0 is the amplitude of the wavevector associated to the atomic dipole transition. The incident electric field is a plane-wave laser beam described by E in = E L exp(ik L • r j -iω L t), characterized by its amplitude E L , polarization ˆ L (with

E L = E L ˆ L ), wavevector k L = k L ẑ (k L ≈ ω 0 /c =
k 0 and ẑ used as a quantization axis), and frequency ω L detuned by ∆ 0 = ω L -ω 0 from the atomic transition. We here use the spherical basis, with unit vectors ê±1 = ∓1/ √ 2(ê x ± iê y ), and ê0 = êz . Throughout this work we use right-hand circular polarization for the laser beam: ˆ L = ê-1 .

Our focus is here on a weak driving field, such that the system presents a linear response to the field, i.e., the linear optics regime. By using rotating-wave approximation (in the atom-drive interaction Hamiltonian) and Markov approximation, the optical response of the system is given by a set of following 3N coupled-dipole equations (CDEs) (for example, see

Refs. [Lehmberg 1970a, Lehmberg 1970b, Manassah 2012, Samoylova 2014, Bienaimé 2011b] for the derivation of these equations):

dβ ζ j (t) dt = i∆ 0 - Γ 0 2 β ζ j (t) -i d ê * ζ • E L exp(ik L • r j ) - Γ 0 2 m =j η G ζ,η (r jm )β η m (t), (II.1) with j, m ∈ 1, N , ζ, η ∈ (±1, 0
) the spherical-basis components, d the electric-dipole transition matrix element (considered to be the same for every atom and transition), Γ 0 = d 2 k 3 0 /3 π 0 the single-atom decay rate and r jm = r jm rjm = r jr m the distance vector between atoms j and m. Note that these equations for β ζ j are written in the rotating frame of the driving field. Moreover, this set of 3N coupled first-order linear differential equations describes the temporal evolution of the amplitudes of the atomic dipoles β ζ j (t). Note that

G ζ,η (r) = 3 2 exp(ik 0 r) ik 0 r δ ζ,η -rζ r * η + δ ζ,η -3r ζ r * η i k 0 r - 1 (k 0 r) 2 , (II.
G(r) = 3 2 exp(ik 0 r) ik 0 r 1 + i k 0 r - 1 (k 0 r) 2 ↔ I + -1 - 3i k 0 r + 3 (k 0 r) 2 r ⊗ r r 2 . (II.3)
As already mentioned, the dipole components β ζ j represent the amplitude of the induced oscillating atomic dipole, and therefore the components of the scattered electric field vector at position R = Rn can be obtained from the emission of the dipoles located at r j :

E ζ sc (R, t) = -i dk 3 0 6π 0 j η G ζ,η (R -r j )β η j (t). (II.4)
In the far-field limit (R

r j , 1/k 0 ), we can approximate |R -r j | ≈ R -n • r j and 1/|R -r j | ≈ 1/R in G ζ,η (R -r j ) [Eq. (II.2)], while we drop terms ∝ 1/|R -r j | 2 and ∝ 1/|R -r j | 3
. Under these approximations, the expression of the vectorial Green's function

reduces to G ζ,η (R -r j ) ≈ 3 2 exp(ik 0 R) ik 0 R (δ ζ,η -nζ n * η ) exp(-ik 0 n • r j )
, where n = R/R is the unit vector in the observation direction. Therefore the components of the scattered electric field [Eq. (II.4)] in the far-field limit are approximated by

E ζ sc (R, t) ≈ -i dk 3 0 4π 0 j η exp(ik 0 R) ik 0 R (δ ζ,η -nζ n * η ) exp(-ik 0 n • r j )β η j (t). (II.5)
Since the electric field in the far-field limit is purely transversal, the intensity of the scattered light in the far field can be computed as

I sc ∝ |E sc | 2 = ζ |E ζ sc | 2 : I sc (n, t) ∝ m,j e -ik 0 n•r jm ζ,η (δ ζ,η -nζ n * η )β η j (t)β ζ * m (t). (II.6)
In the above expression for the intensity of the scattered light, we did not write factor ∝ 1/R 2 before the summations. As we will normalize the obtained scattered intensity, the prefactor is not essential; in our computations of the scattered light in the far-field limit, only the observation direction n determined by the observation angles {θ sc , φ sc } is relevant.

The scalar approximation of the coupled-dipole model, which disregards the vectorial nature of light (i.e., its polarization) and the internal Zeeman structure of the atoms, is obtained by averaging G ζ,η (r jm ) in Eq. (II.2) over random orientations of the pairs of atoms j and m. In the case ζ = η, one obtains rζ r * η = 1/3, while rζ r * η = 0 for ζ = η. Consequently, the near-field terms disappear, as well as polarization coupling, and we obtain the following scalar dyadic Green's function:

G (s) (r) = exp(ik 0 r) ik 0 r . (II.7)
The atomic dipoles are then described by a scalar amplitude β j , whose dynamics is given by the set of N scalar CDEs:

dβ j (t) dt =   i∆ 0 - Γ (s) 0 2   β j (t) - dE L e ik L •r j - Γ (s) 0 2 m =j G (s) (r jm )β m (t) , (II.8)
where the natural decay rate differs by a factor 3/2 from the vectorial one: Γ (s) 0 = (3/2)Γ 0 . The dipole matrix element and the scattering cross section also differ by a factor from the scalar to the vectorial case. The difference between the vectorial and scalar definitions of the natural decay rate stems from the absence of polarization degrees of freedom of light in the scalar model. 1 Note that, in the scalar model, the coupling of atomic dipoles is provided only by the long-range (far-field; ∝ 1/r) dipole-dipole interaction term.

The scattered electric field in the scalar approximation reduces to

E (s) sc (R, t) = -i dk 3 4π 0 j G (s) (R -r j )β j (t), (II.9)
and the scattered intensity can be obtained as

I (s) sc ∝ |E (s) sc | 2
which in the far-field limit reads:

I (s) sc (n, t) ∝ m,j e -ik 0 n•r jm β j (t)β * m (t).
(II.10)

Let us now summarize and add some more details about steps needed to obtain the intensity of the scattered light from CDEs to study the decay dynamics. We have to solve the set of 3N

vectorial CDEs [Eq. (II.1)] or the set of N scalar CDEs [Eq. (II.8)] for β(t). For convenience, we will explain the following steps in the case of the vectorial model, but note that those steps are analogous in the scalar model. First we compute the steady-state solution,

dβ ζ j,st
dt = 0 of Eq. (II.1) by solving the following linear problem:

β ζ j,st = 1 i∆ 0 -Γ 0 /2   i d ê * ζ • E L exp(ik L • r j ) + Γ 0 2 m =j η G ζ,η (r jm )β η m,st   ∀ (j, ζ). (II.11)
Note that by combining this equation and Eq. (II.4), we can rewrite the steady-state CDEs as

β ζ j,st = d (∆ 0 + iΓ 0 /2) ê * ζ • E L exp(ik L • r j ) + E ζ (r j ) , (II.12)
from which we can clearly understand that the amplitude of each atomic dipole j is not only the result of driving by the incident electric field but it is also affected by the electric field

E ζ (r j ) = -i dk 3 0 6π 0 m =j η G ζ,η (r j -r m )β η m (t)
, produced by all the other N -1 dipoles.

After we obtain the stationary amplitudes of atomic dipoles β ζ j,st , we switch-off the driving field (E L = 0) in Eq. (II.1), which gives following set of equations 1 In our vectorial [Eq. II.1] and scalar [Eq. II.8] coupled-dipole simulations, we set Γ 0 = 1 and Γ (s) 0 = 1, respectively. Like that the temporal dynamics for dipoles in vectorial and scalar models is in the units of Γ 0 and Γ (s) 0 , respectively.

dβ ζ j dt = i∆ 0 - Γ 0 2 β ζ j - Γ 0 2 m =j η G ζ,η (r jm )β η m . (II.13)
Then we solve those equations for β ζ j (t) after the switch-off of the incident driving field with the previously obtained steady-state solution as an initial condition β ζ j,st -→ β ζ j (t = 0). In the end, we use the obtained dipole amplitudes β ζ j (t) to compute the scattered intensity given by Eq. (II.6).

In a described way, we compute the scattered intensity for various spatial configurations of atomic positions, i.e., for many realizations of discrete atomic positions. Then we average the intensity over those realizations. Here we choose the number of realizations N r such that the product of N r and the number of atoms N is always the same: N r × N = 60000. As is going to be discussed below, we use a spherical atomic sample, i.e. N atoms confined in a spherical volume. Thus, considering the azimuthal symmetry of such spherical sample (up to the disorder), the obtained intensity is averaged over 51 azimuthal angles φ sc = [0 • , 360 • ). It is then the average intensity (averaged over many realizations and the azimuthal angles) that we use further for our study of temporal decay dynamics.

II.1.2. Spectrum of the system

To interpret the results obtained from the temporal dynamics of the light scattered by atoms, we will look into the distribution of eigenvalues of the collective modes. Let us first introduce the generalized vectorial 3N × 3N matrix G, whose elements are

G η,ζ j,m =        i∆ 0 - Γ 0 2 for j = m, ζ = η - Γ 0 2 G η,ζ (r jm ) for j = m given by Eq. (II.2), (II.14)
such that the vectorial CDEs [Eq. (II.1)] can be rewritten as coupling terms that describe the dipole-dipole interaction via the electromagnetic field. As a result of the coupling between dipoles, collective modes arise which are characterized by

βζ j = m η G η,ζ j,m β η m -i d ê * ζ • E L exp(ik L • r j ). ( II
eigenstates ψ n = ψ 1,1 n , ψ 1,0 n , ψ 1,-1 n . . . ψ j,ζ n . . . ψ N,1 n , ψ N,0 n , ψ N,-1 n
. By diagonalizing the matrix

G, Gψ n = λ n ψ n , we obtain 3N complex eigenvalues λ n = i∆ n -Γ n /2 corresponding to each eigenmode ψ n .
The dispersive imaginary part is associated to the frequency of the eigenmode ω n relative to the frequency of the driving field ω L , i.e., the imaginary part yields the frequency shift of the eigenmode:

∆ n = ∆ 0 + (ω 0 -ω n ) = ω L -ω n = I(λ n ). (II.16)
The dissipative real part of λ n corresponds to the decay rate of the eigenmode:

Γ n = -2 R(λ n ) (II.17)
which characterizes the coupling of the corresponding eigenmode with the vacuum electromagnetic field. Subradiant states, characterized by Γ n < Γ 0 , are coupled weakly to the vacuum EM field compared to the excited state of the single atom, while superradiant states, Γ n > Γ 0 , are coupled strongly.

Another quantity that will be of interest is the inverse participation ratio (IPR), which tells us how many atoms participate in an eigenmode (inversely): it is a measure of 1/N n where N n is the number of atoms that significantly contribute to the mode ψ n . It is given by

IPR n = j,ζ |ψ j,ζ n | 4 ( j,ζ |ψ j,ζ n | 2 ) 2 . (II.18)
In our focus are modes to which many atoms are significantly participating, i.e., modes with

IPR n 1 (Fig. II.1
). For convenience, those are the modes that hereafter we are going to refer to as collective modes, to distinguish them from pair modes, to which two atoms significantly contribute.

Equivalently, in the scalar model N eigenvalues λ (s) n and eigenstates ψ n = ψ 1 n , . . . ψ j n . . . ψ N n can be obtained from the N × N matrix G (s) with elements given by

G (s) j,m =          i∆ 0 - Γ (s) 0 2 for j = m - Γ (s) 0 2 G (s) (r jm )
for j = m given by Eq. (II.7) ,

(II.19)
and the corresponding IPR in the scalar model is provided by

IPR n = j |ψ j n | 4 ( j |ψ j n| 2 ) 2 .
(II.20)

Close pairs of atoms

As it was already pointed out, we are interested in the decay of collective modes, which involve many atoms (N n 1). However, we can see in As they might significantly influence the decay dynamics and even overshadow the collective dynamics, they will need special treatment: we will impose a minimum separation distance on atoms, i.e., an exclusion volume (see the following subsection), to reduce their influence.

The influence of pairs on collective dynamics is discussed in Sec. II.2.3.

To understand the properties of pairs, let us first take a look at the spectrum of a two-atom system. For vectorial waves, the eigenvalues λ pair = i∆ pair -Γ pair /2 corresponding to modes of the atomic pair can be found by diagonalizing the 6 × 6 vectorial Green's matrix [Eq. (II.14)].

For N = 2, it comes down to two eigenvalues of single multiplicity

λ pair,± = - Γ 0 2 1 ± 3 2i e ik 0 r 12 - 2i (k 0 r 12 ) 2 + 2 (k 0 r 12 ) 3 + i∆ 0 , (II.21)
and two eigenvalues of double multiplicity

λ ⊥ pair,± = - Γ 0 2 1 ± 3 2i e ik 0 r 12 1 k 0 r 12 + i (k 0 r 12 ) 2 - 1 (k 0 r 12 ) 3 + i∆ 0 , (II.22)
where r 12 = |r 1r 2 | is the distance between two atoms. λ pair,+ and λ pair,-stand for eigenvalues corresponding to symmetric (superradiant) and anti-symmetric (subradiant) eigenstates, respectively 2 . Let us assume that the atoms are placed along some axis n, such that r 12 = r 12 n. The problem of two-dipole coupling can be decomposed as follows: dipoles oriented along the n axis and dipoles whose orientation is orthogonal to the n axis (for that, there are two mutually perpendicular possibilities). Hence, the eigenvalues of single multiplicity λ pair correspond to dipoles whose orientation is parallel with the alignment axis n. In contrast, double-multiplicity eigenvalues λ ⊥ pair correspond to the case of dipoles oriented orthogonally to the n axis. Note that for superradiant/subradiant modes of the two-atom system, we will refer to as superradiant/subradiant pairs. In Here, eigenvalues were computed for ∆ 0 = 0, such that ∆ pair = ω 0 -ω pair [see Eq. (II.16)].

Note that the decay rates and frequency shifts of pair modes strongly depend on the short separation; the smaller the distance between two atoms, the larger the frequency shift |∆ pair,-correspond to superradiant and subradiant eigenstates, respectively. In the limit kr 12 1, the cooperativity of pairs disappears, i.e., pair decay rates and frequency shifts become equal to those of a single atom: Γ to the long-range dipole-dipole interaction [Guerin 2016].

Let us remark on the inequality between λ pair and λ ⊥ pair eigenvalues. In Fig. II.2(a,b) we can see that for intermediate distances between the two limits k 0 r 12 1 and k 0 r 12 1, the superradiant/subradiant decay rates and frequency shifts of λ pair and λ ⊥ pair are not the same; Γ pair,± (r 12 ) = Γ ⊥ pair,± (r 12 ) and |∆ pair,± (r 12 )| = |∆ ⊥ pair,± (r 12 )|. For subradiant modes, that is the case even in the k 0 r 12 1 limit, where3 (II.24) This difference between and ⊥ modes is an important property of pairs that we are going to exploit to rule out the influence of subradiant pairs on the obtained collective-dynamics results (Sec. II. 2.3) In Fig. II.2(c) we show Γ pair (k 0 r 12 ) vs. ∆ pair (k 0 r 12 ) curves. We will refer to those curves as super-and subradiant pair branches.

Γ ⊥ pair,-= 2Γ pair,-≈ (k 0 r 12 ) 2 5 Γ 0 ∆ pair,-= -2∆ ⊥ pair,-≈ - 3Γ 0 2(k 0 r 12 ) 3 .

II.1.3. Atomic sample

The atomic sample is modeled by assuming a discrete-particle representation of the sample.

We consider a spherical cloud of N identical, motionless atoms with a Gaussian density distribution ρ G (r j ) = ρ exp(-r 2 j /2R 2 ), where r j is the position of the jth atom, R is the RMS radius of the Gaussian cloud and

ρ = N/( √ 2πR) 3 (II.25)
is the peak density of the cloud, i.e., density at the center of the cloud. Throughout this chapter we will use dimensionless density ρλ 3 , with λ = 2π/k 0 the wavelength of the dipole transition. Furthermore, the peak resonant optical thickness of the Gaussian cloud is

b 0 = σ 0 ρ G (0, 0, z) dz = 3N/(k 0 R) 2 , (II.26)
with σ 0 = 6π/k 2 0 the on-resonance scattering cross section. Note that the given expression for the optical thickness b 0 , as well as for the scattering-cross section σ 0 corresponds to the vectorial light. The scattering cross-section for the scalar light differs by the factor of 2/3 from the vectorial-light scattering cross-section:

σ (s) 0 = 2 3 σ 0 .
Consequently, the scalar-light optical thickness is given by b 2 . The same as with the scalar and vectorial definitions of the single-atom decay rate, the difference between the scalar and vectorial scattering-cross section is caused by the absence of polarization degrees of freedom of light in the scalar model.

(s) 0 = 2 3 b 0 = 2N (k 0 R)
Therefore, for a given N and R, and corresponding b 0 and ρλ 3 , we generate N random positions of atoms to form the spherical cloud with the Gaussian density distribution which is determined by the peak density ρλ 3 and RMS size R. We refer to each such event of creating the random set of atomic positions for the same parameters of the Gaussian spherical cloud, as a realization of the atomic positions.

Our interest lies in characterizing subradiance in atomic clouds of different optical thickness b 0 and density ρλ 3 , quantities that both depend on the number of atoms N and size of the sample k 0 R. However, the number of atoms that we can simulate is limited (up to several thousands). We also impose the condition on the sample size R > λ in order to consider macroscopic samples. As a consequence of these limitations on the number of atoms and the sample size, the range of ρλ 3 for a given b 0 (and vice versa) that can be achieved is limited. Therefore, the dependence of the subradiant lifetime on b 0 and ρλ 3 can only be studied piecewise. The range of density and on-resonant optical thickness that we study here is ρλ 3 = [0.8; 40] and b 0 = [2; 72], respectively.

Note that as one increases the density of the cloud, the probability that close pairs of atoms are generated becomes higher. These pairs, introduced in the previous subsection, result in superradiant and subradiant modes that are characterized by strong energy shifts [Stephen 1964] (see Fig. II.2). When driving the system by a laser field with a significant detuning, some of those pairs may be resonant with the field and consequently be strongly excited. They can then play a significant role in the cloud radiation, despite involving few atoms [Fofanov 2021].

However, these pairs are expected to be highly sensitive to atomic motion since they involve very short distances; therefore, the pairs might not be relevant for experiments with thermal clouds. In this work, our main interest lies in the long-lived collective modes, involving many atoms. As a consequence, we implement a hard-sphere radius for atoms, i.e. an exclusion volume to impose a minimal distance r min between the atoms, r ij = |r ir j | ≥ r min ∀(i, j), thus minimizing the influence of pairs. The method of introducing the exclusion volume r min for the atom positions is the following.

We are adding atoms one by one to the sample, i.e., we generate random positions for one atom at a time. After adding each atom i to the sample, we check if distances between just added atom i and all other atoms j that are already in the sample satisfy the exclusion volume condition r ij ≥ r min ∀ j. If the condition is fulfilled, we keep the atom i, i.e., the atomic position r i . On contrary, if r ij ≤ r min for any atom j, the position r i is discarded and we keep drawing new positions of atom i until the exclusion volume criterion is met.

Furthermore, since adding an exclusion volume can lead to an increased sample size R, after generating the atomic sample, we determine the actual RMS size of the cloud and recalculate b 0 and ρλ 3 .

Imposing the exclusion volume, that is relatively large for a given peak density of the Gaussian atomic sample, might cause the distortion of the Gaussian density distribution. However, we have checked that with our choice of the exclusion volume for all atomic clouds with a given (b 0 , ρλ 3 ) studied in this chapter, the shape of the cloud remains Gaussian to a very good approximation.

The results presented in Section II.2.2 have been obtained with a density-dependent exclusion volume defined as r min = ρ -1/3 /π since it allows us to explore high densities without introducing significant positional correlations, while efficiently removing the pairs. In Section II.2.3, we discuss in more details the possible influence of pairs or positional correlations on the results presented in the following section.

II.2. Late-time decay dynamics II.2.1. Temporal evolution of the scattered intensity and the method of obtaining subradiant lifetimes

In Sec. II.1.2, we discussed that the coupling between dipoles yields collective modes, each characterized by its own decay rate and eigenfrequency. However, not all modes contribute significantly to the decay dynamics. The coupling of the external driving field to the mode determines the population (excitation) of the mode [Guerin 2017b], which dictates how relevant, i.e., dominant the mode will be in the decay process. Hence, the temporal behavior of the scattered light intensity is encoded in the population of the modes. In Let us briefly remark on the role of the detuning of the external driving field on the population of modes [Guerin 2017b]. If the driving field is resonant (or near-resonant) with some of the modes, those modes are more strongly coupled to the driving field and hence more strongly excited than the other modes. We here use relatively large detuning ∆ 0 = ±15Γ 0 such that none of the collective modes are near-resonant with the driving field; all collective modes are driven weakly but equally. Note that this does not mean that all collective modes would be equally populated in the case of the large drive detuning. As was discussed in [Guerin 2017b],

there is a geometrical factor independent of the detuning that corresponds to the projection of the driving field on the modes, and it mainly governs the population of the collective modes when the detuning of the drive is large.

Our primary interest is in the late-time decay dynamics of the atomic ensemble, where the decay observed in the scattered light is slower than that of the single atom: the subradiant decay. In order to study the temporal dynamics of the atomic system after the switch-off of the driving field, we use the CDEs to obtain the intensity of the scattered light as a function of time [Eq. (II.6)], as was discussed in the previous section. Note that the obtained intensity is normalized by its value at the instant of the laser switch-off, I 0 = I(t = 0).

To characterize subradiance in atomic samples of different parameters, we need a reliable analysis method of obtained decay curves. Our focus is in determining the subradiant lifetimes τ sub = 1/Γ sub (where Γ sub is the subradiant decay rate). To obtain the subradiant lifetimes from the computed decay curves, we use a fit procedure, as it was used, for example, in [Guerin 2016, Araújo 2018b]. Normalized decay curves are fitted in a given range at late times by the following single exponential function

I(t)/I(0) = A sub exp(-t/τ sub ), (II.27)
with subradiant amplitude A sub and subradiant lifetime τ sub as free parameters.

To explain the choice of the fit range, let us have a look on is quickly decreasing, and then at later times, there is a prolonged decrease of Γ(t). For those later times, we can also clearly identify the well-known cooperativity enhancement of subradiance with b 0 for macroscopic samples: the larger the resonant optical thickness is, the slower the subradiant decay is [Guerin 2016, Araújo 2018b]. We refer to those times as late times relevant for determining subradiant lifetimes and characterizing subradiance in different macroscopic samples.

Because of the slow evolution of the decay at late times, we can find a fit interval of an appropriate length for which the single-exponential fit is a good approximation of the intensity ). The choice in the length of the particular fit range was verified with the R 2 coefficient of the fit (the goodness of the fit). Moreover, the coefficient of variation of the subradiant lifetimes obtained from the exponential fit throughout this chapter is < 1%.

As it was discussed in [Araújo 2018b], the decay rate becomes constant at very late times for which the relative intensity level is very low, I(t)/I(0) 10 -10 , beyond the reach of any current real experiment. Moreover, the final decay rate and the corresponding intensity level at which it emerges depend on the longest-lived collective mode and its population, which might vary from one realization to another. Therefore, to get reliable temporal dynamics for the last collective decay rate, one would need a tremendous amount of realizations. Because of all of that, we are focusing here on higher intensity levels that capture enhancement of subradiant lifetimes with increasing b 0 , although Γ(t) is still slowly evolving in time.

Lastly, let us discuss why we use fit intervals based on the relative intensity level rather than the time intervals. To have a proper comparison of τ sub for different system parameters, we should use the same fit interval at late times for all decay curves of different parameters. As it can be seen on here, would mean that for different optical depths, the intensity levels corresponding to that time-fit interval would be very different. For the lowest optical depths, the corresponding intensity fit levels would be very low. From the experimental point of view, that would be impractical. As we want our analysis approach to be achievable in potential experiments, we use fit intervals based on the intensity level as it is presented in Fig. II.3(a). We have checked that the choice of the fit interval does not qualitatively alter our results presented in the following subsections. To understand how near-field terms affect subradiance, we monitor the scattered light intensity after the switch-off of the driving field, using the vectorial CDEs [Eq. (II.6)]. In this subsection, we consider the total light scattered in a given direction. Here, total means that we collect all polarization components of the scattered light. Furthermore, here we treat the system that is first driven to a steady-state with a plane-wave laser beam with a large detuning

II.2.2. Scaling of subradiance with vectorial light

∆ 0 = -15Γ 0 .
We present in Fig. II.4 the subradiant lifetimes as a function of b 0 for several ρλ 3 . One can see that for lower densities, the sets of data points (b 0 , τ sub ) corresponding to different ρλ 3 collapse on the same line. This is also the case for the scalar model in the dilute regime [Guerin 2016, Araújo 2018b], and it shows that density effects are negligible for the lowest densities. However, for ρλ 3 5, the data sets do not collapse anymore: higher-density samples present shorter subradiant lifetimes for a given optical thickness. This effect of higher densities is even clearer in Fig. II.5(a), where we present the subradiant lifetimes as a function of the density for several values of b 0 . Again, there is no visible effect of the density on the long-lived emission for the lowest densities, but for ρλ 3 5, the late lifetimes become shorter with increasing densities while the on-resonance optical thickness remains constant.

The effect of reduction of subradiant lifetimes with increasing densities is as well seen directly . There is no threshold k 0 R for which the reduction of τ sub starts to occur for all data sets. On the contrary, in Fig. II.5(a) we identified the density threshold ρλ 3 ≈ 5: the value above which τ sub decreases with increasing ρλ 3 , while below that value τ sub is independent of ρλ 3 for a given b 0 . Since there is the unique range of densities (ρλ 3 5) but not of size for which τ sub decreases, we can conclude that the reduction of τ sub is indeed a density effect and it is not caused by the decrease of the sample size.

We have checked that driving the sample with the opposite-sign detuning (∆ = +15Γ 0 ), as well as with a larger detuning (∆ = -30Γ 0 ), yields the same subradiant lifetimes as those in Figs. II.4 and II.5. This excludes density-induced collective shifts [Manassah 2012[START_REF] Javanainen | [END_REF], Zhu 2016, Jenkins 2016, Jennewein 2018] of the atomic resonance as a source of the observed effect. We have also checked that using other late-time fit intervals (for I(t)/I(0) < 10 -4 ) leads to the same conclusion. Although we obtain slightly different subradiant lifetimes with different fit windows, the qualitative behavior is the same: we observe The yellow curve corresponds to the decay of a single atom. For ρλ 3 5 the decay at late times becomes slightly faster as the density is increased. The inset shows the complete decay curves from the moment the laser is switched-off (t = 0). The intensity of the scattered light was computed from the vectorial CDEs [Eq. (II. [START_REF] Cipris | Subradiance with Saturated Atoms: Population Enhancement of the Long-Lived States A[END_REF]]. Other simulation parameters: ˆ L = ê-1 circular polarization of the driving field, ∆ = -15Γ 0 , θ sc = 45 • , r min = ρ -1/3 /π. the decrease of subradiant lifetimes with increasing density for a given optical thickness.

Furthermore, we have verified with several off-axis observation angles θ sc = 45 • , 90 • , 135 • , 180 • that the conclusion reached from Figs. II.4 and II.5 is independent of the observation angle θ, provided it is outside the forward diffraction lobe, where peculiar effects associated to superradiance may occur [Kuraptsev 2017]. Finally, since for higher densities the scaling of subradiant lifetimes with b 0 fails, we have checked (for data from Figs. II.4 and II.5) that the unique scaling with system parameters N , k 0 R, or N/k 0 R does not emerge.

II.2.3. The exclusion volume

The results shown in the previous subsection were obtained with the atomic sample for which we imposed the density-dependent exclusion volume r min = ρ -1/3 /π on each atom to suppress the influence of pairs of very close atoms, while being able to reach high densities without introducing a significant amount of positional correlations. Therefore, we have to check that the choice of density-dependent exclusion volume does not cause the observed density effects on subradiant lifetimes. This includes verifying that with our choice of the exclusion volume, remaining pairs do not play a significant role in the decay dynamics and ensuring that the observed density effect cannot be attributed to the small amount of positional correlations introduced by the exclusion volume.

II.2.3.1. Influence of subradiant pairs on collective dynamics

A safe method to remove the close pairs of atoms from the macroscopic sample is to use an exclusion radius r min = π/k 0 as a minimum inter-particle distance. When we impose such exclusion volume on atoms, the decay rates of pair modes are very close to Γ 0 (Fig. II.2), i.e. strong super-and subradiant pairs are not being created in the sample. However, this is only appropriate for investigating dilute samples [Guerin 2016,Araújo 2018b] since it is not possible to reach densities higher than ρλ 3 ∼ 8 with r min = π/k 0 . Here we thus use a less stringent condition with the density-dependent exclusion volume r min = ρ -1/3 /π [START_REF] Moreira | [END_REF]].

However, for higher densities, with this exclusion volume we do not prevent all the super-/subradiant pairs from being created in the sample, only the closest ones. It is then necessary to check that the remaining pairs are not responsible for the observed effects.

Cut-off of the pair branches

First, let us illustrate how the exclusion volume is limiting the extent of pairs. In For the left panels, the drive is red-detuned ∆ 0 = -15Γ 0 , while for the right panels, it is blue-detuned ∆ 0 = +15Γ 0 . In addition to collective subradiant modes (IPR 1), there are two subradiant branches (and two superradiant branches) that consist of modes of close pairs of atoms (IPR ∼ 1/6 -1/2). Branches labeled as and ⊥ represent pair eigenmodes corresponding to two atomic dipoles that are oriented along the direction of their alignment and orthogonally to their alignment axis, respectively. Two subradiant pair branches are asymmetric: they have different frequencies and lifetimes for a given inter-atomic distance. Without any exclusion volume, the extent of branches is very long, and the driving field is resonant (∆ n = 0) with some of the pair modes (it depends on the sign of the detuning which branch crosses the resonance ∆ n = 0). This results in highly populated modes with ∆ n = 0, and the entire late-time decay dynamics is strongly governed by their decay [see Fig. II.8(a)]. Since we are interested in collective subradiance, we use the exclusion volume to restrict the pair branches' extent to prevent strong coupling of the drive with pair modes. The extent of the branches with the exclusion volume r min = ρ -1/3 /π is illustrated with the black solid lines that correspond to the analytical expressions given by Eqs. (II.21) ( ) and (II.22) (⊥) with the cut-off (the maximum extent) corresponding to the exclusion volume r min . As r min is smaller for larger ρ, the extent of the branches with the imposed density-dependent exclusion volume is getting larger with increasing density. Red solid lines (left panels) and blue solid lines (right panels) highlight the resonance ∆ n = 0. The color code corresponds to the inverse-participation ratio (IPR) given by Eq. (II.18).

super-/subradiant pairs have been created, and therefore, the extent of the pair branches (with IPR∼ 1/6 -1/2) is extremely long. Note that here we show pair branches only in a finite range, since for very short distances of two atoms, the frequency shifts of pair modes are very large. For the red-detuned driving field, the ⊥ subradiant branch crosses the resonance ∆ n = 0, while for the blue-detuned driving field, that is the case for the subradiant

Red-blue asymmetry of the pair branches

Another test which we performed to check the role of pairs is the red-blue asymmetry. As already mentioned, for a given inter-particle distance, the pairs result in two subradiant the subradiant decay is slower for the positive detuning than for the negative detuning.

On the other hand, with density-dependent exclusion volume r min = ρ -1/3 /π, decay curves corresponding to positive and negative detuning are identical; there is no red-blue asymmetry.

Additionally, we show the results obtained with an exclusion volume defined by a fixed minimal distance r min = 0.5/k 0 . For ∆ = -15Γ 0 , the corresponding decay curve is identical to those obtained with r min = ρ -1/3 /π. However, for ∆ = +15Γ 0 , the subradiant decay at very detunings. This absence of red-blue asymmetry for the obtained data with r min = ρ -1/3 /π confirms that there is no influence of pairs on the collective long-lived dynamics.

Polarization of the scattered light

Furthermore, another qualitative difference between collective long-lived modes and subradiant pairs is the polarization of the scattered light. To investigate the polarization properties of the scattered light, the system has been driven to the steady-state with a circularly polarized driving field, and we compute the scattered electric field at θ sc = 90 • . Furthermore, we decompose the scattered field in the far-field limit (where the radial component of the electric field ceases to exist) into two orthogonal components, i.e., polarization channels: parallel and orthogonal to the scattering plane. Then we compute the intensity of the scattered light in these two polarization channels. Note that the scattering plane is defined by the wavevector of the driving field k 0 and the observation direction n.

In the case of a single atom, i.e., single scattering event at θ sc = 90 • , the scattered light has a linear polarization orthogonal to the scattering plane. However, one can see in Therefore, the scattered light coming from subradiant pairs is (partially) polarized, while collective subradiant modes yield fully depolarized scattered light. We have checked that the long-lived decay dynamics is depolarized for all system parameters studied in this chapter with r min = ρ -1/3 /π, further confirming that the corresponding subradiant decay dynamics is indeed due to collective modes and that pairs do not impact the results shown in Figs. II.4

and II.5.

Interestingly, we note that collective long-lived modes yield depolarized light while the superradiant early decay measured at θ sc = 90 • is linearly polarized (orthogonal to the scattering plane). This is fully consistent with an interpretation of superradiance based on a single-scattering event embedded in an effective medium, as discussed in [Weiss 2021]. On the contrary, the observation of the full depolarization of the long-lived decay suggests the multiple-scattering interpretation of subradiance, as discussed in [Fofanov 2021].

On-resonance collective decay dynamics

We also computed the decay dynamics with a drive on resonance, ∆ 0 = 0. In this case, even without any exclusion volume, we do not expect to see the influence of pairs since the driving field is strongly coupled to long-lived resonant collective modes with ∆ n = 0, rather than to The drive is on resonance, ∆ 0 = 0, and the light is collected at θ sc = 45 • . The long-lived dynamics is slower for ρλ 3 = 8 than for ρλ 3 = 16. This confirms that the observed reduction of subradiant lifetimes at increasing densities in vectorial coupled-dipole model with the off-resonant drive is not caused by subradiant pairs or by positional correlations. Here the small discrepancy between the decay curves obtained with and without an exclusion volume can be attributed to the positional correlations that arise when the exclusion volume is imposed on the atoms.

Interplay between collective and subradiant pair modes for ∆ 0 = 0 Together, all the tests above clearly demonstrate that the reduction of the subradiant lifetime with increasing density cannot be attributed to pair physics.

For the sake of providing guidelines on how to remove subradiant pairs from the collective decay dynamics of dense atomic samples (for sample densities even beyond those studied in this chapter), let us summarize our observations on the interplay between collective and pair modes.

As discussed, the temporal decay dynamics results from the competition between the modes of [START_REF] Cipris | Subradiance with Saturated Atoms: Population Enhancement of the Long-Lived States A[END_REF]. With this density-dependent exclusion volume, the minimum imposed inter-particle distance is smaller for larger density; hence, the extent of the pair branches is getting larger with increasing density. Therefore, if one would use the density-dependent exclusion volume r min = ρ -1/3 /π for densities significantly larger than those in this chapter, the drive detuning ∆ 0 should be increased so that the branches with the cut-off do not cross the resonance ∆ n = 0.

By reducing the extent of the pair branches with the exclusion volume before they reach the resonance ∆ n = 0, we make sure that the pairs are not significantly strongly driven compared to the collective modes. Still, there is a competition between collective modes and remaining pair modes, which might result in observing decay dynamics that is partially governed by

II.2.3.2. Influence of positional correlations

Imposing an exclusion volume actually induces some correlations in the atomic positions.

This correlated disorder, in turn, affects the light scattering properties (see, e.g., [Lax 1951, Fraden 1990, Rojas-Ochoa 2004a, Wang 2020]). In particular, increasing the density at fixed exclusion volume or increasing the exclusion volume at fixed density enhances the positional correlations. In Therefore, we can conclude that in our work correlated disorder only has a small marginal role on lifetimes and is not causing the density effects reported in Sec II.2.2.

II.2.4. Comparison with scalar light

The density-induced reduction of the subradiant lifetime observed in

r min = ρ -1/3 /π, θ sc = 45 • , ∆ (s) 0 = -15Γ (s) 0 .
examples of decay curves from the two models for several b 0 and ρλ 3 . At smaller ρλ 3 , the decay curves obtained from the vectorial and scalar CDEs are almost identical for a given onresonant optical thickness. On the other hand, for larger ρλ 3 , we observe a slower subradiant decay in the scalar model than in the vectorial model. Moreover, in panel (b) we show the direct comparison of subradiant lifetimes of the vectorial and scalar models obtained at late times as a function of b 0 for several densities of the sample. As can be observed, there is a qualitative difference between the scalar and vectorial subradiant lifetimes at increasing densities. While in the vectorial case we observe a decrease of τ sub with ρλ 3 , the behavior is opposite with the scalar model: τ sub increases with ρλ 3 for a given b 0 .

This qualitative different behaviour between the scalar and vectorial models clearly demon-strates that the reduction of subradiant lifetimes with the density is due to the near-field part of the dipole-dipole interaction, which is absent from the scalar model.

Note that in the scalar model, both the resonant optical thickness and the natural decay rate differ by a factor 2/3 from their vectorial version [b

(s) 0 = 2N/(k 0 R) 2 = (2/3)b (v) 0 and Γ (s) 0 = d 2 k 3 0 /2 π 0 ]
, so simulations for a given optical thickness and density involve different atom numbers in the scalar and vectorial models. We also restrict ourselves to densities up to ρλ 3 ∼ 20, to avoid the Anderson-localized regime for scalar light [Skipetrov 2016a].

Although until this subsection we used Γ 0 and b 0 notations for the natural decay rate and optical thickness in the vectorial model, here we use Γ and we have checked that the smaller constant exclusion volume r min = 0.5/k gives the same results, as well as the positive detuning ∆ 0 = +15Γ 0 . That, together with the fact that the obtained subradiant lifetimes are longer than the subradiant lifetimes of the remaining pairs with the exclusion volume, confirms that the increase of τ sub with increasing density in the scalar model cannot be attributed to the subradiant pairs. Moreover, a test with the driving field on resonance and without any exclusion volume also yields increased subradiant lifetimes with increasing density for a given b 0 in the scalar model, which rules out the influence of positional correlations.

As was discussed in Sec. II.2.2 for the vectorial model, to increase the density of the sample for a given optical thickness, the size of the sample has to be decreased. Therefore, we should also check for the scalar model, if the observed increase of subradiant lifetimes with increasing density is a density effect or if it stems from the decrease of the sample size. In 2ρλ 3 [panel (b)]. We observe that for a given optical thickness, the scalar subradiant lifetimes are independent of ρλ 3 and k 0 R for smaller densities/larger sample sizes, and then there is an increase of subradiant lifetimes with increasing ρλ 3 /decreasing k 0 R. Contrary to the vectorial model (Fig. II.5) where it was clear that the reduction of subradiant lifetimes is a density effect, here for the scalar model, it is not conclusive: there is no clear k 0 R or ρλ 3 threshold for which the subradiant lifetimes start to increase, as there was in the vectorial case. , in the scalar model that is the case for lower intensity-levels I(t)/I(0) < 10 -6 , i.e., for even later times.

Intermediate times

Let us emphasize here the importance of focusing at sufficiently late times of temporal decay dynamics to be able to identify the b 0 enhancement of subradiant lifetimes, and density effects.

Note that in Figs. II.11-II.12, where we made comparison with the scalar model, subradiant lifetimes have been obtained using the lower intensity-level fit range, i.e., later in time than in the previous subsections concerning only the vectorial model. Moreover, we have already mentioned that in the vectorial coupled-dipole model, we obtain qualitatively the same results

for any appropriate fit range below relative intensity level I(t)/I(0) = 10 -4 . However, that is not the case for the scalar model.

To explain this, let us first recall the importance of fitting the computed decay curves at sufficiently late times. In Note that for these intermediate fit windows, the density effects on subradiant lifetimes are negligible, suggesting that the density effects impact the longest subradiant modes for a given b 0 .

Therefore, one has to make sure that fitting of the decay curves is done at sufficiently late times to obtain the subradiant lifetimes with the characteristic b 0 cooperative enhancement and also to observe (possible) density effects on subradiant lifetimes. In the previous subsection, by comparing scalar and vectorial decay dynamics, we came to the conclusion that the observed decrease of subradiant lifetimes with increasing density for a given on-resonant optical thickness in the vectorial coupled-dipole model is due to the near-field terms (∝ 1/r 2 and ∝ 1/r 3 ) of the dipole-dipole interaction.

II.2.5. Interpretation of the decay-dynamics results

To understand why the near-field terms are detrimental for subradiance, as well as to get further insight into the opposite behaviors for scalar and vectorial light, we turn to analysing the density evolution of the spectrum of eigenvalues (λ n = -Γ n /2 + i∆ n ). Let us recall that for the vectorial coupled-dipole model, we obtain eigenvalues by diagonalizing a 3N × 3N Because of these differences in the spectrum of vectorial and scalar model, we can conclude that in the presence of near-field terms, the increase in density [from panel (a) to (b)] is characterized by a strong broadening of the spectrum, and as a consequence the subradiant tail of the distribution is reduced; the broadening due to the near-field dipole-dipole interaction is at the origin of the reduction of subradiance. While we here focus on the dynamical features of the scattering, it is interesting to note that inhomogeneous broadening, due to the near-field terms of the dipole-dipole interaction, has also been identified as a limiting mechanism for the increase of the refractive index at large densities [Andreoli 2021].

Let us now comment on the analogy of our results with the superradiance of a fully inverted system (beyond the single-excitation limit). Originally studied without accounting for the near-field terms, the superradiant cascade was addressed in two different regimes: for subwavelength and macroscopic clouds [Dicke 1954]. In the case of a subwavelength cloud, it can be assumed that a unique light mode is coupled to the sample, and the cooperativity parameter describing this coupling is the number of particles N . Differently, in the case of a macroscopic cloud, such as studied here, the sample geometry plays a role and the resonant optical thickness b 0 is the relevant cooperativity parameter. The importance of near-field terms in subwavelength dense samples was later recognized, showing that they are detrimental to superradiance [Friedberg 1972, Friedberg 1973, Friedberg 1974, Gross 1982].

Coined "Van der Waals dephasing" due to the 1/r 3 decay they exhibit, the near-field term breaks the symmetry that was central to Dicke's approach, since he would assume the system to decay through a series of symmetric states.

Note that there are two distinct symmetry-breaking effects. The first one occurs when the sample size is increased to become comparable or larger than the optical wavelength: even in the scalar light approximation, the atomic dipoles couple to several optical modes, and the cooperativity parameter is then given by the resonant optical thickness b 0 rather than the particle number N , for macroscopic clouds. This effect was already discussed in the seminal paper by Dicke [Dicke 1954]. The second symmetry-breaking effect is Van der Waals dephasing that corresponds to the strong energy shifts induced by the near-field terms [Gross 1982]. The inhomogeneous broadening resulting from these terms leads to a reduction of the cooperativity. For subwavelength samples, the reduction of the N -factor enhancement characteristic of superradiance stems from the increase of the density rather than a modification of the system size. Although we here address the single-excitation regime and macroscopic clouds, the behaviour of the subradiant lifetimes that we have observed in this chapter is very consistent with the picture developed for the superradiant cascade. For that reason in [Cipris 2021b] we give to our observed density effects on subradiant lifetimes the term Van der Waals dephasing: the effect on cooperativity due to the density increase, i.e., the rise of the near-field terms, beyond the size effects.

As was discussed in the previous subsection, for the scalar model it is not clear whether the observed effect of increase of subradiant lifetimes with increasing density for a given b 0 , is a density or a size effect (Fig. II.12). Therefore, we give two interpretations of this subradiant lifetime increase for scalar light.

When the density is increased for a fixed on-resonant optical thickness, the system size reduces (since R ∝ b 0 /ρ). Since for subwavelength clouds cooperativity scales with N , it could be that upon decreasing the sample size, the number of particles N starts to compete with the on-resonant optical thickness b 0 as the cooperativity parameter. This competition then, could be responsible for the observed increase of the subradiant lifetimes in the scalar model. On the other hand, as we are increasing the density of the sample, we are approaching the localization regime, for which the phase transition is formally reached only for ρ c λ 3 ≈ 20 [Skipetrov 2014, Bellando 2014]. Therefore, weak localization processes could also be responsible for this increase in subradiant lifetimes.

II.3. Conclusion

In conclusion, in this chapter we have reported on the density effect on subradiance in macroscopic atomic clouds in the single-excitation limit. For densities ρλ 3 5, scaling of subradiance with b 0 fails, and subradiant lifetimes are reducing with increasing density for a given on-resonant optical thickness. We have concluded that this reduction of subradiant lifetimes is due to near-field dipole-dipole interaction terms that induce an inhomogeneous broadening which acts against cooperative effects. This Van der Waals dephasing for subradiance presents very similar features as the one discussed for superradiance [Gross 1982] in subwavelength clouds and considering fully-inverted system.

Note that for the densities studied throughout this chapter, we have checked that superradiance is not substantially affected. This difference can be attributed to the very different time scales involved in each phenomenon. Indeed, superradiance is very fast, with time scales shorter than Γ -1 0 , so a very strong broadening is required to affect the dynamics over these short times. Differently, subradiance corresponds to modes with lifetimes of many units of Γ -1 0 , making them much more sensitive to the broadening induced by the near-field terms. This analysis is confirmed by the fact that superradiance is more robust than subradiance against inhomogeneous broadening induced by thermal motion [Weiss 2019, Weiss 2021].

Finally, we have observed that for densities ρλ 3 < 5, the temporal decay dynamics is the same for both models, scalar and vectorial. This implies that for such low densities, the dipole coupling is indeed provided by the long-range (∝ 1/r) dipole-dipole interaction, with negligible density or polarization effects. Therefore, the problem of collective decay dynamics of dilute atomic samples can be addressed with the scalar coupled-dipole model. Moreover, in our experiment the typical density of cold atomic rubidium cloud is ρλ 3 ∼ 10 -2 . Hence, our study here suggests that for such dilute cloud, the density, i.e. near-field, effects are indeed negligible.

Note that the coupled-dipole model introduced in this chapter does not include all the effects, such as light assisted collisions or the temperature, i.e. thermal motion, of cold atoms. Although light assisted collisions might introduce limitations on obtaining high densities of cold-atom sample, high densities ρ ∼ 10 14 cm -3 (ρλ 3 ∼ 10) have been obtained

for Strontium (Sr) [Ferrari 2006] and Ytterbium (Yb) [Takasu 2003] cold atoms which also have the J g = 0 → J e = 1 transition, as the one considered in our coupled-dipole model. Moreover, concerning the thermal atomic motion, our recent experiment demonstrated that subradiance is robust against thermal decoherence for a large temperature range of cold-atom sample [Weiss 2019].

Chapter III

Signature of Anderson localization of light in 3D

Wave propagation in a random medium is a phenomenon that occurs throughout our daily lives; the most common examples are multiple scattering of light in clouds, milk, paper, etc. In fact, the wave transport in a random medium depends on the relation between several characteristic lengths: the wavelength λ, the size of the sample L, and the typical distance between two consecutive scattering events -the scattering mean free path . In the single-scattering regime (ballistic regime), where the sample size is smaller than the mean free path, L , the wave undergoes at most one scattering event before leaving the medium.

The analogous condition of this regime is the optically thin medium with the optical thickness

b(∆ 0 ) = b 0 /(1 + 4∆ 2 0 /Γ 2 0 )
1. On the contrary, when the mean free path is much smaller than the sample size, L, the wave is scattered numerous times before exiting the random medium. This is the so-called multiple-scattering regime. Within this regime, there are two different types of wave transport: diffusive transport ( λ) and localized regime ( λ).

When λ, the interference of multiply scattered waves is negligible; each scattering event can be described by a random walk with a step (incoherent wave propagation) and the wave transport in a random medium is described by the diffusion equation. The random medium in this regime is often called weakly disordered medium. However, upon increasing disorder (randomness) of the medium such that the mean free path becomes comparable to the wavelength, coherent effects in the wave propagation emerge and the diffusive model is not suitable anymore for the description of wave transport. In a sufficiently strongly disordered medium ( λ), there is no diffusive wave transport, i.e. waves become localized within a finite region.

This wave interference phenomenon in a strongly disordered medium that manifests in a breakdown of diffusive wave transport is known as Anderson localization [Anderson 1958]. It is assumed to be a universal wave phenomenon: it should occur for all types of classical and quantum waves. Throughout the years Anderson localization has been observed for many types of waves, both in three-dimensional (3D) and low dimensional (1D and 2D) systems.

However, in spite of decades of immense search, Anderson localization of light in 3D strongly scattering media is still experimentally unobserved. As it is going to be discussed in the following sections, one of the challenges has been identifying an unambiguous signature of light localization. In this chapter, we are going to present the research based on a numerical coupled-dipole model that has led us to identify an unequivocal signature of light localization in 3D in the framework of cold atoms as the disordered scattering medium. Note that for this system that we have been considering: cold-atom sample, the disorder is in the random position of atoms (off-diagonal disorder). In that case the disorder strength, or better to say the scattering strength, is characterized by the spatial density of the sample which is inversely proportional to the scattering mean-free path.

In

III.1. Brief review on Anderson localization

P. W. Anderson has discovered that when treating electrons as quantum particles, i.e., considering their wave character, the diffusive electron transport in a random potential can be significantly altered if the randomness, i.e., the disorder is sufficiently strong. In his seminal work [Anderson 1958], Anderson considered a tight-binding model of an electron with nearest-neighbor hopping in a 3D lattice that contains disorder, i.e., a lattice with random on-site energies. Anderson pointed out that if the degree of randomness (disorder) in the lattice is sufficiently large, localization occurs: the electron is trapped in a deep well of the random lattice potential as its wave function becomes exponentially localized. In other words, Anderson showed that, when the disorder potential exceeds a critical value, the classical diffusive electron transport in a 3D lattice comes to a complete halt: the medium essentially ceases to conduct the charge and becomes an insulator. This disorder-driven phenomenon was coined as Anderson localization. It is also often referred to as strong localization.

Therefore, Anderson predicts two regimes, i.e. phases of wave transport in a medium containing disorder: the diffusive, conducting (metallic) regime and the localized, insulating regime. In the diffusive regime, the states responsible for diffusion are extended and the wave function extends throughout a disordered medium. On the contrary, in the localization regime, where the diffusive wave transport is absent, there are only localized states. In particular, those states are exponentially localized over a typical distance ξ known as the localization length; in the localized regime, the wave function decays exponentially from the origin point r 0 :

ψ(r) ∝ exp(-|r -r 0 |/ξ). (III.1)
As is going to be further discussed in the following subsection, in low dimensional (1D and 2D) disordered systems the waves are always localized no matter how small the disorder is, assuming the waves are noninteracting and the system size is larger than the localization length. However, in three dimensions, a phase transition from diffusive regime to the localized regime can occur upon increasing the disorder of the scattering medium. This transition is often called Anderson transition, while the critical energy that separates the localized regime from the extended regime is called mobility edge.

III.1.1. Anderson localization: disorder-induced wave interference phenomenon

There are a couple of key points to bear in mind when it comes to Anderson localization. For most of these diffusive paths, the interference between probability path amplitudes vanishes on average. b) However, each multiple-scattering path has its reversed counterpart such that wave is scattered by the same set of scatterers but in the opposite order: time-reversal paths. c) In the case of closed loops, when the wave returns to the inital scattering center (A=B), the accumulated phase difference between time-reversed loops is zero and interference effects survive the ensemble average. The probability of return to the inital scattering point A is enhanced due to constructive interference on time-reversed loops. Contribution of such closed loops in wave transport in weak disorder decreases the diffusion coefficient (weak localization). By increasing the disorder, probability for closed loops increases, thus further decreasing the diffusive wave transport. At a critical amount of disorder, when the diffusion coefficient becomes zero, the phase transition from extended-wave regime to localized regime arises. d) Another way for interference of time-reversed loops to survive the configuration average is if k f = -k i , where k i is the wavevector incident on the disordered medium and k f is the outgoing wavevector. These interferences on time-reversed loops give rise to the coherent backscattering, the enhancement of the multiply scattered intensity in the backward direction by a factor of two.

An intuitive way to understand how the localization of waves occurs is the following. We are going to consider waves that are multiply scattered by identical point-like scatterers with random positions (disordered system). The wave propagation from point A to point B can occur along infinitely many different paths, as illustrated on Fig. III.1(a). For each path p there is a probability amplitude a p for the wave transport from A to B. To obtain the total probability P that the wave gets to point B from A one has to sum the probability amplitudes of all possible paths a p and then square the sum:

P = | p a p | 2 = p |a p | 2 + p =p (a p a * p + a p a * p ).
Here, the first term represents the incoherent contribution, i.e. it describes non-interfering paths, while the second term corresponds to the contribution due to interference of path

= (k i + k f ) • (r N -r 1 ), (III.2)
where k i and k f are incident and outgoing wavevectors, respectively, while r

return to A is P = |a 1 | 2 + |a 2 | 2 + 2a 1 a * 2 = 4|a| 2
, twice as big as if only incoherent contribution would be considered. Therefore, the probability of return to the origin point A due to the constructive interference is larger than the probability of propagating away from A (incoherent wave propagation), i.e., the probability of the diffusive wave transport through the sample is decreased. This existence of interference of multiply scattered waves has been initially demonstrated in weakly disordered metallic films [Bergmann 1984].

These interference effects of time-reversed loops in the limit of weak disorder reduce the diffusion coefficient and hence the conductivity. The enhanced return probability due to the formation of time-reversed closed loops can lead to a decrease in diffusive wave transport, since waves stay in the medium for a longer time than what is expected considering only the incoherent diffusion process. This effect that manifests itself in a decrease of the diffusion coefficient in the weak disorder limit is called weak localization and it is considered to be a precursor of strong (Anderson) localization effect. As disorder increases, the probability of formation of time-reversed loops increases; the density of these constructively interfering time-reversed loops increases with decreasing mean free path . Above the critical amount of disorder, when the number of time-reversed loops diverges, i.e. when the return probability reaches unity, the transition to Anderson (strong) localization occurs. In other words, when the formation of time-reversed loops suppresses the occurrence of diffusive paths in strongly disordered medium, the diffusion coefficient becomes zero and the wave becomes trapped, i.e., localized inside the medium [Aegerter 2009].

As already mentioned, the phase transition from conductor (extended waves; diffusive wave transport) to insulator (localized regime) is a characteristic of a 3D disordered medium. In 1D and 2D disordered media, waves can be localized, even for a very weak disorder. This can be understood if we recall that the diffusive wave transport in 1D and 2D is a recurrent process: the diffusive particle returns infinitely many times to the point of origin, forming the time-reversed closed loops with extremely high probability (recurrent random walk: the probability of returning to the origin point is equal to 1 in 1D and 2D [Woess 2000]). In principle, the diffusive wave propagation in 3D disordered medium is not recurrent and the probability of forming constructively interfering closed loops in a weakly disordered 3D medium is much lower than in 1D or 2D. Hence, in 3D the localization does not occur unless the critical amount of disorder is reached such that the probability of returning to the starting point is significantly increased. The diffusion constant can consequently be reduced to zero and the wave is "trapped" on closed loops, i.e. it can not escape from its origin point.

Therefore, in 3D, the coherent wave propagation essentially becomes recurrent in the localized regime [START_REF] Akkermans | [END_REF]].

Another way for waves on time-reversed paths to be perfectly coherent is if the incoming and outgoing wave have the same direction, i.e. k f = -k i [Eq. (III.2)], as shown on Fig. III.1(d).

Then, the constructive interference between waves counterpropagating along a given multiple scattering path leads to the twofold enhancement of the intensity reflected from the disordered sample in the backscattering direction. This phenomena is known as coherent backscattering (CBS). Since coherent backscattering effect can occur for a weak disorder and it lowers the diffusion coefficient, sometimes it is considered as a weak localization phenomena.

Let us note that the explanation provided here on the occurrence of Anderson localization based on the time-reversed closed loops is the simplest and the most intuitive: the problem of Anderson localization in the context of loops is much more involved. Finally, remark has to be made that this intuitive picture of non-vanishing constructive interference for waves propagating on time-reversed paths only allows for an estimation of the lowest order correction (weak localization) to the diffusion constant [Vollhardt 1987]. To study higher order corrections or Anderson transition, other theoretical approaches, which will be mentioned in the following subsection, have to be considered.

III.1.2. Localization criteria

III.1.2.1. Scaling theory of localization

Thouless criterion

In 1970s, D. Thouless and his co-workers introduced several physical ideas important for the formulation of a scaling description of the Anderson localization [Edwards 1972,Thouless 1974[START_REF] Licciardello | [END_REF]. They pointed out that the sensitivity of states to boundary conditions of disordered system is one way1 of distinguishing between the delocalized, extended states, whose wavefunction extends over the whole system, and localized states whose amplitude falls off exponentially away from the localization centre. As boundary conditions reflect the nature of the coupling to the environment, they can play an important role for a finite size system. Any change in boundary conditions significantly alters the energy spectrum of extended states. On the other hand, localized states will not be (notably) sensitive to the change in boundary conditions: their sensitivity to the boundary will be exponentially small since only the exponential tail of the (exponentially) localized state will "feel" the boundary, asuming that the localization centre is not near the boundary.

Let us consider a model of a disordered system as a d-dimensional hypercube of size (2L) d , made by putting small cubes of size L d (with L

). The eigenstates for the particle in the

g(L) = G e 2 = ∆E δW . (III.3)
Thus, via Einstein relation the dimensionless conductance is given by the energy coupling/splitting ratio, the quantity that is very often referred to as Thouless number

g Th = ∆E δW . (III.4)
Therefore, the dimensionless conductance can be used to distinguish between extended and localized states in 3D: g(L) 1 in the diffusive (extended wave) regime, while g(L) 1 for the localized states. Moreover, it is expected that the transition from extended to localized regime occurs when ∆E and δW are of the same order of magnitude, i.e. when g c (L) ∼ 1.

This condition for Anderson transition to take place is known as Thouless criterion.

Although Thouless considered electronic systems, his ideas can be extended to classical waves.

The Thouless number for classical waves can be defined as the ratio of the typical frequency width ∆ω to the average spacing of eigenmodes δω:

g Th = ∆ω δω . (III.5)

One-parameter scaling theory

The work by Thouless and collaborators led to conclusion that the conductance, as a proxy for the disorder in the system, is the single parameter that controls the behaviour of the disordered system, i.e. the localized or extended nature of eigenstates, as the system varies in size. Moreover, Wegner found that for a single-parameter scaling theory to hold, the phase transition from diffusive to localized regime has to be continuous [Wegner 1976]. Based on this ideas, in 1979 the "Gang of Four" (E. Abrahams, P. W. Anderson, D. C. Licciardello and T. V. Ramakrishnan) formulated the scaling theory of localization [Abrahams 1979], a phenomenological theory that interpolates between extended and localized states.

They have considered combining b d cubes of size L d (with the same disorder properties) into a large cube of size (bL) d . Under the assumption that the change in the conductance g(L) with the system size is determined only by the conductance itself, they stated that the conductance of the large system g(bL) has to be the function of the conductance g(L) of the L d cube, which can be written as g(bL) = f (b, g(L)). We can rewrite this as

g(L + bL) -g(L) bL = f (b, g(L)) -g(L) bL , (III.6)
such that b quantifies the system size increment. By taking the limit b → 0 we obtain:

β(g) = d lng(L) d lnL , (III.7) with β(g) = lim b→0 f (b, g(L)) -g(L)
gL . The scaling function β(g) [Eq. (III.7)], which describes how the conductance changes when the system varies in size, only depends on the conductance itself, and not on the size L nor any other microscopic quantity. In other words, the change of conductance when the system size increases is determined by the value of the conductance at the previous length scale.

There are two different asymptotic behaviours of conductance g(L) depending on whether the state is extended (g 1) or localized (g 1). In the case of an extended regime, the conductance is given by Ohm's law:

g(L) = e 2 σL d-2 . (III.8)
On the other hand, in the regime g 1 where the states are exponentially localized with the localization length ξ L, the conductance is described by g(L) ∝ exp(-L/ξ). (III.9)

Hence, the asymptotic form of the scaling function β(g) in these two regimes is

β(g) ∝      d -2 -a/g for g 1, 2
ln(g) for g 1.

(III.10) From Eqs. (III.8) and (III.9), one can see that in a 3D system the conductance increases with the system size when the states are extended, but on the other hand for a localised state, g(L) decreases with L. Therefore, if the scaling function is positive, states are extended, while β(g) < 0 indicates that the states are localized. However, in lower dimensions g(L)

always decreases with L, meaning that in those systems there is no phase transition from diffusive to localized regime, the states always localize.

From the asymptotes of β(g) [Eq. (III.10)] and based on the assumption that β(g) is a monotonous and continuous function, one can sketch the universal curve β(g) for dimensions for different dimensions as a function of the dimensionless conductance g. The scheme was constructed from known asymptotes of β(g) [Eq. (III.10)] in localized and extended regime and by assuming continuity and monotonicity of β(g). Arrows indicate the direction of variation of g as the system size L increases. Such scaling functions imply that for d = 3 there is a phase transition from extended to localize regime, which occurs at critical point g c . In lower dimensions there is no phase transition: states always localize in the thermodynamic limit L → ∞. to the extended regime is reached. Let us now consider a strongly disordered system such that the starting point is g(L) < g c , and thus β(g) < 0. Then, while increasing L, g(L)

decreases and one moves downward on the β(g) curve. At large enough length scales, β(g) corresponds to the scaling function for localized states. Therefore, for a fixed disorder, if g(L) > g c , the system will flow into the limit of extended regime as L increases. On the other hand, if g(L) < g c , the system will flow into the limit of localized regime as L increases.

In other words, depending on the amount of disorder the system is either in the localized or in extended regime in thermodynamic limit L → ∞. The scale invariant critical point g c (β(g c ) = 0) separates those two stable phases: it corresponds to the point of Anderson transition from extended to localized regime. As the direction of variation of β(g) points away (upward and downward) from g c , this critical point is an unstable fixed point.

While the scaling theory predicts the existence of phase transition for 3D systems, that is not the case for lower dimensions. In 1D and 2D, β(g) is always negative and therefore the conductance always goes to its limit of the localized regime [Eq. (III.9)] when the size of the system increases. In other words, in d 2 the system always reaches the localized regime, no matter what the degree of disorder is, as long as the length scale is large enough, i.e. L ξ.

Moreover, d = 2 is the marginal dimension and depends sensitively on the conditions. In the above description, the system is time-reversal invariant with no interactions and the scaling theory predicts that there is no phase transition for d = 2.

Let us now discuss how the effect of the disorder strength is established in the scaling theory of localization [Sheng 2006]. The scaling hypothesis states that the change in conductance with system size represented by scaling function β(g) depends explicitly only on the conductance itself. However, besides having the L dependence, g also depends on the amount of disorder in the system. Then, for the scaling hypothesis to hold, the effect of increasing or decreasing the disorder has to be completely compensated by varying L. A way for this to happen is if

g is a function
ξ(η) = (η -η c ) -ν , η > η c , (III.12) σ(η) = (η c -η) ν , η < η c . (III.13)
Note that since D ∝ σ, the diffusion coefficient has the same power-law behaviour as the conductivity. Therefore, while the localization length on the localized side diverges for η → η c , the diffusion coefficient in the delocalized phase approaches zero for η → η c and vanishes on the other side of the mobility edge in the localized phase, i.e. for η → η c .

The scaling theory of localization does not depend on microscopic properties of the system such as the type of the wave, type of the scatterers (disorder), etc. This is due to the dimensionless conductance which does not depend explicitly on any microscopic quantity of the system, but only implicitly through, for example, the correlation length ζ(η) [Eq. (III.11)].

Therefore, the scaling functions β(g) for d = 1, 2, 3 have a universal character and the scaling theory should hold for both classical and quantum waves no matter of the microscopic properties of disordered sample. However, although the scaling theory predicts the Anderson transition in 3D, it does not guarantee the existence of the β < 0 branch of the scaling function, i.e. the existence of localized states. As was discussed above, the phase transition between diffusive and localized states will arise at the critical disorder strength, which marks the dominance of interference effects in wave transport. It might happen that due to some microscopic properties of the scatterers it is hard to achieve strong enough scattering, or some underlying physical mechanisms could exist which are causing the dephasing of waves in the disordered system and preventing the wave interference effects to dominate the transport, thus prohibiting localization. This is the case for classical waves, particularly light waves for which the observation of strong localization is still lacking. For example, in the case of dielectric scatterers it is hard to produce a strong enough disordered sample to reach the critical point. On the other hand, it has been shown that in the case of point scatterers with degenerate internal structure and light waves with vectorial nature there is no phase transition in 3D [Skipetrov 2014]. The light localization will be discussed in details in Sec.

III.2. The point here is that if the localization is possible, the system should have the same universal scaling function β(g), as in Fig III .3.

The scaling theory has provided an important framework in understanding the localization phe-nomena. For a review on one parameter scaling theory, see Refs. [Lee 1985,Sheng 2006]. Apart from the scaling theory, there are other approaches that allow a description of localization phenomena: the self-consistent theory [Vollhardt 1980a, Vollhardt 1980b, Vollhardt 1982[START_REF] Vollhardt | [END_REF], the random matrix theory [Wigner 1951, Dyson 1962], the supersymmetry nonlinear σ model [Efetov 1983, Efetov 1996], etc. For a review on different theories, see

Refs. [Lee 1985, Evers 2008, Abrahams 2010].

III.1.2.2. Ioffe-Regel criterion

A well-known qualitative criterion as to when the transition from diffusive to localization regime occurs in 3D disordered systems is that of [START_REF] Ioffe | Progress in Semiconductors vol. 4 Chapter: Non-crystalline, amorphous and liquid electronic semiconductors[END_REF]]. It suggests that, when the typical distance between two consecutive scattering events, i.e. the mean-free path , is of the order of the wavelength λ = 2π/k, interference effects between multiple scattered waves become substantial inside the disordered medium and the Anderson localization then takes place when

k (k ) c ∼ 1. (III.14)
Here, (k ) c ∼ 1 represents the critical disorder strength for which the Anderson transition occurs. The localization condition given by Eq. (III.14), that estimates how strongly the disorder medium has to scatter to localize the wave, is known as Ioffe-Regel criterion. This criterion discriminates between diffusive regime (k 1) and localized regime (k < 1). An intuitive way of understanding the Ioffe-Regel criterion is the following. The product of the wave number and the mean free path k quantifies the number of oscillations of a wave between two scattering events. Hence, the condition k < 1 implies that the wave scattered from one scattering center undergoes the next scattering event before it makes one oscillation.

Therefore, due to such strong scattering, the diffusive wave propagation is unsustainable, i.e. the wave cannot be considered to be freely propagating if the scattering mean free path is only a fraction of the wavelength: the classical representation of wave transport has to breaks down.

Ioffe-Regel criterion can be met either by reducing the mean-free path while increasing the scattering strength, i.e. disorder or by reducing the wavenumber k (and with that the energy).

As we will see below, the latter is not achievable for classical waves, but for quantum waves, such as electron waves, it is possible.

Note that the condition for Anderson transition (k ) c ∼ 1 is valid only for 3D disordered systems. As discussed above, in 1D and 2D waves can be localized for any amount of disorder and therefore, for any value of k .

III.2. Strong localization of light

As it was recognized that the interference effects in multiple scattering play an important role in transport properties of electron (see Ref. [Lee 1985] for review on electron localization), i.e.

that the Anderson localization is essentially a wave interference phenomenon, suggestions were made that this disorder-driven localization is a universal wave phenomena that should occur for all types of waves: quantum (electron, phonon, ultracold atoms) and classical (electromagnetic waves such as (visible) light and microwaves, elastic waves such as sound, etc.) [Azbel 1983, Guazzelli 1983, John 1983, John 1984, Kirkpatrick 1985, Anderson 1985, He 1986, John 1987, Flesia 1987] (see Ref. [Lagendijk 1996, Sheng 1990] for review on localization of classical waves). In fact, up to now, Anderson localization has been observed for many types of wave systems: electrons in 1D and 2D [Lagendijk 2009, Dynes 2010] as well as in 3D [Ying 2016],

cold atoms and Bose-Einstein condensates in 1D [Billy 2008,Roati 2008], 2D [White 2020] and 3D [Chabé 2008, Kondov 2011, Lopez 2012, Jendrzejewski 2012, Semeghini 2015], microwaves in 1D [Chabanov 2000] and 2D [Dalichaouch 1991, Laurent 2007], terahertz electromagnetic waves in 1D [Pandey 2017], sound in 2D [Weaver 1990[START_REF] Bretagne | [END_REF]] and 3D [Hu 2008, Aubry 2014,Hildebrand 2014], polaritons in 2D [Sturges 2019,Zhu 2020]. Moreover, there have been many experimental observations of light localization in 1D [START_REF] Topolancik | Experimental Observation of Strong Photon Localization in Disordered Photonic Crystal Waveguides[END_REF],Lahini 2008, Sapienza 2010] and 2D [Schwartz 2007, Riboli 2011]. However, a conclusive experimental evidence on Anderson localization of light in 3D is still lacking, which will be discussed in details below.

III.2.1. Light localization vs electron localization

Although Anderson localization for electrons was suggested in 1958, it was not until the 1980's that the study of Anderson localization was extended to classical waves. Light seemed to be an ideal candidate to study Anderson localization, as photons do not mutually interact, contrary to electrons. Although P. W. Anderson considered non-interacting electrons in his seminal work [Anderson 1958], in reality, they fundamentally interact. This made the unambiguous observation of Anderson localization of electrons challenging, as there is also another type of metal-insulator transition: Mott transition [Mott 1990]. This interaction-driven transition is due to electron interactions, and it has nothing to do with the disorder. Moreover, photons can have an extremely long coherence length, which implied that the localization of light could be observed at room temperature and in macroscopic samples, contrary to electrons, which have a small coherence length, even at low temperatures, due to interactions.

However, there is an important difference between light and electron waves that made the light localization less ideal. Let us, therefore, more closely compare the characteristics of the quantum and classical waves in disordered systems [John 1991, Sheng 2006[START_REF] Akkermans | [END_REF]].

The Schrödinger wave equation for an electron of mass m in a disordered lattice is given by

- 2 2m ∇ 2 + V (r) ψ(r) = Eψ(r). (III.15)
The random potential V (r) that electron feels in the disordered medium is basically a local scattering potential produced by the atomic impurities or defects in the lattice. Since the energy of electron E can be negative, the electron localization can be achieved not only by increasing the density of impurities, i.e. the disorder, but also by making E sufficiently small comparing to the random potential, such that the electron is trapped in a deep negative well of the random potential. In other words, since for quantum waves the disorder can be dominant in the limits of low E, Ioffe-Regel criterion [Eq.(III.14)] might be satisfied by minimizing the energy E.

Let us now consider an electromagnetic wave of frequency ω that propagates in a random medium whose dielectric constant is (r) = + δ . The disorder is described by the fluctuating part δ of the dielectric constant, while is the average value of (r). For simplicity, we will disregard the vector nature, i.e. the polarization of the electromagnetic wave and consider the scalar electric field described by the wave function ψ(r). Then, the analogue of Schrödinger equation Eq.(III.15) for classical waves is the Helmholtz wave equation

-∇ 2 ψ(r) -δ ω 2 c 2 ψ(r) = ω 2 c 2 ψ(r). (III.16)
The quantity ω 2 c 2 plays a role analogous to the electron energy E, while -δ ω 2 c 2 is the analogue of random potential V of Eq.(III.15) and it is responsible for the wave scattering. Since the dielectric constant (r) is real and positive everywhere (we assume a non-dissipative medium), the energy quantity ω 2 c 2 is always positive. Hence, in contrast to electronic systems, it is impossible to have a bound state of photon in a negative potential well. Moreover, having (r) > 0 everywhere means that the photon energy has to be greater than the highest of the potential barriers. Note that for classical waves, both the energy and disorder potential quantities depend on the wave frequency as ∝ ω 2 . Therefore, by decreasing the photon energy, i.e. its frequency ω, the disorder potential δ ω 2 c 2 is decreasing in the same way. Contrary to electron waves, for which localization is enhanced by lowering the electron energy, lowering the photon energy leads to a complete disappearance of scattering, i.e. it destroys the disorder.

Anderson localization of light, hence, cannot take place by minimizing the photon energy, i.e. its frequency.

III.2.2. Scattering regimes for localization of light

Classical scatterers

Despite the problem of low frequencies for Anderson localization of light, in 1984, S.

John [John 1984] suggested frequency regime in which light could be localized in strongly disordered 3D systems. The idea is based on the following. When an electromagnetic wave propagates through the medium containing disorder made of randomly positioned dielectric particles of typical size a, there are mainly three different single scattering regimes [START_REF] Born | [END_REF], Kerker 1969, Bohren 1983, van de Hulst 1980]. In the regime where the wavelength is much larger than the size of the scatterers, i.e. the low energy limit, λ a (Rayleigh scattering), the scattering cross-section is σ s ∝ a 6 /λ 4 , while in the opposite limit λ a (geometrical optics) the scattering cross-section is σ s ≈ 2πa 2 λ 2 . To determine the order of magnitude of the minimal achievable mean-free path ∝ 1/(ρσ s ) in these two regimes, we have to consider the highest density of the scatterers ρ max ∼ 1/a 3 which then yields:

k min ∼      λ 3 /a 3 1 Rayleigh scattering with λ a, σ s ∝ a 6 /λ 4 a/λ 1 geometrical optics with λ a, σ s ≈ 2πa 2 .
(III.17)

Since in these two regimes the mean free path is substantially larger than the wavelength, the Ioffe-Regel criterion [Eq. (III.14)] cannot be fulfilled: the diffusive to localization transition is not achievable even for high densities (close packing of scatterers) as long as λ a or λ a. In other words, in these weak-scattering regimes the photon eigenstates are extended and the wave transport is diffusive.

However, there is an intermediate regime where the wavelength is comparable to the size of the scatterers λ ∼ a (Mie scattering regime). In this regime, due to so-called Mie resonances in the scattering from dielectric spheres, the scattering cross-section can become quite large, σ s ∼ a 2 . Then the mean free path (for ρ max ∼ 1/a 3 ) is

k min ∼ a λ ∼ 1 Mie scattering with λ ∼ a, σ s ∼ a 2 . (III.18)
This suggests that, in the regime λ ∼ a, a sufficiently strongly scattering to localize light is The scattering mean free path (the typical distance between two consecutive scattering events) in the independent scattering approximation [Lagendijk 1996] is given by = (ρσ s ) -1 where ρ is the density of the scatterers and σ s is the total scattering cross section of each individual scatterer. Note that, as the Mie scattering is anisotropic (the scattering is preferential in forward direction), the relevant mean free path in regime λ ∼ a is the transport mean free path * which corresponds to the typical distance over which the wave has lost the memory of its initial direction. The connection between the scattering and transport mean free path is given by * = 1 -cos θ , where cos θ represents the average cosine of the scattering angle θ. In the case of isotropic scattering, such as Rayleigh scattering, * = .

The Ioffe-Regel criterion [Eq. (III.14)] then reads k * 1. Now, let us assume a disordered medium as dielectric particles (spheres) of refractive index n 1 randomly distributed in a surrounding medium of refractive index n 2 . We can rewrite the density of those scattering particles as ρ = f /V , where f is the volume filling fraction of the spheres and

V p = 4πa 3 /3
is the volume of a sphere. Thus, * can be reduced by using a sample of closely packed scatterers and by maximizing the scattering cross section. Although one might think that maximizing the density of the scatterers might lead to small enough * for localization to occur, that does not hold. If the density of scatterers is too high, the entire medium might get inverted: the scatterers take the place of the surrounding medium and vice versa, making the scattering rather weaker than stronger. Thus, not only the intermediate frequency range is required to achieve very small * , but also the intermediate volume fraction of spheres, up to f ∼ 50%. Moreover, if the filling fraction of the scatterers is too high, interparticle correlations and near-field effects might arise [Lagendijk 1996] leading to the increase of the mean free path [Fraden 1990, McNeil 2000]. On the other hand, the mean free path can be reduced by maximizing the scattering cross section. The scattering cross section of Mie resonances is larger for higher refractive index contrast between scatterers and surrounding medium m = n 1 /n 2 . Therefore, to increase the disorder strength, i.e. to minimize , a high refractive index for the dielectric spheres is needed [Sheng 1986]. As is going to be further discussed in the following subsection, apart from generating a disordered sample with an optimum in volume filling fraction, the limitation in high refractive indices of materials creates a difficulty in obtaining strong enough scattering such that the localization condition k * < 1 can be fulfilled.

Resonant scatterers

In the Mie scattering regime, producing such strongly scattering sample for light localization to arise, is very challenging since one has to have control over many parameters (wave frequency and size of the dialectric scatterer, filling fraction of scatterers, index refraction contrast) to sufficiently minimize the mean free path and due to limitations in optimal values of those parameters.

So far, we have been discussing scattering of light on classical dielectric objects. Another way of reaching the localization regime for light waves may consist in using resonant point scatterers, such as atoms, as proposed by D. Sornette [Sornette 1988]. Note that while in the case of Mie scatterers, their Mie resonances are of geometrical origin, the resonant scattering on atoms is due to their internal quantum structure, i.e. internal resonances. As the wavelength of light is typically much larger than the size of an atom, λ a, the spatial variations of electromagnetic field over the size of the atom can be neglected, i.e. one can assume that the electric field E is uniform over the atom. In this long wavelength limit, the atom-photon interaction is simply given by d • E, where d is the atomic dipole moment.

Hence, atoms can indeed be considered as point dipoles, i.e. pointlike scatterers in the long wavelength limit. Moreover, atoms have well-defined transition lines, i.e. very narrow resonances. For the single scattering of near-resonant and weak monochromatic light (such that there is no inelastic scattering) on an atom in the long wavelength limit, the scattering cross section is given by [ de Vries 1998, van Rossum 1999]

σ s = σ 0 1 + 4∆ 2 /Γ 2 ; σ 0 =                3λ 2 2π vectorial light λ 2 π scalar light, (III.19)
where ∆ = ω -ω 0 is the detuning of the light frequency ω from the atomic transition frequency ω 0 and Γ is the natural width of the atomic transition. Since σ 0 ∝ λ 2 , the scattering cross section of light can become very large near resonance. This means that in order to reach the transition criterion ∼ λ/2π, the large density of the scatterers, as in Mie scattering regime for classical scatterers, is not needed. For example, for λ = 780 nm and resonant scattering with σ 0 = 3λ 2 /2π, to achieve the critical amount of scattering k c ∼ 1, a density ∼ 10 14 atom/cm 3 is needed. Another advantage in using atoms is that they are identical scatterers:

they have the same resonant frequency and linewidth of the transition (monodisperse sample).

Cold atoms as resonant point scatterers are promising candidates to study localization of light.

Cold atomic gases are well controlled and characterized systems, and their interaction with external fields is well understood. They are favored over hot atomic vapors since in hot atomic vapors substantial atomic motion would cause significant Doppler and collisional broadening, dephasing and frequency redistribution which reduce the scattering cross section. With advanced laser cooling and trapping techniques [Metcalf 1999], it could be possible to reach sufficient densities of cold atomic gases ∼ 10 14 atom/cm 3 such that Ioffe-Regel criterion is fulfilled. Moreover for weak enough electric field (linear-optics regime), the inelastic scattering and absorption of photons is absent. Although cold atoms seem to be promising candidates for studying coherence effects and achieving light localization, there still might be some phase breaking mechanisms due to residual atomic motion which affect coherent transport properties. Such mechanism is the frequency redistribution due to accumulated Doppler and recoil frequency shifts, which will dictate how many times photon can be scattered before being substantially shifted off resonance [Labeyrie 2006]. There is another major obstacle in light localization with resonant point scatterers: near-field coupling that will be discussed in the subsections below. For a further discussion of localization of light with cold atoms see

Refs. [Kaiser 2000[START_REF] Kaiser | [END_REF], Havey 2009, Kaiser 2009].

III.2.3. Towards the observation of Anderson localization of light in 3D

disordered systems

Coherent backscattering of light

In parallel with the theoretical suggestions of strong localization of light [John 1984, Anderson 1985], the direct experimental observation of coherent backscattering (CBS) of light waves was achieved: a phenomenon often considered to be a precursor of Anderson (strong) localization. It was observed by three different groups around the same time: A. Ishimaru et al. [Kuga 1984, Tsang 1984]; van Albada and Lagendijk [START_REF] Albada | [END_REF]]; Wolf and

Maret [START_REF] Wolf | [END_REF]]. As was discussed in Sec. III.1.1, the multiple coherent-wave scattering in a disordered medium is an underlying mechanism of coherent backscattering, but this phenomena does not require a strong disorder as the Anderson localization in 3D does.

Hence, CBS is considered to be an indicator of weak localization, in which the wave diffusion is lowered due to interferences on time-reversed scattering paths [Akkermans 1986]. The observation of this precursor phenomena implied that the Anderson (strong) localization of light should be possible, and it prompted a fuelled quest of searching for Anderson localization of light.

The mentioned observations of coherent backscattering in the 1980's were achieved with a disordered medium consisting of classical particles, i.e. dielectric spheres with a ∼ λ embedded in a surrounding medium. Observations of CBS phenomena has also been achieved in cold atoms. First experiments reporting observation of coherent backscattering of light were done using a gas of cold Rubidium atoms ( 85 Rb) on the F = 3 → F = 4 transition [Labeyrie 1999, Labeyrie 2000]. However, the reported enhancement factor of the intensity in the backward direction is only ∼ 1.1, much smaller than that observed with classical scatterers (for the systems where there is no decoherence mechanisms, the theoretically predicted enhancement factor is 2, see Sec. III.1.1). That indicated that the interference effects along multiple-scattering paths are more prominent for classical scatterers than for atoms. It has soon been realized that the internal atomic degrees of freedom have a significant influence on interference properties of light scattering: the Zeeman degeneracies of the internal atomic structure in 85 Rb (especially the degeneracy of the ground state) is responsible for the reduction of CBS enhancement [Jonckheere 2000[START_REF] Müller | [END_REF], Müller 2002, Müller 2003].

Because of the Zeeman degeneracy of the ground state, the probability amplitudes of direct and reversed path are generally not equal. Therefore, the imbalance of the amplitudes of the time-reversed paths leads to the partial loss of coherence and reduces the interference3 [Kupriyanov 2003, Labeyrie 2003]. The role of the atomic level degeneracy on coherent backscattering has also been interpreted in terms of dephasing induced by the internal Zeeman degrees of freedom [Akkermans 2002[START_REF] Akkermans | [END_REF]]. It has been shown that this phase breaking mechanism arises only when the ground state is degenerate (J = 0), i.e.

it is absent for non-degenerate ground state J = 0. This has been indeed experimentally confirmed by the observations of coherent backscattering in a gas of cold strontium atoms on J = 0 → J = 1 resonance transition where the factor of two enhancement of the intensity in the backscattered direction was measured [Bidel 2002,Wilkowski 2003,Wilkowski 2004]. Since the Zeeman degeneracy of the atomic transition affects the coherent transport and reduces interference, it is highly advisable to use a non-degenerate transition, such as J = 0 → J = 1, in the quest for reaching the strong localization regime where interferences play a crucial role.

Moreover, it has been observed that, when the transition is degenerate, applying a magnetic field restores the twofold enhancement of the backscattered intensity [Sigwarth 2004]. This came as a surprise, since it was demonstrated with classical scatterers that adding an external magnetic field breaks time-reversal symmetry (due to Faraday rotation effect [Martinez 1994]) and causes amplitude imbalance of the direct and reversed paths leading to diminished interference effects and reduced CBS enhancement [Lenke 2000a]. However, the role of the magnetic field is different with atoms where the amplitude imbalance already exists due to the Zeeman degeneracy of atomic transition. Applying a magnetic field in atomic system lifts the Zeeman degeneracy of the states and for large enough magnetic field, when the Zeeman sublevels are well split ( Γ), the atoms behave as an effective two-level system with no internal degeneracy and hence, the balance of the time-reversed path amplitudes is restored. The effect that applied strong magnetic field has on light scattering from atoms with degenerate internal structure can be associated with the enhancement of the coherence length. It has been also suggested that similar effect of restoring the interference might arise by optical pumping [Kupriyanov 2004]. For further studies on CBS in cold atomic gases, including effects of a thermal motion and non-linear medium, see Refs. [Labeyrie 2004, Chanelière 2004b, Müller 2005, Wellens 2005, Wellens 2006, Labeyrie 2006], while for an overview Refs. [Chanelière 2004a, Havey 2005, Kupriyanov 2006, Labeyrie 2008, Aegerter 2009].

Note that up to date, apart from experiments on coherent backscattering, there has been no reported experiments on Anderson localization of light with cold atomic gases, only with classical particles. Thus, the rest of the discussion in this subsection concerns experiments realized with classical scatterers.

Attempts in observation of strong localization of light

Over the years, there have been many reported experimental observations of Anderson localization of light in low dimensional systems (1D and 2D) [START_REF] Topolancik | Experimental Observation of Strong Photon Localization in Disordered Photonic Crystal Waveguides[END_REF], Lahini 2008, Sapienza 2010, Schwartz 2007, Riboli 2011]. However, despite numerous studies and experimental efforts to observe Anderson localization of light in a three-dimensional disordered system, its unequivocal observation is still lacking up to this date [Skipetrov 2016b].

There have been a couple of reported observations of Anderson localization of light in 3D, but later they got revoked. In 1997, Wiersma et al [Wiersma 1997] claimed this observation with near-infrared light (λ = 1064 nm). They used semiconductor gallium arsenide (GaAs) powders as a strongly scattering medium. The pure GaAs is characterized by a high refractive index (n ≈ 3.5) and extremely small absorption. By varying the grinding time of GaAs crystals, they were able to produce samples of different average particle size (with a finite degree of polydispersity). Their search of Anderson localization of light was based on theoretical predictions on the static transmission of light, i.e. transmission as a function of the sample thickness. In the diffusive regime, the transmission is described by Ohm's law, which predicts the decrease of transmission with the sample thickness as T ∝ 1/L. Near the mobility edge, scaling theory of localization predicts that the diffusion coefficient at the phase transition decreases as 1/L [Abrahams 1979] which leads to a quadratic scaling of the transmission T ∝ 1/L 2 [John 1984, Anderson 1985]. Moreover, deep in the localized regime, where the diffusion coefficient vanishes, the transmission will decrease exponentially as T ∝ exp(-L/ξ) [John 1984, Anderson 1985]. By decreasing the size of the particles, and hence increasing the scattering strength, Wiersma et al observed the change in transmission behaviour from T ∝ 1/L to ∝ 1/L 2 and finally to ∝ exp(-L/ξ). As this observed exponential decrease of the transmission was theoretically predicted for the localized regime, they interpreted it as an evidence for Anderson localization of light waves. However, the observations were questioned [Scheffold 1999, Wiersma 1999

III.2.4. Challenges in observing light localization in 3D

Let us now summarize and discuss some of the main difficulties in achieving and observing strong localization of light with classical scatterers and atoms.

Localization signature

Firstly, one needs to know how to look for the photon localization, i.

Disorder strength

Furthermore, the question is where to look for the light localization, i.e., what disorder (or scattering) strength is sufficient for the states to become localized. The measure of the disorder strength is given by the product of the mean-free path of the photon propagating in the scattering medium and the wavevector k. According to Ioffe-Regel criterion, localization takes place when k < k c ∼ 1. However, k c ∼ 1 represents the qualitative measure of the critical value of the mean free path for the phase transition from extended to localized states; the exact critical value k c is not known and it is very conceivable that it is unique for the specific type of wave and the disordered sample. Related to this is also the challenge in achieving strong enough disorder, i.e. in minimizing sufficiently the mean free path for localization to occur. The mean-free path, as discussed in the previous section, depends on many parameters (density of the scatterers, wavelength, as well as size of the particle and index of refraction contrast in the case of classical scattering particles). Finding an optimum in these parameters to minimize the mean-free path sufficiently for the phase transition to take place is not an easy task, especially since an underlying physical mechanisms might arise that would have an opposite effect, i.e. that would lead to the increase of the mean-free path instead of further minimizing it. An example for that is the one of the tightly packed Mie scatterers where the further increase of their density is not a simple way of decreasing the mean free path since at very high densities the correlations between the scatterers and near-field interactions arise, which prevent the further reduction of the mean free path [Lagendijk 1996].

For resonant point scatterers like atoms, such very high densities as for Mie scatterers are not required to have k ∼ 1, due to their huge near-resonance scattering cross section. However, near-field interactions have an important effect on atoms as well, as it will be discussed below.

A straightforward way of increasing the scattering strength, i.e. of reducing the mean-free path for photons scattered from Mie scatterers, would be to increase the refractive index contrast between the scattering particles and surrounding medium. However, the range of refractive indices of typical optical materials is not very broad: from n ≈ 1 in the air up to n ≈ 4 for germanium [Andreoli 2021]. Therefore, the narrow range of refractive indices does not allow to lower the photon mean-free path immensely. For example, in the experiments mentioned above with titania (n ≈ 2.7) and gallium arsenide (n ≈ 3.5) powders compressed to high densities, the minimal mean-free path that was acquired was down to kl ≈ 3. Even in the failed attempt to observe the light localization using the germanium powder (n ≈ 4), the minimal obtained mean-free path was kl ≈ 3 [Gómez Rivas 2001]. However, this does not mean that the experiments on light localization with classical scatterers of high refractive index reached a dead end. Currently unreached low k could be obtained by, for example, progress in disordered sample production such that the degree of polydispersity is lowered.

Near-field interactions

In the case of a disordered medium consisting of classical scatterers, many studies have shown that upon increasing the density, the positional correlations among the scatterers start to have a significant influence on photon transport: spatial correlation among scatterers yields phase correlation between scattered waves and weakens the effective scattering cross section below that of an individual scatterer, and therefore increases the mean-free path4 [Fraden 1990, Saulnier 1990]. Apart from short range correlations, in the case of very tightly packed dielectric particles, such that the particles are in the near field of each other, nearfield effects arise which reduce the scattering strength and increase the photon mean-free path [McNeil 2000,Peng 2007,Sapienza 2007,Liew 2011,Rezvani Naraghi 2015]. As it has been discussed in Ref. [Rezvani Naraghi 2015], these near-field effects originate from short-range evanescent wave coupling of adjacent scatterers (non-radiative energy transfer between the scatterers): the near-field coupling between particles can be seen as an opening of new transmission, optically connected, channels. Arising of additional transport channels of non-radiative energy transfer due to near-field interaction among the scatterers, effectively increases the photon mean-free path. It has been showed that the near-field effects can be taken into account in the photon transport description by considering that the scattering particle is surrounded by an effective homogeneous dielectric medium [Aubry 2017, Busch 1996].

The detrimental character of near-field interactions on Anderson localization of light has been reported in the case of resonant point scatterers [Skipetrov 2014]. S. E. Skipetrov and I.

M. Sokolov have numerically studied the regimes of light transport considering scalar and vectorial near-resonant light waves propagating in a random discrete 3D system of point-like particles with resonance transition J = 0 → J = 1. They based their analysis on applying the scaling theory, discussed in Sec. III.1.2.1, on the eigenvalues of the Green's matrix (analogous to the one introduced in the previous chapter) from which they obtained the Thouless number g Th . They have shown that with the scalar light waves, for which the near-field dipole-dipole interaction term is absent, the phase transition from extended to localized states, i.e., mobility edges occurs at g Th ≈ 1 [Skipetrov 2016a]. On the other hand, when they assumed the vectorial nature of light, for which the near-field dipole-dipole interaction term is present, they did not identify any mobility edge corresponding to the phase transition (and with that they did not obtained g Th < 1 or the β(g Th ) < 0 branch of the scaling function [Eq. III.7])

suggesting that Anderson localization of vectorial light cannot be achieved [Skipetrov 2014].

Similar results were also obtained by [Bellando 2014]. Later on, they demonstrated that by applying a strong magnetic field, which significantly splits the Zeeman sublevels, the Anderson localization of vectorial light is restored, as the strong magnetic field reduces the contribution of the near-field dipole-dipole interaction terms [START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF]]. The fact that the photon localization is not present for resonant point scatterers in the vectorial wave model is attributed to the near field dipole-dipole interactions. The reason is the following. The propagating electromagnetic waves are transversal by nature, however the excitations from one atom to another are not only transferred by propagating waves, but also non-radiatively through near field dipole-dipole interactions which arise due to the longitudinal component of the radiated electromagnetic field. These near-field interactions open additional non-radiative energy transport channels. The importance of near-field interactions in the excitation transfer increases with increasing density ρ of the scatterers: while they are negligible in dilute samples ρλ 3 1, non-radiative channels become essential for very small inter-particle distances, i.e.

high densities ρλ 3 1 [Nieuwenhuizen 1994]. Therefore, the opening of these non-radiative transport channels due to the near-field interaction prevents the diffusion coefficient to vanish and thus leads to the absence of light localization. That the near-field interactions are indeed preventing the light localization has been furthermore corroborated by considering elastic waves which have propagating longitudinal component and hence, near-field interaction and non-radiative transport are not induced. The results based on the scaling analysis showed the existence of the localized regime, similar to the scalar case for light waves [Skipetrov 2018b].

Moreover, Skipetrov and Sokolov presented the phase diagram for the scalar [Skipetrov 2018c] and vectorial models [Skipetrov 2018a], showing the critical light frequencies for which the phase transition occurs as a function of the density of the sample of the point-like particles.

They identified that the critical density at which scalar light waves exhibits the Anderson transition is ρ c λ 3 ≈ 20, while for the vectorial model with a strong magnetic field, the critical density is ρ c λ 3 ≈ 25.

As it was discussed above, the near-field effects can also impact the transport of light in Recent numerical work [Skipetrov 2016c], in the framework of cold atoms as resonant point scatterers, demonstrated that the late-time decay dynamics of the atomic system, with and without applying strong magnetic field, exhibits slower and faster decay, respectively.

This suggested that the localization signature could be in the prolonged decay of the scattered intensity at very late times, since magnetic field restores the strong localization of light in the system of resonant point scatterers [START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF]]. However, in our recent work [Cipris 2021b], that was described in chapter 2, we showed that near-field terms are detrimental as well for the longest-lived subradiant states. Since the magnetic field reduces the near-field interaction, it is questionable whether this prolonged decay stems from magneticfield resurrected subradiant states or localized states. Whether the localization signature could be found in the late-time decay dynamics is still an open question. Nevertheless, to distinguish Anderson localization from subradiance in the late-time decay dynamics, the analysis based on scaling theory [Abrahams 1979] including different sample sizes, atom number, detunings etc., would have to be done.

Another signature of light localization can be found in the speckle statistics, i.e. variance of the transmitted intensity as demonstrated in [Chabanov 2000] for microwaves in quasi-1D

III.3.1. Thouless number and conductance

In Sec. III.1.2.1 we have introduced the Thouless number g Th = ∆ω/δω, the ratio between the typical linewidth of an eigenmode ∆ω and the average eigenmode spacing δω, which can be interpreted as the average number of eigenmodes in a frequency width ∆ω. The linewidth of an eigenmode (or simply mode), which can be interpreted as the leakage rate of energy from the sample, is closely linked to the sensitivity of modes to changes at the boundary. When the wave is localized within the sample, the amplitude squared of the wave is exponentially small at the boundary and the mode is only weakly coupled to the environment. The mode lifetime is then long and its linewidth correspondingly narrow, so that, ∆ω < δω, i.e. g Th < 1. On the other hand, when the wave is diffusive, modes extend throughout the medium and they couple easily to their surroundings; the energy readily leaks out of the sample and levels are consequently short lived with linewidths greater than the typical spacing between neighbouring modes such that the modes overlap, ∆ω > δω, i.e.

g Th > 1. Thus, the Thouless number, which represents the spectral mode overlap, can be used to distinguish between extended and localized waves, with the localization threshold expected to take place at g Th = 1. However, g Th is not an easily measured quantity (it has been measured only recently [START_REF] Wang | [END_REF]). Nevertheless, since it characterizes transport properties, it can be related to easily measured transport quantities. In the case of electron transport, g Th has been equated to the dimensionless conductance [Eq. (III.

3)] via Einstein's formula [Thouless 1977]. The connection between the Thouless number and the conductance can also be made for classical waves, since the dimensionless conductance can be defined via the Landauer relation as g = a,b T ab = T [Landauer 1970], where T ab is the transmission coefficient between an incoming channel a and outgoing channel b and T is the sum of transmission coefficients over all incoming and outgoing channels, which we will refer to as transmittance. Thus, this relation, which expresses the equivalence between the dimensionless conductance and transmittance, provides an analogy between quantum and classical waves.

Moreover, for diffusive waves, g = T ∝ M L [START_REF] Akkermans | [END_REF], Lagendijk 1996], where is the scattering mean-free path, L is the size of the sample and M is the total number of propagation modes (transport channels) that the sample supports. This relation allows us to view the dimensionless conductance as being related to the sum of contributions of M independent scattering channels each giving on average a contribution of order /L.

Moreover, the interpretation of dimensionless conductance based on transmission eigenvalues was given in terms of open and closed channels [Imry 1986]: most transmission eigenvalues in a disordered system with L are exponentially small (closed channels), however a fraction /L of the total number M of transmission eigenvalues is of order unity (open channels).

Only the open channels M eff = M /L contribute effectively to the conductance. Thus, the dimensionless conductance can be interpreted as the number of open channels or the effective number of active propagating channels M ef f . As the mean-free path decreases for a given L, the effective number of propagating modes and hence, the conductance are decreasing.

As the conductance essentially represents a number of diffusive channels in the sample, the localization is expected to take place when its value falls below unity, which is in agreement with Thouless criterion for localization threshold g Th = 1. In the localized regime the photon transport is exponentially suppressed, as only the tails of the localized wave intensity can leave the sample. Hence, for localized waves g = T ∝ exp(-L/ξ) < 1 with ξ < L, where ξ is the localization length.

III.3.2. Speckle statistics: intensity distribution and intensity fluctuations

Let us now consider the intensity in the transmission when the disordered sample is illuminated by a laser beam. As a result of multiple scattering from disordered medium a complex interference pattern of intensity is produced consisting of many irregularly distributed bright and dark spots known as speckle [Goodman 2009] (see Fig. III.5). This speckle pattern is a result of interferences between light waves multiply scattered on different paths. Each configuration of scatterer's positions produces its own speckle pattern. We can write the intensity detected at point r, as I(r) = I + δI, where I is the average intensity and δI is the fluctuating part of the intensity. This fluctuating part, which is due to interferences, is in fact responsible for appearance of speckle spots, i.e. it is the intensity speckle pattern.

Deeply in the diffusive regime (g 1), under the assumption that the photon scattering paths are uncorrelated, the speckle intensity has an exponential distribution, which is known as Rayleigh distribution [Goodman 2009]

P ( Î) = exp(-Î), (III.20)
where Î = I/ I is the intensity normalized by its average value. The Rayleigh distribution obeys

σ 2 Î = δI 2 I 2 = 1, (III.21)
where

σ 2 Î = δI 2 I 2 = I 2 I 2 -1, (III.22)
is the variance of the normalized intensity Î with δI = I -I . The Eq. (III.21) is the wellknown Rayleigh law which states that the relative fluctuations of the intensity (normalized by its average value) are of order unity.

As mentioned, Rayleigh law assumes that the correlations between different scattering paths in the medium are absent. However, the interferences between random scattering paths in disordered medium lead to such correlations. As the conductance decreases, the interferences become more substantial yielding enhanced degree of correlations. These correlations are encoded in the intensity distribution, i.e. interferences modify the intensity distribution. By including the corrections due to correlations between different random paths in the medium, it has been shown that the intensity distribution is given by 5 [START_REF] Nieuwenhuizen | Intensity Distributions of Waves Transmitted through a Multiple Scattering Medium[END_REF] Rossum 1999]

P ( Î) = i∞ -i∞ dx πi K 0 (2 -x Î) exp (-Φ c (x)) , (III.23a) Φ c (x) = g 1 0 dy y log 1 + xy g + xy g , (III.23b)
with K 0 the modified Bessel function. Note that this intensity distribution function depends on the dimensionless conductance g. The Rayleigh distribution is recovered for g 1 (and as long as g Î). Moreover, note that this intensity distribution is derived considering the diffusive regime far from localized regime, as only the leading corrections due to correlations were considered, in terms of interferences of loopless scattering paths. The interferences on looped paths that are becoming more and more significant as the localization threshold is being approached, and especially in the localized regime, are expected to further modify the intensity distribution. However, this intensity distribution function characterized well the localized microwaves in quasi-1D disordered system [Chabanov 2000], as well as the localized regime for ultrasound in a 3D scattering medium [Hu 2008]. Hence, we are also going to turn to this intensity distribution function to extract the dimensionless conductance from our numerically obtained intensity distributions.

Furthermore, as the interferences between random scattering paths in a disordered medium modify the intensity statistics, consequently the variance of the normalized intensity [Eq.

(III.22)] will exhibit deviations from Rayleigh law, i.e. from unity [Eq. (III.21)]. For a disordered waveguide, the variance of the normalized intensity, when taking into account the correlation corrections, is given by [START_REF] Akkermans | [END_REF] 5 For an incoming beam with a Gaussian profile.

σ 2 Î = 1 + 4 3g + 2 15g 2 .
(III.24)

The first term expresses the Rayleigh law and corresponds to uncorrelated channels, while the following two terms describe correlations between channels. Note that Eq. (III.24) is essentially derived for the diffusive regime when taking into account the corrections due to interferences on loopless scattering paths. Nevertheless, it give us an insight into how the variance is modified by the interferences: it implies that the relative intensity fluctuations are of order unity deeply in the diffusive regime (g 1), while they become large in the localized regime (g < 1). In addition, the first term of Eq. (III.24) is known as the short-range correlation, while the second and third term as long-and infinite-range correlations, respectively. In fact, these short-, long-and infinite-range correlation terms are often defined as the leading terms in the variance of the intensity, total transmitted intensity and conductance, respectively, as each of the correlation terms have higher order corrections in terms of 1/g and 1/g 2 [van Rossum 1999]. For example, the weak-localization correction would yield 1/g contribution to the short-range correlation. Such corrections were studied in [Garcia 1993]. Moreover, the strong deviations of the normalized intensity variance from unity were observed in microwave quasi-1D experiments in the localized regime, indicating that the variance of the intensity is a mean to distinguish between extended-wave and localized regimes [Chabanov 2000].

III.3.3. Numerical model

As in the previous chapter, we use the coupled-dipole model, scalar and vectorial, to obtain

dβ j (t) dt = i∆ 0 - Γ 0 2 β j (t) - dE L exp(ik L • r j ) - Γ 0 2 m =j G(r jm )β m (t) , (III.25a) with G(r) = exp(ik 0 r) ik 0 r , (III.25b)
while the response of the system in the vectorial model is given by the following set of coupled-dipole equations for amplitudes β ζ j of atomic dipole j on transition to m e = ζ:

dβ ζ j (t) dt = i∆ 0 - Γ 0 2 β ζ j (t) -i d ê * ζ • E L exp(ik L • r j ) - Γ 0 2 m =j η G ζ,η (r jm )β η m (t), (III.26a) with G ζ,η (r) = 3 2 exp(ik 0 r) ik 0 r δ ζ,η -rζ r * η + δ ζ,η -3r ζ r * η i k 0 r - 1 (k 0 r) 2 .
(III.26b)

In Let us recall that the first term of (scalar and vectorial) coupled dipole equations describes the dynamics of an individual atom, the second one describes the interaction of an atomic dipole with the incident electric field, while the third term is the atom-atom (i.e. dipole-dipole) coupling term which describes the interaction of the atomic dipole with the electric field radiated from other dipoles. This third term takes into account that each atomic dipole is affected by other dipoles via the radiation field that they produce: coupling between atom i and atom j is mediated by the excitation transfer from atom i to atom j. Therefore, this third term accounts for multiple scattering and collective effects. The dipole-dipole interaction is encoded in the dyadic Green's function: the scalar Green's function G(r ij ≡ r ir j )

I(R) ∝ E L (R) exp(ik L R) -i dk 3 0 4π 0 exp(ik 0 R) ik 0 R j β j,st exp(-ik 0 nr j ) 2 (III.27)
in the scalar couple-dipole model, and

I(R) ∝ E L (R) exp(ik L R) -i dk 3 0 4π 0 exp(ik 0 R) ik 0 R j ζ,η β η j,st (δ ζ,η -nζ n * η ) exp(-ik 0 n • r j ) 2 (
E L (R) = E L (X, Y, Z) = E 0 w 0 w(Z) exp - X 2 + Y 2 w(Z) 2 exp i k L X 2 + Y 2 2R(Z) -arctan Z z R , (III.29)
where

E 0 = E L (0, 0, 0), w(Z) = w 0 1 + (Z/z R ) 2 is the beam radius at position Z, w 0 = w(0) is the beam waist corresponding to the minimum beam radius, R(Z) = Z [1 + (z R /Z) 2 ]
is the radius of curvature of wavefronts and z R = k L w 2 0 /2 is the Rayleigh range which determines the length over which the beam can propagate without diverging substantially.

Note that in the scalar model the incident electric field is given by E

in (R) = E L (R) exp(ik L R),
while in the vectorial model we take into account the polarization of the electric field ˆ L so that the incident electric field is

E in (R) = E L (R) exp(ik L R) = E L (R)ˆ L exp(ik L R).
We consider that the atoms are randomly distributed inside of the cube of side length L or inside of the cylinder of diameter L and height L. Furthermore, we consider a uniform density distribution, so that the density of atoms is ρ = N/L 3 and ρ = 4N/(πL 3 ) for cubic and cylindrical atomic sample, respectively. Note that in the case of cylindrical sample, we assume that the axis of the cylinder is along the propagation axis of the electric field.

Moreover, contrary to the numerical study in the previous chapter, here we do not impose any minimal distance, i.e. exclusion volume on atom positions. Furthermore, we compute the transmitted intensity around the forward direction, at a given angle ({θ obs , φ obs } = n) from the laser propagation axis, such that we avoid the Mie-scattering lobe (∝ 1/k 0 L; which is due to the finite size L of the disordered medium, i.e. the atomic sample). Moreover, to study the statistics of the transmitted light, we compute the intensity at the given observation point for many spatial configurations (i.e. realizations) of atomic positions. Note that computing the intensity for different positional configurations implies that only the atomic positions, within the given sample geometry and uniform density, are different from one realization to another, while all the other parameters of the system are the same.

III.3.4. Signature in the scalar coupled-dipole model

Let us first discuss the scalar-wave case where the internal Zeeman structure, as well as the light polarization, are discarded [Eq. (III.25a)]. The atoms are uniformly distributed within a cube of a side length k 0 L = 32.4. The transmitted intensity is obtained from Eq. (III.27) for the observation angles (θ obs , φ obs ) = (75 • , 0 • ) from the laser axis ẑ in the far-field at distance R = 250L from the sample. Moreover, the radiated intensity is obtained for many realizations of random atomic positions.

In the very dilute limit of an atomic sample, i.e. the limit of a weakly disordered sample corresponding to the diffusive regime, where the photon scattering paths are uncorrelated, the transmitted intensity has a probability distribution function that obeys Rayleigh law and its variance σ 2

Î is equal to one. Such behaviour is illustrated in Fig. III.6(a) for an atomic sample with a low density, below the localization threshold (ρ < ρ c ≈ 20/λ 3 ). However, for higher densities above the localization threshold, deviations from Rayleigh law appear, as can be observed in Fig. III.6(b). Here, the numerically obtained intensity distribution is well described by the intensity distribution given by Eq. (III.23) which takes into account correlations. The dimensionless conductance g is obtained from fitting the numerically obtained intensity probability distribution to Eqs. (III.23). In the diffusive limit, i. [Skipetrov 2018c]. Indeed, when the Lorentz-Lorenz shift (local-field effect) is accounted for in the evaluation of the mean-free path, it leads to the following critical laser detuning ∆ c to meet the criterion k 0 < α [Kaiser 2000, Kaiser 2009]: The very good agreement between the results of the scaling analysis [Skipetrov 2018c] and the present ones suggests that the intensity variance is a suitable observable of the localization transition, and that the changes in the conductance g are indeed associated with that transition. Note that the variance of the intensity falls back to unity within the V-shaped area characterized by large intensity variance. This implies that the strong deviation of intensity variance is an indicator of the transition from diffusive to localized regime, rather than the localized regime itself (this will be discussed more in Sec. III.4, together with the fact that in the localized regime far from the transition we also observe an increase of the dimensionless conductance).

∆ c = ρλ 3 8π 2 ± 1 2 3α ρλ 3 4π
Moreover, to circumvent the possible role of finite-size effects, it has been checked that our results hold both at fixed system size kL while varying N , and at fixed atom number N while varying kL. Similar results were as well obtained by using a cylindrical geometry of the atomic sample, instead of cubic [Cottier 2019a].

Note that the results presented in this section ( Î ≈ 1 for all considered detuning and densities: we do not obtain significant deviations of intensity variance from unity, as is the case with the scalar model (Fig.

III.7

), which indicates that there is no strong localization of light in 3D with the vectorial coupled-dipole model. This is in agreement with theoretical works based on the scaling analysis [Skipetrov 2014, Bellando 2014], where it has been shown that the light localization is absent with the vectorial model due to the near-field dipole-dipole interaction, which opens up additional non-radiative photon transport channels [Skipetrov 2014[START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF] 1 atoms. For high enough densities in the scalar model, in the bottom of the distribution a tail of modes of very long lifetimes with large IPR emerges which corresponds to the localized modes. Localized modes with large IPR and very long lifetimes that appear only at sufficiently high densities are to be discriminated from subradiant pairs which are present at any density. On the contrary, in the vectorial model such localized states do not appear for higher densities in the bottom of distribution, signaling the absence of light localization [Bellando 2014, Skipetrov 2014].

As already mentioned, it has been shown that applying a strong magnetic field restores strong localization of light within the vectorial model [START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF]]. Hence, we now focus on intensity fluctuations with the vectorial coupled-dipole model when the magnetic field is applied. We consider a static and spatially uniform magnetic field aligned with the z axis, B = B ẑ. Having atoms subjected to an external magnetic field lifts the degeneracy of Zeeman substates, such that the resonant frequency of substates of excited state |J e = 1, m e = 0, ±1 is ω me = ω 0 + m e ∆ B , where ∆ B = g e Bµ B / is the Zeeman shift of substates, with g e the Lande factor of excited state and µ B the Bohr magneton. The set of vectorial couple-dipole equations with the applied external magnetic field are the following For m e = 0, the intensity is detected at (θ obs , φ obs ) = (15 • , 90 • ) 8 and we average over many realizations of atomic positions. The obtained variance of the transmitted intensity σ 2 Î as a function of the atomic density ρλ 3 and laser detuning ∆ me = ω L -ω me for m e = 0, ±1 is shown on Fig. III.12. For m e = ±1, we observe significant deviations of intensity variance from unity for large enough densities (ρλ 3 25), similarly to observations in the scalar model. However, for m e = 0 we do not observe such strong deviations. In addition, the strong intensity fluctuations observed for m e = ±1 appear in the narrow frequency range for positive laser detunings ∆ me ∼ 1Γ 0 (V-shaped area of increased intensity fluctuations, similar to the scalar case). This corresponds to our observations made on the eigenvalue distribution discussed above. Furthermore, our results are in agreement with the theoretical work based on scaling analysis of the Thouless number obtained from eigenvalues, where it has been demonstrated that light localizes for m e = ±1, but not for m e = 0 [START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF]]. In addition, the set of parameters (ρλ 3 , ∆ me ) for which we observe large intensity fluctuations is in agreement with the localized regime identified by the eigenvalue scaling analysis (see phase diagram on Fig. 3 of Ref. [Skipetrov 2018a]). The overall good agreement of our intensity variance observations in the vectorial model without and with applied strong magnetic field (Fig. III.8 and III.12) with the scaling analysis [Bellando 2014, Skipetrov 2014[START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF], Skipetrov 2018a] furthermore corroborates that the strong deviations of intensity variance from unity, i.e. strong intensity fluctuations, is a suitable observable for detecting the strong localization of light.

dβ ζ j (t) dt = i(∆ 0 -m e ∆ B δ ζ,me ) - Γ 0 2 β ζ j (t)-i d ê * ζ •E L exp(ik L •r j )- Γ 0 2 m =j η G ζ,η (r jm )β η m (t), (
G η,ζ j,m =        -im e ∆ B δ ζ,me - Γ 0 2 for j = m, ζ = η - Γ 0 2 G η,ζ (r jm ) for j = m; G η,ζ ( r
Let us now discuss why light localizes for m e = ±1, but not for m e = 0. Firstly, let us make some additional observations. We have discussed that the strong magnetic field splits the atomic system into three decoupled scattering subsystems. This can indeed be seen on Fig. field is applied and therefore the substates are degenerate. In that case, the amplitudes of each m e = 0, ±1 have approximately equal non-vanishing values, as the photon scattered from one of the transitions m e = 0, ±1 can be rescattered from any transition m e = 0, ±1 due to the degeneracy of sublevels. However, when the strong magnetic field is applied 8 To drive the m e = 0 transition, we need the laser beam with its linear polarization parallel with the applied magnetic field, i.e. along the direction of quantization axis ẑ. Then, we set the laser propagation direction and cylinder axis along ŷ direction. We want the observation angle relative to the laser propagation direction to be the same as for m e = ±1: (θ obs , φ obs ) = (15 • , 90 • ) for m e = 0 with k L ŷ corresponds to θ obs = 75 • and any φ obs angle for m e = ±1 with k L ẑ. Note that θ obs is the polar angle defined from positive z-axis, while φ obs is the azimuthal angle in the xy-plane defined from positive x-axis. ]. This indeed shows that the atoms subjected to strong magnetic field effectively behave as a non-degenerate two-level system. Thus, we can rewrite the vectorial coupled-dipole equations given by (III.31), with the corresponding dyadic Green matrix given by (III.26b), as a set of N coupled equations of each decoupled (transition-wise) subsystem associated with corresponding m e as: 3 for m e = ±1

dβ j (t) dt = i∆ me - Γ 0 2 β j (t) -i dE L exp(ik L • r j ) - Γ 0 2 m =j G (r
G (r) = 3 2i exp(ik 0 r)×            1 2 [1 + cos 2 θ] 1 k 0 r - 1 2 [1 -3 cos 2 θ] i (k 0 r) 2 - 1 (k 0 r)
[1 -cos 2 θ] 3 for m e = 0, (III.34)

1 k 0 r + [1 -3 cos 2 θ] i (k 0 r) 2 - 1 (k 0 r)
where θ is the angle between r and z-axis9 . We can see that the Green's matrix G for is applied [Bellando 2014, Skipetrov 2014[START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF], Skipetrov 2016a]. Furthermore, the set of parameters (ρλ 3 , ∆ 0,me ) for which we observe strong intensity fluctuations, for both scalar and vectorial models, are in an excellent agreement with the identified mobility edges and localized regime from the scaling analysis (compare our Figs. III.7 and III.12 with Figs. 3 of Ref. [Skipetrov 2018c] and [Skipetrov 2018a], respectively). This striking similarity between our results and the scaling-analysis based results, regarding the set of parameters for which we observe strong intensity fluctuations, furthermore corroborates that the significant deviations of intensity variance from unity are to be attributed to the phenomenon of light localization, and hence confirming that the variance of the intensity is a suitable observable of light localization.

There are, however, some open questions. One thing that has not been discussed yet is [van Rossum 1999]). Although this intensity probability function was able to characterize the phase transition in our system (as well as the localized regime of microwaves in quasi-1D sample [Chabanov 2000] and of ultrasound in 3D [Hu 2008]), it is plausible that deeply in the localized regime, where closed loops are expected to modify the intensity distribution, this function is not suitable to properly characterize the transport properties of the system. Moreover, in addition to localized modes in our system there are also other collective modes, like superradiant and subradiant modes. The question is whether there are correlations between localized modes and these collective modes, especially the subradiant ones, which as well essentially cause the light to be trapped inside of the sample, and how that affects the localization and intensity statistics. In fact, recently the intensity statistics at the mobility edges has been studied inside of the scattering medium [Skipetrov 2019], contrary to our case where we studied the intensity statistics outside of the disordered medium. Nevertheless, observations similar to ours have been made at the mobility edges:

in the scalar case the intensity variance at the mobility edge corresponding to the larger laser detuning is much larger than the intensity variance at the other mobility edge [Fig.

III.7(a)], while in the vectorial case with the applied magnetic field the intensity variance at the two mobility edges is similar [Fig. III.12]. The suggestion was made that this difference of intensity fluctuations at the mobility edges could be attributed to the existence of subradiant pair modes, depending on how strongly bounded they are to survive multiple scattering (for more details see Ref. [Skipetrov 2019] and its Supplemental Material). The fact that we observe large intensity fluctuations inside of the localized regime with the vectorial model with applied magnetic field, whereas with the scalar model only at the transition, could possibly be explained by the different properties of subradiant pair modes in those two models (as it was mentioned, the subradiant pairs in the vectorial model with applied magnetic field depend on the angle between the magnetic field axis and the pair axis and they are covering a large part of the complex plane of eigenvalue distribution: see for details, for example, Fig.

S5 of Supplemental Material of Ref. [Skipetrov 2019]).

Although there are some open questions yet to be answered, our results indicate that the variance of the intensity is a suitable indicator of the phase transition to the localized regime.

Experimentally, this could be achieved, for example, through configuration averages with along the frequency axis becomes less broad once the magnetic field is applied and it is also slightly less broad for m e = ±1 than for m e = 0, which could reflect the inhomogeneous broadening due to near-field interaction, which is then more suppressed for m e = ±1 once the magnetic-field is applied. Moreover, in addition to the localized modes that appear in the bottom of the eigenvalue distribution on Fig. III.11(c,e) for m e = ±1 when the strong magnetic field is applied, there are also long-lived collective modes with low values of IPR (yellow) emerging that are not present in the case when atoms are not subjected to the external magnetic field [Fig. III.11(a)]. This implies that the longest-lived subradiant modes could be resurrected by applying the strong magnetic field, which is expected to happen since the near-field interaction is at the origin of their absence. Note that we have not checked this by observing the late-time decay dynamics with applied magnetic field as it was out of the scope of our study. However, it would be most definitely interesting to study this, as switching the magnetic field on/off would introduce a way of switching on/off the longest-lived subradiant states, which would be of interest for potential applications of subradiant states. It would be as well interesting to study more closely the potential signature of Anderson localization in late-time decay-dynamics [Skipetrov 2016c], where a scaling analysis involving different sample sizes, atom number, laser detuning, etc. would have to be done to discriminate between the decay of subradiant and localized modes.

Conclusion

In this thesis, we have studied different properties of the long-lived collective modes which arise due to the interaction of light with ensemble of cold atoms, with the general aim of expanding our understanding of those modes in different regimes, mainly beyond the linear-optics and dilute atomic sample regimes.

In Chapter I we have shown experimentally that increasing the saturation parameter s(∆) of the driving field significantly enhances the population of the subradiant states: we demonstrated that the subradiant population increases by a factor of ∼ 200 in the experimentally explored range of s(∆). More importantly, we demonstrated the superlinear increase of the subradiant population as P sub ∝ s β with the power-law exponent β = 1.5, while the final subradiant lifetimes that were observed remained unaffected by the strength of the driving field. Since in the linear-optics regime the subradiant population scales linearly with s (P sub ∝ s), the observation of superlinear scaling demonstrates that, with the atomic samples of the experiment, we are in the regime beyond linear-optics for the subradiant modes, even with the saturation parameter as low as s(∆) ∼ 10 -3 . This superlinear scaling of subradiant population has been interpreted as a process similar to optical pumping, such that besides being directly populated by the driving laser, the subradiant states additionally gain in population through the decay processes from multi-excitation superradiant states.

Moreover, the comparison of experimental data with numerically obtained results indicated that this behavior is related to the generation of quantum correlations in the system. This suggests that subradiance might be a resource of quantum correlation or entanglement, and that increasing the driving field intensity is a way to overcome the weak coupling of those states in order to exploit them more efficiently, as for example in quantum memory devices for information storage.

The study presented in the Chapter II adresses (numerically) the question of how the near-field dipole-dipole interaction affects subradiance. For that, we have used the vectorial coupled-dipole model and we observed that the near field terms are detrimental for subradiance, which was seen in the reduction of the subradiant lifetimes with increasing density. In particular, we observed that the near-field effects have a negligible role on subradiant lifetimes in the low-density regime, however for higher densities we observed that subradiant lifetimes decrease with density for a given on-resonance optical thickness. These results have been benchmarked with the results of scalar coupled-dipole model (which does not consider the near-field interaction) for which the decrease of subradiant lifetimes with density has not been observed, indicating that in the vectorial model the decrease indeed comes from the near-field interaction. The detrimental character of near-field terms can be interpreted as being due to an effective inhomogeneous broadening induced by the near-field interaction, which effectively reduces the optical thickness as the density increases, which yields the reduction of subradiant lifetimes.

The detrimental character of near-field interaction had also been previously numerically The studies on Anderson localization of light obtained with the vectorial coupled-dipole model consider atoms with non-degenerate ground state, i.e. a J g = 0 → J e = 1 transition.

Our experiment, on the other hand, uses rubidium atoms, which have a more complex atomic structure (with hyperfine structure and Zeeman degeneracy in the ground states).

Although the coupled-dipole model could be extended to atoms with degenerate ground state [Pellegrino 2014, Jenkins 2016, Jennewein 2016, Facchinetti 2016, Sutherland 2017], the calculations are much more complex and they have not been performed to address the problem of light localization, so we do not now whether the results obtained for atoms with non-degenerate ground state can be simply extended to the case of atoms with degenerate ground state. Therefore, our group is currently building a new cold-atom experiment based on ytterbium atoms, which have a characteristic J g = 0 → J e = 1 transition, with the ultimate goal of observing Anderson localization of light in 3D cold-atom samples.

Some experimental challenges are still ahead. One is the impact of the residual thermal motion of cold atoms. Although we have shown that subradiant modes in dilute samples are robust against the atomic motion, another work addressing denser samples suggests that it may be detrimental for AL modes [Kuraptsev 2020]. Moreover, thermal motion will also be Another open question is the effect of the nonlinearity, which is hard to avoid if light is localized in a very small volume within the sample, as a few photons might be enough to produce a nonnegligible saturation parameter. In particular, our study on subradiance suggests that collective modes with long lifetimes have a lower saturation threshold than shorter-lived ones. We do not know if this conclusion still holds for localized modes. This will be the subject of further investigations.

Appendices

1 2 3

 123 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter I -Subradiance beyond the linear optics regime . . . . . . . I.Experimental setup and data acquisition. . . . . . . . . . . . . . . . . . . . . I.1.1 Cold atomic samples and probe beam pulses . . . . . . . . . . . . . . I.1.2 Experimental cycles and measurement protocol . . . . . . . . . . . . I.1.3 Calibration of the saturation parameter of the probe laser . . . . . . I.1.4 Spurious effects due to the high saturation parameter: Heating and pushing of the atomic cloud . . . . . . . . . . . . . . . . . . . . . . . I.Experimental subradiant data . . . . . . . . . . . . . . . . . . . . . . . . . . . I.2.1 Temporal dynamics of the scattered intensity: long lived decay and saturation parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . I.2.2 Subradiant lifetime and relative amplitude . . . . . . . . . . . . . . . I.2.3 Population of subradiant states . . . . . . . . . . . . . . . . . . . . . I.Interpretation and comparison with theory . . . . . . . . . . . . . . . . . . . I.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.3.2 Experimental vs. numerical results . . . . . . . . . . . . . . . . . . . I.3.3 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii xiv TABLE OF CONTENTS Chapter II -Subradiance beyond the dilute regime . . . . . . . . . . . 47

Figure 1 :

 1 Figure 1: Collective Dicke states of sample of n two-level atoms. Here I 0 represents the emission rate of an single atom. Adapted from Ref. [Dicke 1954].

Figure 2 :

 2 Figure 2: Time evolution of the normalized excited state population (black solid curve) after switching off the laser (linear-optics regime). At first, the population decreases faster (superradiance) than the single atom decay (black dashed line) and then slower (subradiance). The inset shows the emission diagrams at short times (blue) and long times (red) after the laser switch-off. Taken from Ref. [Bienaimé 2012].

Figure 3 :

 3 Figure 3: (a) Experimentally obtained fluorescence decay from cold-atom sample of different on-resonance optical depths b 0 , with a detuning δ = -3.15 (in units of natural linewidth) of the driving laser. All decay curves are normalized to the fluorescence level right at the switch off of the laser. For comparison, the theoretical single atom decay exp(-t/τ at ) is also shown (dash-dotted line). The longest subradiant lifetime τ sub is obtained from an exponential fit (dashed lines) at late time. (b) Measured subradiance lifetimes τ sub /τ at as a function of b 0 . All measured points collapse on a single line, independently of the detuning. The linear scaling of τ sub with b 0 is stressed by the linear fit (solid line). Taken from Ref. [Weiss 2018].

Finally

  , we have numerically studied the Anderson localization of light in 3D atomic samples with the focus on finding an unambiguous signature of light wave localization. By using scalar and vectorial coupled-dipole models we identified this signature to be the enhanced fluctuations of the scattered light. The part of the study with the scalar CD model was done by Florent Cottier [Cottier 2019a], while I did the study involving the vectorial CD model. This work led to the article [Cottier 2019b] and will be presented in Chapter III after a detailed introduction on the phenomena of Anderson localization, in particular of light waves.

Figure I. 1 :

 1 Figure I.1: Scheme of 87 Rb D 2 transition hyperfine structure. The transitions of interested for generating and probing cold cloud of 87 Rb atoms.

Figure I. 2 :Figure I. 3 : 4 :

 234 photons, i.e. the fluorescence signal, are recorded with a multichannel scaler (MCS6A by FAST ComTec) with a time bin of 1.6 ns while averaging over the cycles. Note that averaging over the experimental cycles corresponds to averaging over various spatial configurations of atomic positions. After applying the series of 12 probe pulses and simultaneous fluorescence measurement, the MOT is switched on again (beginning of a new cycle) so that most of the atoms are recaptured. This allows us to have a fast experimental cycle (the total duration of ∼ 110 ms), so that the measured fluorescence signal can be averaged over a large number of experimental cycles, typically ∼ 5 × 10 5 , which corresponds to more than 12 hours of data acquisition in total for only one set of parameters, such as probe saturation parameter or probe detuning. Note that for a couple of lowest probe saturation parameters, averaging over more cycles was required (up to several days).

  , the protocol for measuring the parameters of the cloud (size, number of cold atoms, optical thickness, density) simultaneously to the data acquisition is the following. The absorption imaging is performed instead of the fluorescence measurement once every 250 cycles. The goal is to characterize the cloud at times corresponding to each of the applied 12 pulses, but with taking into account the effect pulses have on the cloud. Hence, as illustrated on Fig. I.3, to measure the optical thickness and other parameters of the cloud corresponding to the nth pulse, we apply (n -1) probe pulses after which follows the absorption imaging at the time corresponding to the nth pulse. This protocol enables us to have a good calibration of b 0 and other parameters while probing our sample with laser pulses of different intensities.

  We used several approaches to properly calibrate the saturation parameter. The first calibration method is based on the measurement of the fluorescence level. The atomic cloud (b 0 = 37 ± 1) was illuminated by the probe beam with ∆ = -4Γ and the fluorescence signal was recorded by a photodetector. The fluorescence measurement was done as a function of s (∆) of the probe beam [Fig. I.5(a)]. Since the total scattering rate is ∝ s/(s + 1), we fitted the measured fluorescence levels by f = Bs /(Cs + 1) and we obtained the correction factor for the saturation parameter C = 0.36 ± 0.04.

I. 2 .

 2 3b)where p 21 = 1/2 is the decay probability of F = 2 → F = 1 transition, s 22 is the saturation parameter for the F = 2 → F = 2 transition with g 22 = 1/6 the degeneracy factor of the Zeeman sublevels of the transition F = 2 → F = 2, ∆ 22 = 44Γ + ∆ the detuning of the probe beam from that transition and s 0 = I/I sat the on-resonance saturation parameter (the same as in Eq.(I.2)). Therefore, by measuring the hyperfine depumping rate Γ depump , we can obtain the on-resonance saturation parameter s 0 from Eqs. I.3a and I.3b, and then the detuning-dependent saturation parameter s(∆) for the transition F = 2 → F = 3 from Eq.Because of the depumping losses, the number of atoms that can undergo a transitionF = 2 → F = 3 is decreasing as N (t) = N (t = 0)exp(-Γ depump t).As the fluorescence level is proportional to the number of atoms (at sufficiently large detuning), depumping rate can be obtained by measuring the steady-state fluorescence level as a function of the laser duration. In fact, since we illuminate the atomic cloud by a series of 12 probe pulses each of duration 10 µs, we added up the durations of subsequent pulses. By fitting an exponential to the measured fluorescence levels [Fig.I.5(b)], we obtained the depumping rate, from which we determined the saturation parameter s(∆) = (1.95 ± 0.09) • 10 -2 . Therefore, with s (∆) = (6.07 ± 0.01) • 10 -2 this method yielded the correction factor C = s/s = 0.32 ± 0.01. Another way of obtaining depumping rate is by measuring the optical thickness, which is proportional to the number of atoms (b 0 ∝ N/kR 2 ), as a function of the laser duration. We measured the optical thickness without and with repumping laser, which re-excites atoms back from the dark F = 1 state. Then, by fitting an exponential to the ratio of measured b 0 without and with repumper laser [Fig.I.5(c)], we extracted the depumping rate, from which we finally obtained s(∆) = (1.89 ± 0.07) • 10 -2 . Since s (∆) = (5.56 ± 0.01) • 10 -2 , with this calibration method we obtained the correction factor C = 0.34 ± 0.01. The correction factors of the saturation parameter obtained using the three described methods are in good agreement and with similar relative uncertainties. Therefore, we use the average correction factor C = 0.34 ± 0.02. It should also be noted that the precise calibration of the probe detuning ∆ has been done by measuring the optical thickness as a function of the detuning b

values of b 0

 0 , corresponding to different times of flight and thus different numbers of applied pulses, yield similar results (see Fig.I.14).

Figure I. 5 :

 5 Figure I.5:Different calibration methods of the saturation parameter of the probe beam. Black dots represent experimental data, while green solid lines are fit curves. a) Measured fluorescence level as a function of the saturation parameter s (∆) evaluated from the measurement of the peak intensity of the probe beam, together with the fit function f = Bs /(Cs + 1). Experimental data was obtained with a cloud of resonance optical thickness b 0 = 37 ± 1, a probe beam detuning ∆ = -4Γ and a pulse duration of 5 µs. b) Measured fluorescence level as a function of the total interaction time with the probe beam, fitted by an exponential decay, from which the depumping rate Γ depump = (1.51 ± 0.07) × 10 -4 Γ was extracted. The measurement was done with ∆ = -10Γ and a pulse duration of 10 µs. c) Ratio between measured b 0 without and with repumping as a function of the total interaction time with the probe beam. The exponential fit to the experimental data gave Γ depump = (8.37 ± 0.03) × 10 -5 Γ. The data were obtained with ∆ = -8Γ and a pulse duration of 10 µs.

  Figs. I.6 and I.7 we show examples of measured fluorescence decay for various values of the probe saturation parameter. It should be noted that we directly compare decay curves of different s(∆) corresponding to a given probe detuning ∆ and on-resonance optical thickness b 0 since it has been shown that the lifetime of collective modes in large extended cold atomic clouds scales with the cooperativity parameter b 0 ∝ N/(kR) 2 [Roof 2016, Guerin 2016, Cipris 2021b]

Figure I. 6 :

 6 Figure I.6: Temporal dynamics of light emitted after the laser switch off from the cloud of optical thickness b 0 ≈ 50 (a-b), b 0 ≈ 80 (c-d) and b 0 ≈ 160 (e-f) driven by the probe laser with the detuning ∆ = -4Γ for different values of the saturation parameter s(∆) (color code). Panels on the left (a,c,e) show the fluorescence decay normalized with the steady-state fluorescence level, while the curves on the right panels (b,d,f) correspond to the same, but non-normalized decay curves. Note that the y-axis of the plots showing non-normalized decay curves (labeled I) corresponds to the photon count rate, i.e. total photon count per second (and per experimental cycle). The black dashed line on left panels stands for the single-atom decay I(t)/I 0 = exp(-Γt).

  Figure I.7: The same as in Fig. I.6 but for ∆ = -2Γ (a-d) and ∆ = -8Γ (e-f). The optical thickness is b 0 ≈ 50 (a,b,e,f) and b 0 ≈ 80 (c-d). Panels on the left show decay curves normalized with the steady-state intensity, while panels on the right non-normalized decay curves for different values of the saturation parameter s(∆).

Figure I. 8 :Figure I. 9 :Figure I. 10 :

 8910 Figure I.8: Subradiant lifetimes τ sub (a,c,e) and corresponding relative subradiant amplitudes Ãsub (b,d,f) obtained from the fit I(t)/I 0 = Ãsub exp(-t/τ sub ) of the normalized decay curves [Fig. I.7(a,c)] over different time windows Γt Fit . The on-resonance optical thickness of the cold-atomic cloud is b 0 ≈ 50 (a,b), b 0 ≈ 80 (c,d) and b 0 ≈ 160 (e,f), while the detuning of the probe beam is ∆ = -4Γ.

  Figure I.11: (a) Temporal dynamics of the scattered intensity after the switch-off of the driving field (at t = 0), for different saturation parameters s(∆) as well as for ∆ = -2Γ (solid line) and ∆ = -4Γ (dashed line). The decay curves for both detunings collapse, which indicates that the temporal dynamics is independent of the detuning. The intensity curves are normalized by s. Even after this normalization, the amplitude of the slow decay (corresponding to the long-lived collective modes) increases with s, which illustrates well the super-linear growth of the long-lived mode populations. Inset: Steady-state scattering rate Γ st as a function of the saturation parameter s for b 0 ≈ 50 and ∆ = -2Γ. The measured Γ st is represented by full circles, while the solid line corresponds to the fit function Γ st ∝ s/(s + 1). (b) Subradiant amplitude A sub as a function of s(∆) for detunings ∆ = -2Γ (• symbol), ∆ = -4Γ ( * symbol) and ∆ = -8Γ ( symbol) as well as for different fit intervals t Fit represented by different colors as indicated on the legend. The inset shows the corresponding lifetimes τ sub . A sub and τ sub were obtained by fitting I(t) = A sub exp(-t/τ sub ) to the non-normalized decay curves (which are for example shown on Fig. I.7(a)) over a given time window t Fit . The on-resonance optical thickness is b 0 ≈ 50 for the results shown on both (a-b) panels.

Figure I. 14 :

 14 Figure I.12: Normalized population P sub = A sub,tot τ sub /N of long-lived collective modes as a function of the saturation parameter s(∆) for the cold atomic cloud with optical thickness b 0 = 54 ± 2, as well as for laser drive with different detunings ∆. Amplitudes and lifetimes of subradiant modes are obtained from the fluorescence decay curves over the time windows t Fit ∈ [150; 250]/Γ (a) and t Fit ∈ [20; 40]/Γ (b). The black solid lines correspond to a fit P ∝ s β , with β = 1.49 ± 0.16 (a) and β = 1.53 ± 0.14 (b). Insets: Corresponding subradiant lifetime τ sub as a function of s(∆).

  On Fig. I.15(a), we show the direct comparison of experimentally and numerically (from QPC simulations) obtained normalized population of subradiant modes [Eq. (I.4)]. Although there is some discrepancy in the values of the normalized population (which might be due toP

Figure I. 15 :

 15 Figure I.15: (a)Normalized population P sub of long-lived states as a function of the saturation parameter. Experimental data (symbols, left axis) acquired for b 0 = 54±2, with error bars describing the 95% confidence bounds (statistical uncertainty only). QPC simulations (lines, right axis) realized for b 0 ≈ 5 and N = 100, averaging over 40 realizations (error bars of order 1%, not shown here). The red dash-dotted line stands for a fit P sub ∝ s β , with β = 1.007 ± 0.010 (fit on numerical data), indicating the linear-optics (LO) regime for subradiant states. The blue dotted line corresponds to a power-law fit P sub ∝ s β with β = 1.49 ± 0.23 (fit on experimental data), which indicates the regime where the population of subradiant states is being enhanced by a process similar to optical pumping (OP) via superradiant states. Insets: Corresponding lifetimes of long-lived states as a function of the saturation parameter obtained from the experiment (top panel) and from the simulations (bottom panel). (b) The decay dynamics of the population P n of the n-excitation states for a cloud of N = 9 atoms obtained with exact simulations. The atomic cloud is initially driven to steady-state by a strong resonant electric field with Ω 0 = 5Γ. We observe that the lower the excitation number n, the slower the decay. The numerical results presented here were obtained by N. A. Moreira and T.S. do Espirito Santo.

  photon, i.e. single-excitation eigenstates |± = (|eg ± |ge )/ √ 2 (only one excited atom), in addition to the ground |gg and double-excited |ee states (both atoms are excited) [see Fig. I.16(b)]. The symmetric collective state, i.e. superradiant state, |+ presents decay rate Γ + = Γ(1 + sin(kr 12 )/(kr 12 )) (Γ + > Γ for r 12 λ) and energy shift ∆ + = Γ cos(kr 12 )/2kr 12 , while the decay rate and energy shift associated to the antisymmetric (subradiant) state are Γ -= Γ(1 -sin(kr 12 )/(kr 12 )) (Γ -< Γ for r 12 λ) and ∆ -= -Γ cos(kr 12 )/2kr 12 , respectively. The pump couples efficienty to the superradiant |+ state, while very weakly to the long-lived, i.e. subradiant |-one. Moreover, we introduce the effective Rabi frequency for each single-excitation mode |+ and |-: Ω + = √ 2 cos(kr 12 /2)Ω 0 and Ω -= √ 2 sin(kr 12 /2)Ω 0 (up to a phase), which can be identified by rewriting the driving Hamiltonian in terms of the states |+ and |-.

  Figure I.16: Dicke space for N 1 two-level atoms, where the downward arrows depict the decay processes. The processes in green refer to the Dicke superradiant cascade through symmetric states [Dicke 1954], and the red ones to the linear-optics (LO) processes (i.e., the single-excitation regime). (b) Energy levels for N = 2 coupled atoms. (c) Computed population of the long-lived |-state for N = 2 coupled atoms. The red dash-dotted curve corresponds to P -= s -/2 ∝ s, while the dashed blue one scales as s 2 (labelled OP). Simulation realized for N = 2 atoms distant of kr = 0.05, aligned with the pump axis, with a detuning ∆ = -500Γ. Simulation has been done by R. Bachelard.

  sample of N 2 atoms [seeFig. I.16(a)]: although weakly coupled to the external drive and consequently being weakly populated directly by the drive, the long-lived modes can be efficiently populated via the decay processes from the multi-excitation superradiant states, which are well coupled to the external drive. In other words, the experimentally and numerically obtained super-linear growth of subradiant population P sub ∝ s β , with β ∼ 1.5, can be attributed to a mechanism analogous to optical pumping via superradiant states, such that the long-lived modes, in addition to being directly populated by the probe laser, are gaining population through the decay from multi-excitation superradiant states. Let us note that it is not clear why we observe the power-law scaling P sub ∝ s β , with β ∼ 1.5, in contrast to P sub ∝ s 2 for the simple case of N = 2 atoms. The difference in the power-law exponent might be due to a vastly larger Hilbert space for N 2 atoms, and thus much larger number of collective states with different decay rates, i.e. lifetimes.Using numerical simulations, the study of the dynamics of the many-excitation states reveals that the longest lifetimes are found in the subradiant modes with fewer excitations, toward which the system quickly decays, i.e. the single-excitation long-lived modes have the longest lifetimes, while the multi-excitation long-lived modes have shorter lifetimes [see Fig. I.15(b)].

[

  Pellegrino 2014], and it prevents Anderson localization[Bellando 2014, Skipetrov 2014, Máximo 2015, Máximo 2019], as well as large refractive indices[Andreoli 2021].

  ζ |β ζ j (t)| 2 gives the probability of atom j being excited at time t, while |β ζ j (t)| 2 represents the probability that atom j is excited to state |J e = 1, m e = ζ at time t. The first term of the CDEs [Eq. (II.1)] describes the single-atom dynamics. The second one represents the driving by the external electric field. The last term is the coupling term that corresponds to the dipole-dipole interaction, and it gives rise to collective effects like superradiance, subradiance, and multiple scattering. The dipole-dipole interaction mediated by the electric field produced by radiating dipoles is characterized by the dyadic vectorial Green's function given by

2 )

 2 with rζ = êζ • r the component of the unit vector r = r/r along the direction ζ = 0, ±1. It is the imaginary part of the Green's function that is responsible for collective energy shifts, while its real part gives rise to collective decay rates, of which some are superradiant Γ sup > Γ 0 and other subradiant Γ sub < Γ 0 . Note that the vectorial Green's function [Eq. (II.2)] contains both far-field (1/r) and near-field contributions (1/r 2 and 1/r 3 ), as well as polarization coupling (∝ rζ r * η for ζ = η) through both near-and far-field terms. Considering G ζ,η as a spatial component of a symmetric tensor ↔ G, Eq. (II.2) can be written in more compact form: ↔

  .15) Diagonal terms of matrix G describe the single-atom dynamics, whose real part corresponds to the single-atom decay rate Γ 0 = -2 R(G ζ,ζ j,j ), while the imaginary part corresponds to the atomic frequency in the rotating frame of the external driving field, i.e., it corresponds to the detuning ∆ 0 = ω L -ω 0 = I(G ζ,ζ j,j ). On the other hand, the off-diagonal terms are

Figure II. 1 :

 1 Figure II.1: (a) Eigenvalue distribution of vectorial coupled-dipole modes for an atomic sample with b 0 = 30, ρλ 3 = 15 (and N = 1102). The number of realizations is N r = 25. Each point represents one of the eigenmodes with corresponding frequency shift ∆ n = ω L -ω n [Eq. (II.16)] and decay rate Γ n [Eq. (II.17)]. (b) Eigenvalue spectrum corresponding to the scalar coupled-dipole model [Eq. (II.19)] for an atomic sample with b (s) 0 = 30, ρλ 3 = 15 (and N = 3721). The number of realizations is N r = 22. The color code in both panels represents the inverse-participation rate (IPR) [Eqs. (II.18) and (II.20) for vectorial and scalar model, respectively]. Data points with low IPR (yellow) correspond to the collective modes in which N n 1 atoms are significantly participating (IPR 1). Red data points are modes of very close pairs of atoms created in the sample, i.e, modes to which 2 atoms are significantly contributing (IPR ≈ 1/2 for pair modes in scalar model; IPR ≈ 1/6 -1/2 for pair modes in vectorial model due to 3 Zeeman substates of excited state). Black spirals, i.e., pair branches correspond to the analytical expressions of pair eigenvalues [given by Eqs. (II.21) and (II.22) for vectorial model and Eq. (II.23) for scalar model]. The exclusion volume for atoms has not been used here. The detuning of the external driving field is ∆ 0 = -15Γ 0 (equivalently in the scalar model ∆ 0 = -15Γ (s) 0 ).

Figure II. 2 :

 2 Figure II.2: Complex eigenvalues λ pair = i∆ pair -Γ pair /2 corresponding to the modes of atom pair (N = 2) in the vectorial coupled-dipole model [Eqs. (II.21)-(II.22)]. Decay rates Γ pair (a) and frequency shifts (with respect to the atomic resonant frequency ω 0 ) ∆ pair = ω 0 -ω pair (b) of pair modes as a function of the distance between two atoms k 0 r 12 . (c) Pair branches, i.e., Γ pair (k 0 r 12 ) vs. ∆ pair (k 0 r 12 ) representation of (a) and (b) combined.

  Fig. II.2 we show the real Γ pair = -2 R(λ pair ) and the imaginary ∆ pair = I(λ pair ) parts of pair eigenvalues [Eqs. (II.21) and (II.22)] as a function of the distance between two atoms.

  . In the limiting case k 0 r 12 1, the superradiant decay rate is Γ ,⊥ pair,+ ≈ 2Γ 0 , while the subradiant lifetime is τ ,⊥ pair,-→ ∞. This corresponds to a well-known Dicke limit[Dicke 1954]. Moreover, the frequency shift is |∆ ,⊥ pair,± | Γ 0 for k 0 r 12 1. In Fig. II.2(b-c), we show only a short range of ∆ pair since the frequency shift is extremely large for k 0 r 12 1. On the other hand, when the separation between two atoms is k 0 r 12 1, super-and subradiant pairs cease to exist, i.e., decay rates and frequency shifts become equal to those of a single atom: Γ ,⊥ pair,± ≈ Γ 0 and |∆ ,⊥ pair,± | ≈ 0, respectively. Note that this is only the case for N = 2 atoms; for macroscopic samples with N 1, super-and subradiant collective effects still prevail for k 0 r ij 1, due 2 In the scalar model, the pair eigenvalues λ are obtained by diagonalizing the 2 × 2 scalar Green's matrix [Eq. (II.19)] yielding λ

  pair,± ≈ ∆ 0 . On the other hand, if two atoms are very close, kr12 1, pair decay rates and frequency shifts are Γ (s) pair,+ ≈ 2Γ

Figure II. 3 :

 3 Figure II.3: (a) Decay curves, i.e. normalized intensity of the scattered light I(t)/I 0 (all polarization channels), with I 0 = I(t = 0), for several b 0 and ρλ 3 obtained from the vectorial CDEs [Eq. (II.6)] after the switch-off of the driving field (t > 0). (b) Instantenous decay rate Γ(t) of the corresponding decay curves obtained from Eq. II.28. Black solid lines are fit curves [Eq. (II.27)] in the fit range I(t)/I 0 = [10 -6 , 5 × 10 -6 ] [panel (a)] and the corresponding decay rates [panel (b)]. The black dotted line presents the decay curve of the single atom I(t)/I 0 = exp(-Γ 0 t) in panel (a) while in (b) it shows the corresponding single-atom decay rate Γ 0 . The incident electric field has ˆ L = ê-1 circular polarization, while its detuning is ∆ 0 = -15Γ 0 . The scattered light is collected at θ sc = 45 • . We used the exclusion volume r min = ρ -1/3 /π.

  Fig. II.3, where panel (a) shows decay curves for several b 0 and ρλ 3 , and panel (b) shows the instantaneous decay rate Γ(t) which corresponds to the decay curves in panel (a). The instantaneous decay rate has been obtained from the following equation Γ(t) = -d[ln I(t)] dt (II.28) by assuming that at each moment the decay curve is an exponential function with the instantaneous decay rate: I(t) ∝ exp(-Γ(t)t). As it can be seen in Fig. II.3(b), light emission is superradiant (Γ(t) > Γ 0 ) at early times after the switch-off of the driving field. As for long-lived light emission (Γ(t) < Γ 0 ), there are intermediate times for which the decay rate

  decay within that interval (see examples of the fit curves in Fig. II.3(a)-(b)

  Fig. II.3(b), the instant of time when late-time dynamics starts is different for different parameters of the sample: for larger b 0 late-time dynamics appears later in time (for example, t ∼ 20τ 0 for b 0 = 6 and t ∼ 80τ 0 for b 0 = 64). Therefore, to have the same time-fit interval at late times for all decay curves of different b 0 = [2; 72] that we are studying

Figure II. 4 :

 4 Figure II.4: Subradiant lifetime τ sub as a function of the on-resonance optical depth b 0 for several densities of the sample ρλ 3 , computed from the vectorial coupled-dipole model. The scattered light is collected at θ sc = 45 • . The lifetimes τ sub are obtained from an exponential fit of the total scattered light intensity (all polarizations) in the fit range I/I 0 = [10 -6 , 5 × 10 -6 ]. For ρλ 3 > 5, the linear scaling of τ sub with b 0 (independently of ρ) fails, as data corresponding to the same b 0 but different ρλ 3 do not collapse onto a single curve. Other system parameters: ˆ L = ê-1 circular polarization of the driving field, ∆ 0 = -15Γ 0 , r min = ρ -1/3 /π.

  Fig. II.5(b). For b 0 5, the subradiant lifetimes are approximately constant in the entire range k 0 R = [8; 24]. On the other hand, for b 0 20, the subradiant lifetimes are decreasing with decreasing k 0 R for all k 0 R = [7; 15]. There is no threshold k 0 R for which the reduction

Figure II. 5 :

 5 Figure II.5: Subradiant lifetime τ sub for several b 0 as a function of (a) density ρλ 3 and (b) corresponding size of the sample k 0 R = (2π) 3/2 b 0 3ρλ 3 . The insets in (a) are close-ups of the three lowest b 0 data sets. The lifetimes τ sub were extracted from an exponential fit of the scattered light intensity after the laser switch-off in the fit range I/I 0 = [1 • 10 -6 , 5 • 10 -6 ]. c) Subradiant part of the scattered intensity I(t)/I 0 as a function of time for b 0 = 4.9 (blue set of curves), b 0 = 13.9 (green curves), b 0 = 39.4 (purple curves) and for different ρλ 3 which are represented by different color shades and different types of lines. Those decay curves correspond to some of data shown in panels (a)-(b).The yellow curve corresponds to the decay of a single atom. For ρλ3 5 the decay at late times becomes slightly faster as the density is increased. The inset shows the complete decay curves from the moment the laser is switched-off (t = 0). The intensity of the scattered light was computed from the vectorial CDEs [Eq. (II.6)]. Other simulation parameters: ˆ L = ê-1 circular polarization of the driving field, ∆ = -15Γ 0 , θ sc = 45 • , r min = ρ -1/3 /π.

  Fig. II.6 we show the eigenvalue distribution of collective modes [Eqs. (II.17) and (II.16)] for several b 0 and ρλ 3 . Since we did not impose any exclusion volume on atoms in this case, very strong

Figure II. 6 :

 6 Figure II.6: Eigenvalue distribution, Γ n vs. ∆ n = ω L -ω n [Eqs. (II.17) and (II.16)], of the coupled-dipole modes computed without any exclusion volume and for b 0 = 16, ρλ 3 = 5 (a,b), b 0 = 30, ρλ 3 = 20 (c,d) and b 0 = 68, ρλ 3 = 40 (e,f).For the left panels, the drive is red-detuned ∆ 0 = -15Γ 0 , while for the right panels, it is blue-detuned ∆ 0 = +15Γ 0 . In addition to collective subradiant modes (IPR 1), there are two subradiant branches (and two superradiant branches) that consist of modes of close pairs of atoms (IPR ∼ 1/6 -1/2). Branches labeled as and ⊥ represent pair eigenmodes corresponding to two atomic dipoles that are oriented along the direction of their alignment and orthogonally to their alignment axis, respectively. Two subradiant pair branches are asymmetric: they have different frequencies and lifetimes for a given inter-atomic distance. Without any exclusion volume, the extent of branches is very long, and the driving field is resonant (∆ n = 0) with some of the pair modes (it depends on the sign of the detuning which branch crosses the resonance ∆ n = 0). This results in highly populated modes with ∆ n = 0, and the entire late-time decay dynamics is strongly governed by their decay [see Fig.II.8(a)]. Since we are interested in collective subradiance, we use the exclusion volume to restrict the pair branches' extent to prevent strong coupling of the drive with pair modes. The extent of the branches with the exclusion volume r min = ρ -1/3 /π is illustrated with the black solid lines that correspond to the analytical expressions given by Eqs. (II.21) ( ) and (II.22) (⊥) with the cut-off (the maximum extent) corresponding to the exclusion volume r min . As r min is smaller for larger ρ, the extent of the branches with the imposed density-dependent exclusion volume is getting larger with increasing density. Red solid lines (left panels) and blue solid lines (right panels) highlight the resonance ∆ n = 0. The color code corresponds to the inverse-participation ratio (IPR) given by Eq. (II.18).

  Fig. II.6 with the density-dependent exclusion volume, the collective modes have the longest subradiant lifetimes, not the pair modes. We have checked that this observation holds for all system parameters corresponding to data shown in Figs. II.4-II.5. Moreover, those subradiant lifetimes, in Figs. II.4-II.5, obtained from the exponential fit at late times, are longer than the lifetimes of the remaining pairs with the r min cut-off [see Fig. II.7(a) for the maximum subradiant pair lifetimes for given exclusion volume and density]. This indeed suggests that the obtained subradiant lifetimes correspond to the collective subradiant modes.

Figure II. 7 :

 7 Figure II.7: Maximum subradiant pair lifetimes τ pair,-= 1/Γ pair,-(a) and maximum pair frequency shifts ω 0 -ω pair,-(b) of and ⊥ subradiant pair branches (green and purple lines, respectively) as a function of the exclusion volume k 0 r min which defines the minimum inter-particle distance.The subradiant lifetimes and frequency shifts were obtained from real and imaginary part of pair eigenvalues given by Eq. (II.21) ( ) and Eq. (II.22) (⊥) for the inter-particle distance defined by the exclusion volume. The subradiant pair lifetimes and (absolute) pair frequency shifts for a given k 0 r min are larger for branch than for ⊥ branch. Eigenvalues marked with different symbols correspond to the density-dependent exclusion volume r min = ρ -1/3 /π for several densities. Inset of panel (a) shows the exclusion volume r min = ρ -1/3 /π as a function of the density. The frequency shift of pairs in panel (b) is expressed relative to the atomic transition frequency ω 0 . To obtain the frequency shift relative to the frequency of the driving field ω L , the detuning of the drive ∆ 0 has to be taken into consideration: ∆ pair = ω L -ω pair = ω 0 -ω pair + ∆ 0 .

Figure II. 8 :

 8 Figure II.8: Influence of pairs on collective dynamics in the vectorial model for the atomic sample with b 0 = 16 and ρλ 3 = 10. (a) Decay of the scattered light (all polarizations) collected at θ sc = 45 •for opposite sign detunings and different exclusion volumes. Without any exclusion volume, the entire long-lived dynamics is influenced by pairs, with a noticeable difference between detunings of opposite signs (green lines). For r min = ρ -1/3 /π, the decay dynamics is independent of the detuning sign. With r min = 0.5/k 0 , only the early and intermediate dynamics is independent on the sign of the detuning, and for very late times, the decay for ∆ = +15Γ 0 starts to be slower than for ∆ = -15Γ 0 , showing the influence of pairs. (b) Decay of the light scattered at θ sc = 90 • decomposed into two polarization channels: parallel (I par ) and perpendicular (I perp ) to the scattering plane. Decay curves of both channels are normalized by the total scattered intensity (I = I par + I perp ) at t = 0. Without any exclusion volume, both superradiant and subradiant parts are partially polarized. With the exclusion volume r min = ρ -1/3 /π, the superradiant part is polarized while the subradiant part is fully depolarized. c) Two polarization channels of light scattered at θ sc = 90 • for the exclusion volume r min = 0.5/k 0 . For negative detuning, the entire subradiant part is depolarized, whereas, for positive detuning, light is depolarized at intermediate times and at late times; at very late times, the light becomes polarized due to the influence of pairs.

  Fig.II.8(a) where we concluded that, with r min = 0.5/k 0 and positive detuning, late-time dynamics is affected by subradiant pairs.

  the subradiant pair modes which are frequency shifted [Fig.II.9(a)]. In Fig.II.9(b) we show decay curves obtained with ∆ 0 = 0 and without any exclusion volume for two different ρλ3 of the atomic cloud while b 0 is the same. One can see that the decay is slower for ρλ 3 = 8 than for ρλ 3 = 16. Therefore, we again observe the same density effect (the reduction of subradiant lifetimes with increasing density) as in Sec. II.2.2 with off-resonant driving field and density-dependent exclusion volume.In Fig.II.9(b) we also present decay curves obtained with the drive on resonance but with the exclusion volume r min = 0.5/k 0 . The observations of the density effects are the same as those without any exclusion volume. Decay curves of the same ρλ 3 , but obtained with and without the exclusion volume, are similar. The slight discrepancies between those curves are not due to the pairs since there is no influence of pairs with the drive on resonance whether or not the exclusion volume is used (this was confirmed by checking the polarization properties of the scattered light; we observed fully depolarized subradiant decay dynamics with or without exclusion volume). The discrepancies are due to positional correlations that will be discussed in the following subsection.

  Figure II.9: (a) Eigenvalue distribution, Γ n vs. ∆ n = ω L -ω n [Eqs. (II.17) and (II.16)] of coupleddipole modes in vectorial model with the drive on atomic transition resonance ω L = ω 0 and without any exclusion volume. As the drive is strongly-coupled to the collective ∆ n ≈ 0 modes, the long-lived collective modes are significantly more populated comparing to subradiant pair modes which are frequency shifted. As a consequence, even without any exclusion volume, the late-time decay dynamics is not influenced by pairs. Only at extremely late-times, when the collective modes lose the population, the pair decay dynamics would emerge. Black solid lines correspond to the pair branches obtained from the analytical expression of pair eigenvalues [Eqs. (II.21) and (II.22)] with the cut-off given by the exclusion volume r min = 0.5/k 0 . The colorcode corresponds to the inverse-participation ratio given by Eq. (II.18); collective modes are characterized by IPR 1, while pair modes by IPR∼ 1/6 -1/2. Here b 0 = 20, ρλ 3 = 8 (b) Temporal decay dynamics of the scattered light (all polarization channels) obtained from the vectorial CDEs [Eq. (II.6)] after the switch-off of the driving field for b 0 = 20 and two different ρλ 3 , with and without the fixed exclusion volume r min = 0.5/k 0 . The drive is on resonance, ∆ 0 = 0, and the light is collected at θ sc = 45 • . The long-lived dynamics is slower for ρλ 3 = 8 than for ρλ 3 = 16. This confirms that the observed reduction of subradiant lifetimes at increasing densities in vectorial coupled-dipole model with the off-resonant drive is not caused by subradiant pairs or by positional correlations. Here the small discrepancy between the decay curves obtained with and without an exclusion volume can be attributed to the positional correlations that arise when the exclusion volume is imposed on the atoms.

  different populations and different lifetimes and frequency shifts. When using an off-resonant driving field ∆ 0 = 0 without imposing any exclusion volume on the atoms, some of the pair modes are resonant with the drive since pairs created in the sample have large frequency shifts (Fig.II.6). This strong coupling of the driving field with resonant pairs then results in a significantly larger population of pair modes than collective modes, and the decay dynamics is dominated by pairs (Fig.II.8) [Fofanov 2021]. Therefore, to remove the significant influence of pairs for a given ∆ 0 , one has to prevent the strong coupling of the driving field with any of the pair modes. The way to do that is to use such an exclusion volume that cuts the extent of the pair branches before they reach the resonance ∆ n = ω L -ω n = 0. With ∆ 0 = ±15 and with the exclusion volume r min = ρ -1/3 /π for the densities explored in this chapter, we indeed cut the extent of the branches before reaching the resonance ∆ n = 0 (Fig.II.

Figure II. 10 :

 10 Figure II.10: Comparison of subradiant lifetimes with different types of exclusion volume in vectorial coupled-dipole model.(a) Subradiant lifetimes τ sub as a function of the on-resonant optical thickness b 0 for several densities obtained with the density-dependent exclusion volume r min = ρ -1/3 /π (× symbols) and with the fixed exclusion volume r min = 0.5/k 0 (♦ symbols). The driving field is red detuned ∆ 0 = -15Γ 0 . (b) Comparison of subradiant lifetimes obtained with exclusion volumes r min = ρ -1/3 /π (+ symbols) and r min = 0.5/k 0 (• symbols) for the blue detuned driving field ∆ 0 = +15Γ 0 . The color code in panels (a-b) corresponds to the density of the sample (legend) and it is the same for both types of exclusion volume. (c) Subradiant lifetimes as a function of the minimum separation distance between atoms, i.e., exclusion volume r min , with b 0 = 16, ρλ 3 = 10 and a red-detuned driving-field ∆ 0 = -15Γ 0 . The observed increase of τ sub with increasing exclusion volume is due to the positional correlations. For all panels, subradiant lifetimes have been obtained from the exponential fit of the scattered light intensity (all polarization channels) collected at θ sc = 45 • in the fit range I/I(0) = [10 -6 , 5 × 10 -6 ].

  Fig. II.10(a) we observed that the lifetimes obtained with two different types of exclusion volume, the density-dependent exclusion volume and the smaller fixed exclusion volume, lead to very similar results. Apart from the lowest b 0 's, where we observe the influence of pairs for r min = 0.5/k 0 , the slight difference between the results obtained with the two types of exclusion volume can be attributed to a very small influence of the correlations introduced by the exclusion volume. The influence of correlations can be better characterized by plotting the subradiant lifetime as a function of r min [Fig.II.10(c)], where we observe that the subradiant lifetime slightly increases with increasing exclusion volume. Note that this behaviour does not come from pairs: pairs would produce the opposite effect since the increase of minimum inter-particle distance imposed by the exclusion volume reduces the subradiant pair lifetimes. Therefore, this implies that an increase of the positional correlations yields an increase of the subradiant lifetime. Since the positional correlations are also enhanced with increasing density for the fixed exclusion volume, we would expect to see the ehancement of subradiant lifetimes with increasing density if there were an influence of positional corelations. However, we observe the same effect of reduction of subradiant lifetimes with increasing density for the fixed exclusion volume as with the density-dependent exclusion volume [Fig.II.10(a)]. As a consequence, the subradiant lifetime reduction observed for increasing density cannot be attributed to positional correlations.Moreover, the test done with the drive on resonance and without any exclusion volume (hence, without any positional correlations induced by the exclusion volume), yields as well the decrease of τ sub with ρλ 3 for a given b 0 [Fig.II.9(b)], confirming that positional correlations are not at the origin of the reduction of τ sub with ρλ 3 .

Figure II. 11 :.Figure II. 12 : 3 .

 11123 Figure II.11: Comparison between the vectorial and scalar coupled-dipole models. (a) Temporal evolution of the intensity of the scattered light for several b 0 and ρλ 3 = 1.5 (set of red lines) and ρλ 3 = 15 (blue set of lines). Dotted lines represent the decay curves obtained from the vectorial CDEs [Eq. (II.6)], while solid lines correspond to the scalar CDEs [Eq. (II.10)]. The color code is the same for both vectorial and scalar curves and it corresponds to different b 0 . (b) Comparison of subradiant lifetimes of scalar and vectorial model plotted as a function of the on-resonant optical thickness for several densities of the sample. Lifetimes have been obtained from the exponential fit of decay curves in the fit range I/I(0) = [2 × 10 -7 , 10 -6 ]. For both panels and for both, vectorial and scalar, models: r min = ρ -1/3 /π, θ sc = 45 • , ∆ (v,s) 0 = -15Γ (v,s) 0 . In the vectorial model we use the vectorial definitions of the single-atom decay rate Γ (v) 0 and optical thickness b (v) 0 , while in the scalar model we use their scalar versions: Γ (s) 0 = (3/2)Γ (v) 0 and b (s) 0 = (2/3)b (v) 0 .

  the sake of clarity with scalar and vectorial definitions. To eliminate subradiant pairs and positional correlations as a source of increase of τ sub with increasing density for a given b 0 in the scalar model, we have done similar tests as in Sec. II.2.3 for the vectorial model. The results presented in Fig. II.11 are obtained with the density-variable exclusion volume r min = ρ -1/3 /π and the negative detuning ∆ 0 = -15Γ 0

Fig

  subradiant lifetimes obtained using the scalar CDEs as a function of ρλ3 [panel (a)] and as a function of the corresponding sample size k 0 R =

Figure II. 13 :

 13 Figure II.13: Comparison of the subradiant lifetimes obtained using different fit windows in (a)vectorial coupled-dipole model and (b)scalar coupled-dipole model. Different fit windows are represented by different symbols, while the color code corresponds to different densities and it is the same for both panels. While in the vectorial model we observe b 0 enhancement of the subradiant lifetimes and the density effects with the fit windows I(t)/I(0) < 10-4 , in the scalar model that is the case for lower intensity-levels I(t)/I(0) < 10 -6 , i.e., for even later times.

  Fig. II.3 we have shown the decay curves from the vectorial model for several b 0 and the corresponding temporal dependence of instantaneous decay rate. At late times, there is a slow temporal evolution of the lifetimes and more importantly, the well-known b 0 enhancement of the subradiant lifetimes. However, prior to this late-time decay dynamics, i.e. at intermediate times, there is a fast temporal evolution of subradiant lifetimes and interestingly, an overlap region where the system corresponding to different b 0 exhibits the decay with the same lifetimes at the same time. In fact, the larger b 0 is, the overlap with the decay curve corresponding to the subsequent b 0 lasts longer. This can be understood if we recall that for larger b 0 , the subradiant modes with the longer lifetimes are created (see, for example, Fig. II.6). Hence, for intermediate times, where we observe the overlap of the decay curves of different b 0 , modes common to the atomic clouds of different b 0 decay at the same time. Subsequently, at late times, when those modes lose the population, we observe the decay of the longest subradiant modes (which are longer for larger b 0 ). In Fig. II.13(a) we show subradiant lifetimes from the vectorial model obtained using different fit windows. We can see that for fit windows I(t)/I(0) > 10 -4 , the plateau of subradiant lifetimes emerges as b 0 increases, since for those fit windows, there is the overlap of decay curves of different b 0 corresponding to decay of modes that are common to systems with different b 0 .However, in the scalar model the overlap lasts longer than in the vectorial model for the same system parameters (at higher densities). In Fig.II.13(b) for the scalar model, we can see that the plateau in b 0 evolution of subradiant lifetimes appears for lower intensity-level fit windows I(t)/I(0) > 10 -6 , i.e., for even later times, than in the vectorial model. The reason for this is that in the scalar model at higher densities the lifetime extent of the collective modes is longer than in the vectorial model (for the atomic samples of the same b 0 and ρλ 3 in both models; see, for example, Fig. II.1).

Figure II. 14 :

 14 Figure II.14: Eigenvalue distribution for two different densities ρλ 3 = 6 (a,c) and ρλ 3 = 16 (b,d) obtained using the vectorial (a,b) and scalar (c,d) models. The resonant optical thickness is b (v,s) 0 = 20 for all panels and the number of realizations is N r = 40 (a), N r = 285 (b), N r = 36 (c) and N r = 255 (d) such that the number of eigenvalues is the same in all cases. Here there is no exclusion volume and the detuning of the drive is ∆ 0 = -15Γ 0 . The color code corresponds to the logarithm of the eigenstate density.

Figure II. 15 :

 15 Figure II.15: Eigenvalue distribution for two different densities ρλ 3 = 20 (a) and ρλ 3 = 42 (b), obtained using the vectorial coupled-dipole model. The resonant optical thickness is b 0 = 50 for both panels and the number of realizations is N r = 28 (a), N r = 130 (b), such that the number of eigenvalues is the same in all cases. Here there is no exclusion volume and the detuning of the drive is ∆ 0 = -15Γ 0 . The color code corresponds to the logarithm of the eigenstate density.

  Sec. III.1 we present the mechanism behind the localization of waves, as well as the discussion on the critical amount of disorder needed for the localization to occur (localization criteria). In Sec. III.2 we discuss the localization of light waves, mainly the challenge in observing it experimentally in 3D scattering systems. Note that in this chapter we introduce only the key topics and concepts of Anderson localization which are relevant for the study that we have conducted on identifying the signature of localization of light. The research done in the last 60 years on localization of various types of waves is immense and more detailed overviews can be found for example in [Sheng 2006, Abrahams 2010]. Moreover, in Sec. III.3 we present the numerical results based on coupled-dipole model, which verify the use of statistics of the scattered light as an observable of the localization transition for light in 3D (the part involving scalar light was done by Florent Cottier, while I did the part with vectorial light). Finally, in Sec. III.4, we further discuss the obtained results.

Firstly,Figure III. 1 :

 1 Figure III.1: Scheme illustrating multiple scattering of waves in disordered media. (a) When considering wave propagation from point A to point B, there are many possible diffusive paths.For most of these diffusive paths, the interference between probability path amplitudes vanishes on average. b) However, each multiple-scattering path has its reversed counterpart such that wave is scattered by the same set of scatterers but in the opposite order: time-reversal paths. c) In the case of closed loops, when the wave returns to the inital scattering center (A=B), the accumulated phase difference between time-reversed loops is zero and interference effects survive the ensemble average. The probability of return to the inital scattering point A is enhanced due to constructive interference on time-reversed loops. Contribution of such closed loops in wave transport in weak disorder decreases the diffusion coefficient (weak localization). By increasing the disorder, probability for closed loops increases, thus further decreasing the diffusive wave transport. At a critical amount of disorder, when the diffusion coefficient becomes zero, the phase transition from extended-wave regime to localized regime arises. d) Another way for interference of time-reversed loops to survive the configuration average is if k f = -k i , where k i is the wavevector incident on the disordered medium and k f is the outgoing wavevector. These interferences on time-reversed loops give rise to the coherent backscattering, the enhancement of the multiply scattered intensity in the backward direction by a factor of two.

  amplitudes. As the paths have different lengths, the amplitudes a p accumulate different phases: the phases of the probability amplitudes are random. Hence, the sum of the interference terms vanishes on average and the probability for propagation is determined by the incoherent-contribution term. Although this might lead to the conclusion that the interference contribution is generally irrelevant, there is an exception to such assumption which is crucial for localization of waves to occur: the time-reversed paths. For each diffusive path P consisting of N scatterers such that the wave is scattered by those scatterers in the order 1 → 2 → ... → N , there exists one reversed counterpart P consisting of the same N scatterers but with propagation in the opposite order N → ... → 2 → 1. Reversed path P can be regarded as direct path P propagating backward in time, hence the name time-reversed paths (Fig.III.1(b)). To avoid any confusion, note that throughout this chapter we are refering to the pair of direct and its reversed path as time-reversed paths. The accumulated phase difference between the time-reversed paths is given by ∆φ

1

 1 and r N are positions of initial and final scattering centers along the diffusive path. When ∆φ = 0, timereversed paths are perfectly coherent and the interference term of path amplitudes survives the configuration averaging. This occurs in two different ways: when the multiply scattered wave returns to the inital scatering center, r 1 = r N , forming a closed loop (Fig.III.1(c)), orwhen k f = -k i . As it is going to be discussed below, the former pertains to the recurrent scattering inside of a medium and it is relevant for Anderson localization. The latter, on the other hand, concerns the studies of incoming and outgoing waves, such that the wave incident on the medium is then being multiply scattered in the backward direction, and it is relevant for the phenomenon of coherent backscattering.As depicted in Fig.III.1(c), in the case of closed paths (loops), the starting point A and end point B coincide (r A = r B ). The probability to go from A to B in the case when points A and B are identical represents the probability that the wave is scattered back to the scattering center from which it originated, i.e. the return probability to the starting point A. For each closed path there exists an identical closed path which is traversed in the opposite direction such that system posses time-reversal symmetry. When propagating along those two time-reversed closed paths, waves acquire exactly the same phase and in the absence of time-reversal symmetry breaking effects, the probability amplitudes of two time-reversed loops are equal a 1 = a 2 = a. Hence, the time-reversed loops constructively interfere. Because of this non-vanishing interference between the time-reversed closed loops, the probability of

Figure III. 2 :

 2 Figure III.2: Schematic energy spectrum of a quantum particle in the disordered system. (a) If the width of energy levels ∆E is much larger than the average level spacing δW , the eigenstates are extended. (b) If δW ∆E, the eigenstates are localized.

Figure III. 3 :

 3 Figure III.3: Schematic picture of the scaling function β(g) = d lng(L) d lnLfor different dimensions as a function of the dimensionless conductance g. The scheme was constructed from known asymptotes of β(g) [Eq. (III.10)] in localized and extended regime and by assuming continuity and monotonicity of β(g). Arrows indicate the direction of variation of g as the system size L increases. Such scaling functions imply that for d = 3 there is a phase transition from extended to localize regime, which occurs at critical point g c . In lower dimensions there is no phase transition: states always localize in the thermodynamic limit L → ∞.

d = 1

 1 , 2, 3, as shown on Fig. III.3. It can be seen from Fig III.3. that the behaviour of the scaling function strongly depends on the dimensionality of the system. In d = 3 there is a critical point g c ∼ 1 given by β(g c ) = 0 at which β(g) goes through a sign change. Suppose a weak disorder such that g(L) > g c , and hence β(g) > 0. Then, as L increases, g(L) increases which makes β(g) more negative. Hence, one moves upward on the β(g) curve from the starting point g(L) > g c . The arrow on the curve in Fig. III.3 thus indicates the direction of variation of g(L) as L increases. At large enough length scales L, the asymptotic behaviour of β(g) ∝ d -2 = 1 corresponding

  of a single parameter ζ/L: g = f (ζ(η)/L).(III.11)The characteristic length scale ζ(η) is a function of the disorder strength η and it corresponds to the localization length ξ in the localized regime, while in the extended regime it can be associated to the inverse of the conductivity ∝ σ -1 . For a fixed amount of disorder (and consequently fixed ζ(η)), the domain of values for the parameter ζ(η)/L can be completely covered by varying L, hence compensating for the effect of the disorder. However, note that in 3D it would be impossible to reach the critical point g c and to crossover from extended to localized regime (and vice versa) by simply varying the sample size for a given amount of disorder, since g c is independent on L and there are two qualitatively different behaviours of β(g) depending on in which regime the system is. As stated many times, the Anderson transition takes place when the strength of disorder is sufficient for interference effects to dominate the wave transport. Hence, the disorder strength (which depends on the microscopic properties of the system) will dictate the transition and decide whether the system is on the positive or negative branch of the scaling function β(g). Moreover, the scaling theory predicts the behaviour of the localization length ξ and the conductivity σ near the critical point. Since g c does not depend on L, Eq. (III.11) implies that the correlation length ζ(η) has to diverge at the transition. The diverging behaviour of ζ(η) can be described by a power law with the critical exponent ν [Lee 1985] as ζ(η) ∝| η -η c | -ν , where η c is the so-called mobility edge, i.e. the critical disorder strength at the phase transition. Since ζ(η) ≡ ξ(η) in the localized regime and ζ(η) ∝ σ -1 in the extended regime, it follows:

Figure III. 4 :

 4 Figure III.4: Scheme illustrating different multiple scattering regimes. The Ioffe-Regel criterion states that the mean free path has to be approximately equal to λ/2π (black curve) for Anderson transition to take place. In the geometrical-optics and Rayleigh regime, scattering is always weak and k 1, hence localization is not possible. However, in the intermediate frequency regime λ ∼ a, Anderson localization can occur if the mean-free path is sufficiently reduced, which requires large refractive index contrast between the scatterers and surrounding medium m. Here, m > m .

  ], since in the presence of absorption (by possible impurities or defects in the scattering sample) the transmission decreases exponentially T ∝ exp(-L/L a ), where L a is absorption length. Therefore, it is hard to discriminate between localization and absorption using static transmission measurements since they both yield exponential decrease of the transmission with the length of the sample. The observations of claimed light localization were finally disproved in [van der Beek 2012], as time-of-flight experiments demonstrated that what they interpreted as the localization signatures can indeed be attributed to a weak absorption of light due to strains and defects in the sample resulting from grinding of GaAs. Then during 2006-2014, in several publications [Störzer 2006, Aegerter 2006, Aegerter 2007b, Aegerter 2007a,Sperling 2013,Sperling 2014], Maret et al claimed the observation of Anderson localization of visible light in 3D titania (TiO 2 ) powder samples (with high refractive indexn ≈ 2.7), where they observed deviations from diffusive-wave transport. The time-dependent measurements[Lenke 2000b] allowed them to distinguish between absorption and considered localization signatures[Aegerter 2009]: while in the absence of localization, the time-resolved transmission decay exponentially due to both scattering and absorption, deviations from this single exponential behavior points to Anderson localization[Skipetrov 2006]. Although their observations could not be attributed to the spurious absorption of light by the defects or impurities, it was questioned whether the observed localization signature was confused with nonlinear effects[Scheffold 2013, Maret 2013]. Later on, their claim on observing localization of light was revoked[Sperling 2016], as it was demonstrated that a weak fluorescence due to impurities in the disordered samples could explain the observed deviations from the diffusion model.

  e., how to unquestionably identify the photon localization such that other effects cannot be mistaken for the localization signature. The experiments in 3D disordered systems mentioned above, claimed that the localization signature was the deviation of the stationary or time-resolved light transmission from the diffusion theory. However, it turned out that the transmission measurements can only indicate localization if the spurious effects, such as fluorescence and absorption, can be measured independently and if they do not overshadow the localization effects on transmission. Therefore, we need to look for a more reliable localization signature. This is precisely what we discuss in the following sections: we demonstrate an unambiguous signature of light localization based on the statistical properties of the light, particularly the intensity fluctuations [Cottier 2019b].

III. 3 .

 3 disordered medium consisting of classical scatterers by effectively increasing the photon mean free path. Although it had been elusive whether the conclusion of absence of Anderson localization of vectorial light for resonant point scatterers can be readily extended to classical scatterers, the recent study indeed suggests so, i.e. that the near-field coupling may be responsible for the still unobserved Anderson localization of light with classical scatterers[Cobus 2021]. Moreover, recently it has been studied how to decrease the impact of longitudinal electromagnetic fields in the case of light scattering from dielectric spheres[Escalante 2017], which opens a way of observing strong localization of light with classical scatterers in the case that non-radiative transport channels prohibit it. Statistical signature of strong localization of light in3DAs mentioned, in the previous attempts to observe the light localization, the experimental observable of Anderson localization of light was assumed to be found in the decay of stationary[Wiersma 1997] or time-resolved transmission[Störzer 2006] as deviations from classical diffusion. However, static measurements of transmission are not promptly suited for observing effects of strong localization of light since they do not allow to distinguish between absorption and localization of light [van der[START_REF] Van Der | [END_REF]. Time-resolved measurements allow that. However, deviations from classical diffusion in time-resolved measurements can be a signature of strong localization only if there are no other spurious effects emerging that could explain the deviations and overshadow the localization, such as the fluorescence due to impurities in the disordered sample[Sperling 2016].

  disordered sample. Prompted by that research, we numerically studied the speckle statistics of light waves in the framework of cold atoms as disordered medium. Here we will present our results that indicate that the unambiguous signature of localization of light in 3D can be indeed identified in the large intensity fluctuations. Such statistical signature of light localization could be relatively easily measured in the experiments and it does not require the exhaustive scaling analysis.

  the light intensity scattered by an ensemble of atoms. Here we will briefly reintroduce the coupled-dipole model and describe the numerical procedure, as well as the atomic sample (Fig.III.5).We consider an ensemble of N identical, motionless, point-like atoms located at random positions r j , j = 1, ..., N , within a given sample geometry and driven by an external monochromatic electric field. Like in the previous chapter, we consider the low-intensity limit of the driving field, i.e., linear optics regime. In the vectorial model, we consider 4-level atoms with the ground state |g = |J g = 0, m g = 0 and the triple-degenerate 6 excited state |e with substates |e me = |J e = 1, m e = 0, ±1 , coupled to the electric field by the transition dipole moment. The substates of excited state |J e = 1, m e = 0, ±1 correspond to the projection m e = 0, ±1 of total angular momentum J e on quantization axis z. On the other hand, in the scalar model the polarization of the electric field is discarded and atoms are considered to be "pure" two-level point scatterers without Zeeman internal structure. The scalar coupled-dipole equations (CDEs) for the amplitude of the atomic dipole β j are given by

  the above scalar and vectorial CDEs, j, m ∈ 1, N , d is the electric-dipole transition matrix element considered to be equal for every atom and every transition, Γ 0 is the decay rate of an isolated, individual atom, ∆ 0 = ω L -ω 0 is the detuning of the frequency of the incident electric field ω L with respect to the atomic resonance ω 0 , k L is the wavevector of the incident field with k L ≈ ω 0 /c = k 0 and r jm = r jm rjm = r jr m is the distance vector between atoms j and m. In the vectorial CDEs, δ ζ,η is Kronecker delta, rζ = êζ • r is the component of the unit vector r = r/r along the direction ζ = 0, ±1 and ζ, η ∈ (±1, 0) are the spherical-basis vector components whose unit vectors are ê±1 = ∓1/ √ 2(ê x ± iê y ), and ê0 = êz .

Figure III. 5 :

 5 Figure III.5: Scheme of the numerical model. Atoms are considered to be point-like dipoles, randomly positioned within the given 3D sample geometry with uniform density. In the scalar model, the atoms are treated as two-level resonant scatterers with resonance frequency ω 0 . Contrary to the scalar model, the internal (Zeeman) level structure and polarization of light are considered in the vectorial model, where we assume ground state |J g = 0 and excited state |J e = 1 of atoms.The excited state is triple degenerate with substates |J e = 1, m e = 0, ±1 . Atoms are driven by an external electric field which is considered to be a Gaussian beam with waist smaller than the size of the sample. The frequency of the external electric field ω L is detuned by ∆ 0 = ω L -ω 0 from the atomic resonant frequency ω 0 . Our focus is on the statistics of the transmitted intensity, which has a characteristic speckle pattern. We detect the intensity radiated in a given direction in the far-field limit outside of the sample. To obtain the variance of the transmitted intensity we average the intensity over many configurations of atomic positions.

  III.28) in the vectorial model, where nζ = êζ • n the component of the unit vector n = R/R along the direction ζ = 0, ±1. Note that the first term in the above expressions corresponds to the incident electric field, while the second one to the electric field radiated by the atomic dipoles in the sample.As we want to avoid strong single-scattering contribution from the boundary of the sample (strong single-scattering, which obeys Rayleigh statistics, can overshadow the effect of localization on transmitted light statistics), instead of using a plane wave, we use a Gaussian beam whose waist w 0 is smaller than the size of the sample L: w 0 = L/4. For a monochromatic Gaussian beam propagating along the z axis (k L = k L ẑ), the electric field amplitude at R = Rn = X x + Y ŷ + Zẑ position is given by

  e. in the absence of Anderson localization, the extracted conductance yields arbitrarily large values g 1, corresponding to the divergent number of active propagating optical modes, i.e. diffusive channels (Sec. III.3.1). However, for higher densities and particular laser detunings (the detuning dependence will be discussed below) small values of g, close to unity (down to g = 0.27 ± 0.02 as for ρ = 44/λ 3 , ∆ 0 = 1Γ 0 ), are obtained, implying the emergence of light localization. Prompted by the research based on eigenvalue scaling analysis which identified the phase transition in a given range of densities (ρλ 3 20) and laser frequencies (∆ 0 ∼ 1Γ 0 ) [Skipetrov 2018c], we monitor the dimensionless conductance g and the normalized intensity variance σ 2 Î in the range of densities ρλ 3 = [5 : 44] and laser detunings ∆ 0 = [-2, 2]Γ 0 . As presented in Fig. III.7(b), the conductance collapses fromvery large values in the diffusive regime to low ones, close to unity, in a V-shaped area corresponding to ∆ 0 ∼ 1Γ 0 and ρ 20/λ 3 . This V-shaped area, i.e. its corresponding ρλ3 and ∆ 0 for which we observe the collapse to small values of g, thus implying the localization, is in very good agreement with the localization phase diagram obtained from the eigenvalue scaling analysis (see Fig.3in Ref.[Skipetrov 2018c] where the mobility edges are identified).Moreover, in the diffusive regime for which the Rayleigh intensity distribution was observed with large g 1, we obtain σ 2 Î ≈ 1. On the other hand, the collapse of the dimensionless conductance to small values close to unity is accompanied by a significant increase of the intensity variance [see Fig. III.7(a)]: the parameters (ρλ 3 , ∆ 0 ) for which low values of g are found, correspond to increased intensity fluctuations, i.e. strong deviations of σ 2 Î from unity, as expected (for example, the obtained intensity distribution presented in Fig.III.6(b) with g = 0.27 corresponds to σ 2 Î ≈ 6.8). The V-shaped area, for which we observe low values of g and correspondingly strong deviations of σ 2 Î from unity, corresponds to the Ioffe-Regel criterion for localization k 0 < k 0 c ∼ 1

Figure III. 6 :

 6 Figure III.6: Probability distribution function of the normalized intensity Î = I/ I in (a) the diffusive regime (ρ = 5/λ 3 , ∆ 0 = 0, g 20) and (b) the localized regime (ρ = 44/λ 3 , ∆ 0 = 1Γ 0 , g = 0.27 ± 0.02). The black dashed curves correspond to Rayleigh intensity distribution [Eq. (III.20)], the blue crosses correspond to the numerical results obtained from Eq.(III.27) for 10 4 realizations and the red dash-dotted curve was computed from Eq. (III.23). Intensity is detected at the angle (θ, ϕ) = (75 • , 0 • ) from the laser propagation axis, from a system of size kL = 32.4 and (a) N = 684 or (b) N = 6066. The results were obtained by Florent Cottier [Cottier 2019a].

Figure

  Figure III.7: Phase diagram for the (a) variance of the intensity σ 2 Î and (b) conductance in the (∆ 0 /Γ 0 , ρλ 3 ) plane. Simulations are realized for a homogeneous cubic cloud of side length kL = 32.4, using 800 realizations and an observation angle θ = 75 • . The value of the conductance is saturated at the arbitrary value of 20, as it diverges for non-localized samples. The black curve corresponds to Eq. (III.30) for α = 0.5. The results were obtained by Florent Cottier [Cottier 2019a].

  the threshold k c = 0.5 provides a good approximation of the critical region where we observe large intensity variance (see Fig. III.7(a)), confirming the validity of the Ioffe-Regel criterion (α = 0.5 is of order unity) as an indicator of the localization transition in our system. This also gives an insight into the existence of two mobility edges (laser frequencies) for a given density where the localization transition takes place: for a given ρλ 3 , there are two laser frequencies for which the same critical photon mean free path is obtained [Eq. (III.30)].

  Fig. III.8. As can be seen, σ 2Î ≈ 1 for all considered detuning and densities: we do not obtain significant deviations of intensity variance from unity, as is the case with the scalar model (Fig.

  Figure III.8: Variance of the normalized transmitted intensity σ 2 Î as a function of the density of the cylindrical atomic sample (of length L and diameter L) ρλ 3 and the laser detuning ∆ 0 = ω L -ω 0 in the vectorial model, where we consider |J g = 0 ground state and three-fold degenerate |J e = 1, m e = 0, ±1 excited state for atoms. The size of the sample is k 0 L = 21.5. The ground state is coupled to the excited state by the circularly σ + polarized laser beam (E L ⊥ ẑ; ˆ L = ê-1 ) whose wavevector is in ẑ direction (k L ẑ). The axis of cylindrical atomic sample is parallel to the propagation direction of the incident laser beam ẑ. The radiated intensity is collected at θ obs = 75 • and we average over 400 positional realizations and 51 values of azimuthal angle φ obs = [0 • , 360 • . We do not observe strong fluctuations of the transmitted intensity, i.e. significant deviations of the intensity variance from unity, as it was the case in the scalar model (Fig. III.7), implying the absence of localization of light with the vectorial model, in agreement with theoretical studies involving scaling analysis [Bellando 2014, Skipetrov 2014].

Figure III. 9 :

 9 Figure III.9: Complex eigenvalue distribution of coupled-dipole modes, λ n = i∆ n -Γ n /2, obtained from the scalar model [Eq. II.19] (a,b) and vectorial model [Eq. II.14] (c,d) with a cylindrical sample of ρλ 3 = 5, N = 1260 (k 0 L = 43) (a,c) and ρλ 3 = 40, N = 1260 (k 0 L = 21.5) (b,d). The imaginary part of the mode eigenvalue corresponds to the mode frequency shift ∆ n = ω 0 -ω n , where ω n is the mode eigenfrequency, while the real part corresponds to the mode decay rate Γ n = 1/τ n , which is the inverse of the mode lifetime τ n . The number of realizations in the scalar case is N r = 60, while in the vectorial case N r = 20. The color code corresponds to the inverse participation ratio (IPR). The black solid lines correspond to the pair branches obtained from the analytical expression of pair (two atoms positioned very closely) eigenvalues, given by Eq. (II.23) for scalar and Eqs. (II.21) and (II.22) for vectorial model with ∆ 0 = 0. The large IPR (red/dark red) characterizes the localized modes and modes of close pairs of atoms (super-and subradiant pairs), while the low IPR (yellow) corresponds to the collective (extended) modes of N n1 atoms. For high enough densities in the scalar model, in the bottom of the distribution a tail of modes of very long lifetimes with large IPR emerges which corresponds to the localized modes. Localized modes with large IPR and very long lifetimes that appear only at sufficiently high densities are to be discriminated from subradiant pairs which are present at any density. On the contrary, in the vectorial model such localized states do not appear for higher densities in the bottom of distribution, signaling the absence of light localization[Bellando 2014, Skipetrov 2014].

  III.31) where G ζ,η (r jm ) is given by Eq. (III.26b). Note that the difference between these coupleddipole equations and those without applied magnetic field [Eq. (III.26a)] is the Zeeman shift that appears in the first term of the equation above with ∆ 0 -m e ∆ B = ω L -ω me = ∆ me which signifies the detuning of the laser from the resonant frequency of the corresponding substate m e . Let us first have a look at the eigenvalue distribution of coupled-dipole modes in the vectorial model with applied magnetic field to have a better insight into what to expect when studying intensity fluctuations. Similarly to the studies in Sec. II.1.2, we consider the 3N × 3N matrix G 7 , whose elements are

  Figure III.10: Eigenvalue distribution obtained by diagonalizing the G matrix given by Eq. (III.32) for a cylindrical atomic sample of size k 0 L = 21.5 and uniform density ρλ 3 = 5 in the case of degenerate excited state (no magnetic field is applied; ∆ B = 0) (a) and in the case of atoms subjected to a strong magnetic field that produces a Zeeman shift ∆ B = 1000Γ 0 (b-e). When the strong magnetic field is applied, the eigenvalues are well split into three groups centered around ω 0 -ω n = -m e ∆ B , thus each group corresponding to a transition to a given m e substate (b). Panels (c-e) are the close-ups of panel (b) on an eigenvalue group corresponding to m e = +1 (c), m e = 0 (d) and m e = -1 (e). The dashed black lines on panel (a) corresponds to the analytical solution of subradiant and superradiant pair eigenvalues given by Eqs. (II.21) and (II.22). When the strong magnetic field is applied, the eigenvalues associated to the close pairs of atoms are as well split into three groups corresponding to a given transition m e = 0, ±1. Contrary to the case of degenerate atoms, in the case of applied magnetic field, the pairs eigenvalues depend on the angle θ 12 between the axis along which the atomic pairs are positioned and the axis given by the direction of applied magnetic field ẑ. The gray and black solid lines on panels (c-e) correspond to the pairs eigenvalues for θ 12 = 0 • and θ 12 = 90 • , respectively, obtained by diagonalizing the G matrix [Eq. III.32] for N = 2 atoms and their varying separation distance with θ 12 = 0 • and θ 12 = 90 • . For all other angles θ 12 , the eigenvalues fall in between those two lines. The colorcode corresponds to the inverse participation ratio (IPR). The eigenvalues are plotted for N r = 300 realizations of atomic positions.

Figure III. 11 :Figure

 11 Figure III.11: Same as Fig. III.10 but with density ρλ 3 = 50 and N r = 30 realizations of atomic positions.. With the applied magnetic field, the long-lived tail containing modes with large IPR appears in the bottom of the eigenvalue distribution for m e = ±1 (c,e), but not for m e = 0 (d). These localized modes emerging for m e = ±1 at high densities in the narrow frequency window with very long lifetimes are to be distinguished from subradiant pair modes with large IPR that are present in the distribution group of every m e for any density (see Fig. III.10). The insets in panels (c) and (e) are close-ups of the long-lived tail that appears for m e = ±1 when the magnetic field is applied.

  III.13 where we show for each atomic dipole i = 1, ..., N the steady state amplitudes |β ζ i | 2 for ζ = m e = 0, ±1 which characterize the probability of the excitation of atom i to substate |J e = 1, m e = ζ . On Fig. III.13(a) we show those amplitudes in the case when no magnetic

Figure III. 13 :

 13 Figure III.13:The dipole amplitude squared |β me i | 2 of an atomic dipole i = 1, ..., N and a given transition to m e = 0, ±1 for k 0 L = 21.5 and ρλ 3 = 50 obtained from Eq. (III.31) in steady-state. For atomic dipoles with degenerate excited state (∆ B = 0), |β me i | 2 for each m e = 0, ±1 have approximately equal non-vanishing values (the external driving field is circularly σ + polarized with detuning ∆ 0 = 0), since excitation transfer among dipoles can happen through any transition corresponding to m e = 0, ±1 (a). When the strong magnetic field is applied (∆ B = 1000Γ 0 Γ 0 ), and the external electric field drives the transition to a given m e substate, the photons can be rescattered only through that particular transition m e due to the huge Zeeman shift; hence only the dipole amplitude corresponding to that transition m e has a finite non-vanishing value, while the dipole amplitudes associated to other transitions are negligible, implying that the atoms subjected to a strong magnetic field act as a non-degenerate two-level system (b-d). Panels (b,c,d) correspond to the case of external electric field driving the transition to m e = +1, -1, 0, respectively, with detuning ∆ me = 0. The atomic sample-external field configuration for a given m e transition is the same as those depicted on left side of Fig. III.12. The dipole amplitudes of each dipole are averaged over 20 realizations.

  , the dipole-dipole coupling via radiated electric field is achieved only through one specific m e transition. For example, on Fig. III.13(b) we show the case when the atomic transition to m e = +1 is driven by the external electric field with detuning ∆ +1 = 0 and σ + polarization. Then, due to the large Zeeman shift of substates, only the transition to m e = +1 can be coupled by the electric field radiated by dipoles (or in terms of multiple scattering, the photon can be multiply scattered only from the given m e = +1 transition). Hence, only the component m e = +1 of dipole amplitudes |β me i | 2 ∼ 1 while other two are negligibly small, |β m e i | 2 |β me i | 2 with m e = m e . The same holds when the transition to m e = 0 or m e = 1 is driven by an external electric field [ see Fig.

  r) 2 with ζ = m e . (III.33b)Although the form of these CDEs looks very much alike to scalar CDEs given by Eq. (III.25a) as we got rid of polarization (transition) coupling terms, it is important to note that there is still a significant difference: the near-field interaction (∝ 1/r 2 , 1/r 3 ) encoded in the Green's matrix G , which is absent in the scalar model[Eq. (III.25b)]. Let us now rewrite the Green's matrix G by using r±1 = ∓ sin θ √ 2 exp(±iφ) and r0 = cos θ:

III. 4 .

 4 photon to escape from the system. Our observations, both for the scalar and vectorial models, made on the fluctuations of the transmitted intensity are in agreement with the results based on scaling theory with the Thouless number obtained from the modes eigenvalues, where it has been demonstrated that the light indeed localizes when the scalar model is considered, while the localization is absent in the vectorial model with degenerate Zeeman sublevels, and moreover that the localization is resurrected for m e = ±1 when the strong magnetic field

  the difference between the strong intensity fluctuations obtained from the scalar model and from the vectorial model: in the scalar model the strong intensity fluctuations are observed only for the mobility edges, i.e. at the localization phase transition (Fig. III.7), while in the vectorial model with applied magnetic field we also observe strong intensity fluctuations inside of the localized regime (Fig. III.12). On one hand, it is not surprising that the strong intensity fluctuations appear at the localization phase transition as strong fluctuations are characteristic for other types of phase transitions. On the other hand, if it were only a characteristic of the phase transition, why are the strong intensity fluctuations identified as an indicator of the strong localization in quasi-1D systems[Chabanov 2000] for which there is essentially no phase transition (for ξ < L, the system is always localized in 1D)? In our case of the scalar model, one important aspect is that the obtained dimensionless conductance g inside of the localized regime has large values g > 1 [Fig.III.7(b)] (let us recall that g essentially characterizes the number of diffusive channels and the localization transition is expected to take place when g = 1, hence one would expect to observe g < 1 inside of the localized regime). Note that the dimensionless conductance has been obtained from the fit of the numerically obtained intensity distribution to Eq. (III.23). However, the intensity distribution function given by Eq. (III.23) has been obtained for the diffusive regime by taking into account corrections due to correlations among loopless scattering paths, without considering paths forming a loop, which are important for the localization phase transition and the localized regime itself (for more details on this intensity distribution function see, for example, Ref.

  cold atoms that are motionless over the timescale of the measurement. The most promising candidates for observing light localization are Strontium and Ytterbium cold atoms, which have the characteristic |J g = 1, m g = 0 → |J e = 1, m e = 0, ±1 transition.Finally, let us discuss the role of near-field dipole-dipole interactions on the light localization and its connection to the near-field interaction effects on other physical phenomena. It has been discussed that the near-field dipole-dipole interactions introduce additional transport channels (non-radiative) that prevent the localization of photon in the strongly disordered system[Skipetrov 2014], as the opening of these non-radiative transport channels essentially increases the photon scattering mean-free path. Hence, the suppression of the near-field interactions, and with that of the opening of non-radiative transport channels, by applying a strong magnetic field, recovers the light localization when the full vectorial model is considered. Furthermore, recently the detrimental role of near-field interactions for the maximum obtainable index of refraction of optical materials has been identified: the nearfield interaction among randomly distributed atoms yields an inhomogeneous broadening of the atomic resonance frequencies, which puts a limitation on the maximum obtainable refractive index[Andreoli 2021]. Moreover, we have demonstrated in the previous chapter the detrimental role of near-field interactions on subradiance: above a critical density of the atomic sample, the scaling of subradiance with the on-resonance optical thickness b 0 breaks down and the subradiant lifetime for a given b 0 decreases with increasing density (see Figs.II.4 and II.5). We have attributed this to the inhomogeneous broadening due to near-field interactions, as we have observed in the eigenvalue distribution that with increasing density for a given b 0 the longest-lived subradiant modes disappear while the spread of eigenvalues along the frequency axis becomes broader (see Fig. II.14). In fact, the plausible explanation is that the inhomogeneous broadening stemming from near-field interactions results in the effective decrease of b 0 (and increase of the photon mean-free path) and it modifies b 0 in a more stronger fashion as the density increases, yielding the observed decrease of subradiant lifetimes. This interpretation of the inhomogeneous broadening originating from near-field interactions could be extended to the problem of Anderson localization of vectorial light in terms of increased scattering mean-free path as the result of the broadening. In fact, we can see in the eigenvalue distribution shown on Fig. III.11 that the spread of eigenvalues

  identified for Anderson localization of light waves: the scalar coupled-dipole model predicts Anderson localization of light waves, but not the vectorial one[Skipetrov 2014], since the near-field terms open up additional non-radiative transport channels. However, adding a strong magnetic field allows one to recover localized modes for vectorial light[START_REF] Skipetrov | Magnetic-Field-Driven Localization of Light in a Cold-Atom Gas[END_REF]].In Chapter III we numerically studied a possible signature of Anderson localization of light in cold atoms, based on the measurement of intensity fluctuations of the scattered light. We have shown that at the transition to the localized regime, the fluctuations are enhanced. In particular, for vectorial light, we have demonstrated the absence and recovery of enhanced intensity fluctuations without and with the applied magnetic field, which further corroborates that the intensity fluctuations could serve as an experimental way to probe the localized states.

a

  source of intensity fluctuation of the scattered light [Eloy 2018], which might complicate the analysis. Nevertheless, preliminary computations by R. Bachelard with atomic motion have shown that the enhanced fluctuations at the transition in the scattered intensity are still present.

Note that the numerical model represented in the following two chapters is based on this mean-field model with the additional assumption of the low excitation limit (linear-optics regime) such that the ground-state population is considered to be unchanged

Note that in the rotating frame of the driving field with ∆ 0 = 0, the frequency shifts are ∆ pair,-=ω L -ω pair,-≈ -3Γ 0 2(k 0 r 12 ) 3 + ∆ 0 and ∆ ⊥ pair,-= ω L -ω ⊥ pair,-≈ 3Γ 0

4(k 0 r 12 ) 3 + ∆ 0 .

For example, the participation ratio (or the inverse participation ratio), which tell us how many scatterers are involved in the eigenstate could also be used to discriminate between localized and extended state, since in the localized state only few scatterers are involved (it is a local phenomenon), while for the extended state many scatterers contribute. In the most limiting case, when the state is localized on one scatterer, the participation ratio is equal to unity. On the other hand, for the extended state the participation ratio is equal to N , when all N scatterers are involved in the state.

Here we included the weak localization correction going as a/g (a > 0 is a constant) to the scaling function for g g c , since β(g) = d -2 = 0 for d = 2 (see Ref.[Lee 1985]).

Note that there is more to the story in terms of light polarization: the observed CBS enhancement depends also on the polarization channel of the outgoing light. For more details on that, see one of the mentioned references.

Positional correlations (short range order) can also lead to the enhancement of the effective scattering cross section as it has been showed in experiments in colloidal liquids[Rojas-Ochoa 2004b]

The discussion of the vectorial coupled-dipole model when the degeneracy of the excited state is lifted by applying a strong magnetic field is in Sec. III.3.5.

In comparison to Eq. (II.14), note that here we have the additional -im e ∆ B δ ζ,me term (corresponding to the Zeeman shift of substates) for the diagonal elements of the G matrix and we dropped the laser detuning term ∆ 0 , i.e. we consider the frequency shifts of modes in respect to the atomic resonance: ∆ n = ω 0 -ω n .

Since in Eq. (III.33a) G (r ij ) with r ij = r i -r j , θ then corresponds to the angle between rij and z-axis.

Acknowledgements

this work. I am deeply grateful to him for welcoming me to Brazil, for wonderful collaboration throughout the years of my PhD research, for many fruitful discussions, as well as for his support and advice.

the eigenvalues are split into three well-separated groups [Fig. III.10(b) and Fig. III.11(b)].

Those eigenvalue groups are concentrated around ω 0 -ω n ≈ -m e ∆ B . In other words, each of the eigenvalue groups is associated with the transition to the corresponding m e = 0, ±1 substate of excited state. Therefore, by applying the strong magnetic field the degeneracy of excited state is lifted such that the Zeeman substates m e = 0, ±1 are well split (∆ B Γ 0 ) and atoms effectively behave as a non-degenerate two level system (by assuming that the external electric field couples the transition to a particular m e substate, the coupling of atomic dipoles via the electric field produced by all other dipoles can be achieved only through the given m e transition while the coupling through other transitions is negligible, since the other substates are "unreachable" for the photon scattered from a given m e transition due to large Zeeman shift). Furthermore, we observe that in the low-density case, the distributions of Let us now focus on intensity fluctuations with vectorial model when the strong magnetic field is applied. We use a cylindrical atomic sample (diameter L and length L) whose axis is along the direction of propagation of external electric field and size is given by k 0 L = 21.5. We solve the coupled-dipole equations given by Eq. (III.31) for β ζ j in steady state, with ∆ B = 1000Γ 0 . Then from obtained steady-state dipole amplitudes we calculate the radiated intensity in the far-field (R = 250L) given by Eq. (III.28). We study the variance of the intensity for