Spatial correlations of the spontaneous decay rate as a probe of dense and correlated disordered materials

We study theoretically and numerically a new kind of spatial correlation for waves in disordered media. We define $C_{\Gamma}$ as the correlation function of the fluorescent decay rate of an emitter at two different positions inside the medium. We show that the amplitude and the width of $C_{\Gamma}$ provide decoupled information on the structural correlation of the disordered medium and on the local environment of the emitter. This result may stimulate the emergence of new imaging and sensing modalities in complex media.


I. INTRODUCTION
Sensing and imaging are key applications of the study of light propagation in strongly scattering media [1]. Optical coherence tomography [2] is one of the most emblematic example but is limited to small optical thicknesses where the single scattering regime takes place. Because of multiply scattered light, sensing and imaging deeply inside a strongly disordered system is very challenging and has been a matter of intense study in the last two decades. Important breakthroughs were achieved recently by "learning" the system using wavefront shaping techniques [3,4], by using multimodal approaches such acousto-optics [5], or by taking advantage of particular features of light scattering in complex environments such as the memory effect [6], to cite a few examples.
Another possibility consists in using fluorescent emitters embedded inside the scattering medium. It is well known that the spontaneous decay rate Γ of such an emitter strongly depends on the local environment [7]. More precisely, this decay rate is proportional to the Local Density of States ρ (LDOS) at the position of the emitter [8,9]. This makes this quantity strongly non-universal which is of great interest in terms of imaging [10], sensing [11][12][13] and control [14]. By performing statistics, signatures of the local order around the emitter [15] and of transport regimes [16] are revealed.
Interestingly, the fluctuations of the LDOS are encoded in the spatial intensity correlation function (speckle correlation) measured outside the medium . More precisely, LDOS fluctuations generate an infinite-range contribution to the speckle correlation function denoted by C 0 [15,17,18]. This contribution is a feature of speckle patterns produced by a point source (e.g. a fluorescent emitter) located inside the medium. Measuring C 0 amounts to measuring LDOS fluctuations. In the optical regime, this can be achieved by measuring fluctuations of the spontaneous decay rate of fluorescent emitters [19][20][21]. In acoustics, direct measurements of C 0 from speckle * Electronic address: remi.carminati@espci.fr correlations have been reported [22].
In this paper, we introduce and study a new type of spatial correlation function, denoted C Γ , and defined as the correlation function of the spontaneous decay rate of a single emitter measured at two different positions inside the disordered medium. As will be shown, this correlation function generalizes the usual C 0 contribution. We demonstrate that the amplitude and the width of C Γ provide decoupled information on the structural correlation of the disordered medium and on the local environment of the emitter, which makes this correlation function particularly interesting for sensing and imaging in complex media.

II. DECAY RATE STATISTICS AND SPECKLE CORRELATIONS
The normalized correlation function of the intensity measured in the far-field is defined as where I(u) is the intensity in direction u and . . . denotes an average over all possible configurations of disorder. This correlation function is usually splitted into three components [23] C(u, The C 1 term is usually the predominant short-range term and gives typically the size of the speckle spot. C 2 and C 3 are long-range terms with smaller amplitudes. When the speckle pattern is produced by a point source embedded inside the scattering medium, an additionnal term of infinite-range exists, and is denoted by C 0 [17]. The C 0 contribution to the correlation function is related to the normalized fluctuations of the LDOS ρ(r 0 ) [15,18] where r 0 is the position of the emitter. For a fluorescent emitter in the weak-coupling regime, the spontaneous decay rate Γ(r 0 ) is proportionnal to the LDOS [9], and C 0 can be rewritten as Several studies of LDOS fluctuations in disordered media, or equivalently of C 0 , have been reported [9,11,16,19,20,[24][25][26]. It was shown that C 0 originates from near-field interactions with the nearby scatterers, providing a non-universal behavior that is particularly relevant for sensing and imaging [15,19]. Moreover, C 0 is also expected to be influenced by structural correlations in the disorder [13,15]. In this article, we show that the new correlation C Γ carries enough information to extract signatures of both the near-field interactions and the structural correlations, thus providing a potentially useful extension of the usual C 0 correlation function.

A. Methodology
To get insight into the behavior of the new correlation function, we begin with a numerical study. The system of interest is depicted in Fig. 1. N point-dipole scatterers are lying between two coaxial cylinders of radii R 0 and R respectively and of longitudinal size 2R. The inner region with radius R 0 corresponds to the region within the medium in which the fluorescent source is free to move. The optical properties of the scatterers are described by an electric polarizability where ω is the emission frequency, ω 0 the resonant frequency of the scatterers, γ the linewidth and c the speed of light in vacuum. From the polarizability α, we can compute the scattering (σ s ) and the extinction (σ e ) crosssections of one scatterer. They are given by where k 0 = ω/c = 2π/λ. The optical theorem is correctly fulfilled by the polarizability model, and in a nonabsorbing medium such as the one considered here, we have σ e = σ s . Defining the density of scatterers by is the volume of the scattering system, we have also access to the scattering mean-free path ℓ B in the limit of an uncorrelated system (Boltzmann mean-free path). Its expression is To generate disorder correlations, a fictitious exclusion volume of diameter a is forced between scatterers. This mimic a hard sphere potential. By increasing the value of a, one increases the correlation level. Instead of using the parameter a to characterize the level of structural correlation, we use the effective volume fraction f defined as where V 0 = πa 3 /6 is the exclusion volume around each scatterer. Note that the effective volume fraction f has to be understood as a correlation parameter, that is changed by changing V 0 only (the real density of scatterers N is constant throughout the study).
PSfrag replacements 1: Sketch of the system. The strongly scattering medium lies between two cylinders of radii R0 and R, respectively, and of length 2R. The inner region with radius R0 corresponds to the region within the medium in which the fluorescent source is free to move, from position r0 to position r ′ 0 . To mimic hard spheres correlations, a minimum distance a is forced between scatterers.
The emitter lies initially at r 0 , the center position (see Fig. 1), that can be changed to another position r ′ 0 along the cylinder axis. The distance R 0 corresponds to the minimum distance forced between the source and all scatterers. In other words, it parameterizes the near-field environment of the source (proximity of scatterers). The spatial correlation function C Γ studied in this paper is defined as where ∆ = |r 0 − r ′ 0 |. For ∆ = 0, this expression coincides with the definition of C 0 in Eq. (4).
To compute C Γ , we have first to solve Maxwell's equations for a point dipole source. For that purpose, we use the coupled dipoles method [27]. It consists in calculating first the exciting field on each scatterer given by a set of N linear equations: where µ 0 is the vacuum permeability and p the source dipole. G 0 is the Green function in vacuum. For vector waves in three dimensions, it is given by where PV, I, ⊗, δ are the Cauchy principal value operator, the identity tensor, the tensor product operator and the Dirac delta function respectively. We have used the notations R = r − r 0 and R = |R|. Once the exciting fields on each scatterer are known, the field at any position can be computed using By varying the orientation of the source dipole p, the Green function of the full system G can be obtained from the relation E(r) = µ 0 ω 2 G(r, r 0 , ω)p. The LDOS averaged over all orientations of the source dipole can be deduced by [9] ρ(r 0 , ω) = 2ω where Tr denotes the trace of a tensor. We usually prefer to deal with normalized quantities. Defining the vacuum LDOS as ρ 0 = ω 2 / π 2 c 3 and the vacuum decay rate of the emitter by Γ 0 , we have Repeating the operation for another position r ′ 0 of the source and averaging over disorder configurations leads to an estimate of C Γ .

B. Numercial results
We have performed numerical simulations on a system with parameters such that k 0 ℓ B = 19 (strength of the disorder) and b B = 2R/ℓ B = 1.25 (optical thickness). The parameters are given in the caption of Fig. 2. Thus the numerical simulations are performed in a dilute system and close to the single-scattering regime. The simulations are performed by varying the level of structural correlation of the disorder, measured by the correlation parameter f . The results are shown in Fig. 2 (a). We clearly see that the width of the curves is independent on f , suggesting that the width depends essentially on the microscopic length scale R 0 that measures the proximity of scatterers around the emitter. The dependence of f is encoded in the amplitude of the correlation function C Γ for ∆ = 0, as shown in Fig. 2 (b). Note that this amplitude corresponds to C 0 , that is known to depend on f [15].
For weak structural correlations of the disorder (small values of f ), the exclusion distance a between scatterers is small, and more than one scatterer can lie in the near field of the emitter. This implies that the emitter can interact with two or more scatterers, inducing a strong dependence of the amplitude of C Γ on the correlation parameter f . For a high level of structural correlations (large values of f ), the exclusion distance a is large enough to exclude the possibility of interaction with more than one scatterer. For that reason, the amplitude of C Γ is almost independent on f in this regime. This qualitatively explains the shape of the curve in Fig. 2 (b).

IV. ANALYTICAL THEORY
To get physical insight, we support the numerical data by a theoretical analysis. This has also the advantage to provide simple analytical formulas that could be useful in practice. As the optical thickness b B is close to unity, the system operates in the single scattering regime. In that case, the correlation function can be computed analytically, at the price of a few crude but nevertheless controlled approximations. In the single-scattering regime, the decay rate is given by Using the expression of the vacuum Green tensor [Eq. (11)], we find where R j = |r j − r 0 |. As ∆ = |r 0 − r ′ 0 | is on the order of R 0 ≪ λ, we consider that the most important contribution is given by the scatterers lying in the near field of the emitter. Thus we now consider a subset Ω of scatterers located in the vicinity of the source inside a volume V ′ . This subset is defined by and under this near-field approximation, the decay rate becomes where α ′′ (ω) = Im α(ω). The computation of the correlation C Γ [Eq. (9)] requires the computation of the first two statistical moments of the decay rate: where P ({r j }) is the probability density of having the scatterers at positions {r j }. We denote by N ′ the average number of scatterers in Ω, a quantity that depends essentially on the exclusion radius R 0 around the emitter and on the correlation parameter f . Using this notation, the average decay rate is given by where P (r j ) is the probability density of finding one scatterer at position r j . The integral involves a fast decaying function in space, meaning that the integration volume V ′ can be replaced by V without changing the result. Then P (r j ) = V −1 and we obtain The second moment can be obtained in a similar way, leading to where R ′ k = |r k − r ′ 0 | and P (r j , r k ) is the probability density of having two scatterers at positions r j and r k . It is given by with h the pair correlation function. This leads to from which the following expression of the correlation function C Γ is readily deduced: To compute the integrals analytically, we consider infinite cylinders (i.e. R ≫ ∆) and a small radius for the inner cylinder (i.e. R 0 ≪ R). We obtain 0 , the second term in Eq. (25) can be neglected and the correlation C Γ reduces to the first and the last terms. The bulk pair correlation function is considered such that it only depends on the distance between the two points r and r ′ . We also consider a small correlation level (f ≃ 0.1) such that the pair correlation function can be approximated by (for higher structural correlation levels, a refined model is required [28]) This leads to where V 0 (r) is the exclusion volume around the scatterer centered at position r. For dilute media, the quantity 1/|r ′ −r ′ 0 | 6 is slowly varying and can be replaced by 1/|r− r ′ 0 | 6 , so that Eq. (28) reduces to Finally, the correlation function C Γ is given by This expression provides a theoretical basis to the qualitative discussion presented at the end of section III. First, it shows that the width of the correlation function C Γ depends only on the exclusion radius R 0 , i.e. on the minimum distance between the fluorescent source and the nearest scatterers. More precisely, the full width at half maxium can be approximated by ℓ 1/2 = 0.89R 0 . Second, as far as the amplitude is concerned, two important regimes can be identified. For large values of f (i.e. f > 0.02), the source is chiefly interacting with one scatterer, so that N ′ = 1. From Eq. (30), one sees that the amplitude of C Γ does not depend on f in this regime, as already seen in Fig. 2 (b). Conversely, for small values of f , the source can interact with more than one scatterer, and the amplitude of C Γ depends on f through both N ′ and V 0 . It is interesting to compare precisely the numerical results with the approximate analytical model. The comparison is shown in Fig. 3. It can be seen that Eq. (30) describes very well the dependance of the correlation function C Γ on ∆ for all values of f , showing that the analytical model provides a very accurate description of C Γ .

V. CONCLUSION
In summary, we have highlighted a new type of spatial correlation function based on the decay rate of a fluorescent emitter measured at two different positions inside a disordered medium. A numerical and analytical study has revealed that this correlation function contains more information than the usual C 0 intensity correlation function. In particular, by measuring its width and amplitude, it is possible to decouple the effect of the nearfield interactions (proximity effects) and of the structural correlation of disorder. This opens new perspectives for imaging and sensing in complex media, and in particular in correlated media whose interest in photonics is growing up.