# Light in correlated disordered media

Kevin Vynck<EMAIL_ADDRESS>Univ. Bordeaux, Institut d’Optique
Graduate School, CNRS, Laboratoire Photonique Numérique et Nanosciences
(LP2N), F-33400 Talence, France  Univ. Claude Bernard Lyon 1, CNRS, Institut
Lumière Matière (iLM), F-69622 Villeurbanne, France  Romain Pierrat  ESPCI
Paris, PSL University, CNRS, Institut Langevin, F-75005 Paris, France  Rémi
Carminati<EMAIL_ADDRESS>ESPCI Paris, PSL University, CNRS,
Institut Langevin, F-75005 Paris, France  Luis S. Froufe-Pérez  Physics
Department, University of Fribourg, CH-1700 Fribourg, Switzerland  Frank
Scheffold<EMAIL_ADDRESS>Physics Department, University of
Fribourg, CH-1700 Fribourg, Switzerland  Riccardo Sapienza  Imperial College
London, Blackett Laboratory, London, SW7 2AZ, United Kingdom  Silvia Vignolini
Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge,
CB2 1EW, United Kingdom  Juan José Sáenz  Donostia International Physics
Center (DIPC), 20018 Donostia-San Sebastian, Spain  IKERBASQUE, Basque
Foundation for Science, 48013 Bilbao, Spain

###### Abstract

The optics of correlated disordered media is a fascinating research topic
emerging at the interface between the physics of waves in complex media and
nanophotonics. Inspired by photonic structures in nature and enabled by
advances in nanofabrication processes, recent investigations have unveiled how
the design of structural correlations down to the subwavelength scale could be
exploited to control the scattering, transport and localization of light in
matter. From optical transparency to superdiffusive light transport to
photonic gaps, the optics of correlated disordered media challenges our
physical intuition and offers new perspectives for applications. This article
reviews the theoretical foundations, state-of-the-art experimental techniques
and major achievements in the study of light interaction with correlated
disorder, covering a wide range of systems – from short-range correlated
photonic liquids, to Lévy glasses containing fractal heterogeneities, to
hyperuniform disordered photonic materials. The mechanisms underlying light
scattering and transport phenomena are elucidated on the basis of rigorous
theoretical arguments. We overview the exciting ongoing research on mesoscopic
phenomena, such as transport phase transitions and speckle statistics, and the
current development of disorder engineering for applications such as light-
energy management and visual appearance design. Special efforts are finally
made to identify the main theoretical and experimental challenges to address
in the near future.

###### Contents

  1. I Introduction
  2. II Theory of multiple light scattering by correlated disordered media
    1. II.1 Theoretical framework
      1. II.1.1 Average field and self-energy
      2. II.1.2 Refractive index and extinction mean free path
      3. II.1.3 Multiple-scattering expansion
      4. II.1.4 Average intensity and four-point irreducible vertex
      5. II.1.5 Radiative transfer limit and scattering mean free path
      6. II.1.6 Transport mean free path and diffusion approximation
    2. II.2 Media with fluctuating continuous permittivity
      1. II.2.1 First insights on structural correlations
      2. II.2.2 Lorentz local fields
      3. II.2.3 Average exciting field
      4. II.2.4 Long-wavelength solutions for strongly fluctuating media: Bruggeman versus Maxwell-Garnett
      5. II.2.5 Expressions for the scattering and transport mean free paths from the average intensity
    3. II.3 Particulate media
      1. II.3.1 Expansion for identical scatterers
      2. II.3.2 Extinction mean free path and effective medium theories
      3. II.3.3 Scattering and transport mean free paths for resonant scatterers
    4. II.4 Summary and further remarks
  3. III Structural properties of correlated disordered media
    1. III.1 Continuous permittivity versus particulate models in practice
    2. III.2 Fluctuation-correlation relation
    3. III.3 Classes of correlated disordered media
      1. III.3.1 Short-range correlated disordered structures
      2. III.3.2 Polycrystalline structures
      3. III.3.3 Imperfect ordered structures
      4. III.3.4 Disordered hyperuniform structures
      5. III.3.5 Disordered fractal structures
    4. III.4 Numerical construction of correlated disordered media
    5. III.5 Fabrication of correlated disordered media
      1. III.5.1 Jammed colloidal packing
      2. III.5.2 Thermodynamically-driven self-assembly
      3. III.5.3 Optical and e-beam lithography
    6. III.6 Measuring structural correlations
  4. IV Modified transport parameters
    1. IV.1 Light scattering and transport in colloids and photonic materials
      1. IV.1.1 Impact of short-range correlations: first insights
      2. IV.1.2 Enhanced optical transparency
      3. IV.1.3 Tunable light transport in photonic liquids
      4. IV.1.4 Resonant effects in photonic glasses
      5. IV.1.5 Modified diffusion in imperfect photonic crystals
    2. IV.2 Anomalous transport in media with large-scale heterogeneity
      1. IV.2.1 Radiative transfer with non-exponential extinction
      2. IV.2.2 From normal to super-diffusion
  5. V Mesoscopic and near-field effects
    1. V.1 Photonic gaps in disordered media
      1. V.1.1 Definition and identification of photonic gaps in disordered media
      2. V.1.2 Competing viewpoints on the origin of photonic gaps
      3. V.1.3 Reports of photonic gaps in the litterature
    2. V.2 Mesoscopic transport and light localization
    3. V.3 Near-field speckles on correlated materials
      1. V.3.1 Intensity and field correlations in bulk speckle patterns
      2. V.3.2 Near-field speckles on dielectrics
    4. V.4 Local density of states fluctuations
  6. VI Photonics applications
    1. VI.1 Light trapping in layered media
    2. VI.2 Random lasing
    3. VI.3 Visual appearance
      1. VI.3.1 Photonic structures in nature
      2. VI.3.2 Synthetic structural colors
  7. VII Summary and perspectives
    1. VII.1 Near-field-mediated mesoscopic transport in 3D high-index correlated media
    2. VII.2 Mesoscopic optics in fractal and long-range correlated media
    3. VII.3 Towards novel applications
  8. A Green functions in Fourier space
    1. A.1 Dyadic Green tensor
    2. A.2 Dressed Green tensor
  9. B Derivation of Eq. (26)
  10. C Configurational average for statistically homogeneous systems
    1. C.1 Particle correlation functions
    2. C.2 Fluctuations of the number of particles in a volume
  11. D Local density of states and quasinormal modes

##  I Introduction

Correlated disordered media are non-crystalline heterogeneous materials
exhibiting pronounced spatial correlations in their structure and morphology.
The topic has bloomed in the context of optics and photonics, gradually
unveiling the considerable impact of structural correlations on the
scattering, transport, and localization of light in matter. These findings let
us envision novel types of materials with unprecedented optical
functionalities and raise a number of challenges in theoretical modelling,
material fabrication and optical spectroscopy. This article aims to provide an
overview of this emerging research field, starting from the basic principles
of light interaction with heterogeneous media to the most recent and still
actively debated topics.

The scattering of light by heterogeneous media has a long and venerable
history, which started more than a century ago with pioneering studies on the
refractive index of fluids of atoms or molecules Lorenz (1880); Lorentz
(1880), the electromagnetic scattering by particles Rayleigh (1899); Maxwell
Garnett (1904); Mie (1908), and the phenomenon of critical opalescence in
binary fluid mixtures Smoluchowski (1908); Einstein (1910); Ornstein and
Zernike (1914). The foundations of a rigorous theoretical treatment of
multiple light scattering were built in the 1930s Kirkwood (1936); Yvon (1937)
to take its full dimension a few decades later with various important
contributions Foldy (1945); Lax (1951, 1952); Twersky (1964); Keller (1964).
These early works already pointed out the key role played by structural
correlations on light scattering, a nice illustration of this being the
transparency of the cornea resulting from short-range correlations in
ensembles of discrete scatterers Maurice (1957); Hart and Farrell (1969);
Benedek (1971); Twersky (1975).

A new branch of research exploiting light waves to study mesoscopic phenomena
in disordered systems emerged in the 1980s, prompted by experimental
demonstrations of weak localization Tsang and Ishimaru (1984); Van Albada and
Lagendijk (1985); Wolf and Maret (1985) and theoretical predictions for the
three-dimensional Anderson localization of light John (1984); Anderson (1985).
The advent of photonic crystals Yablonovitch (1987); John (1987), wherein
photonic band gaps are created by a periodic modulation of the refractive
index in two or three dimensions, gave additional momentum to research by
greatly stimulating the development of nanofabrication techniques for high-
index dielectrics López (2003). The following decade witnessed a flourishment
of studies on periodic dielectric nanostructures Joannopoulos  _et al._ (2011)
and disordered media made of resonant (Mie) scatterers Fraden and Maret
(1990); Saulnier  _et al._ (1990); van Albada  _et al._ (1991); Busch and
Soukoulis (1995); Lagendijk and Van Tiggelen (1996) from two overlapping
communities Soukoulis (2012).

The importance of short-range structural correlations on light transport in
disordered systems Fraden and Maret (1990); Saulnier  _et al._ (1990) and of
random imperfections on light propagation in periodic systems Sigalas  _et
al._ (1996); Asatryan  _et al._ (1999); Vlasov  _et al._ (2000) was recognized
quite early. Research on correlated disordered media in optics however really
took off in the mid-2000s with experimental studies showing that disorder
could be engineered to harness light transport Rojas-Ochoa  _et al._ (2004);
García  _et al._ (2007); Barthelemy  _et al._ (2008). The surprising
observation of photonic gaps in disordered structures with short-range
correlations Edagawa  _et al._ (2008); Florescu  _et al._ (2009); Liew  _et
al._ (2011), reports of mesoscopic phenomena in imperfect photonic crystals
Conti and Fratalocchi (2008); Toninelli  _et al._ (2008); Garcia  _et al._
(2012) and the prospects of new generations of photonic devices like random
lasers Gottardo  _et al._ (2008), thin-film solar cells Vynck  _et al._
(2012); Oskooi  _et al._ (2012); Martins  _et al._ (2013) and integrated
spectrometers Redding  _et al._ (2013), greatly contributed to the emergence
of the research field. Figure 1 illustrates some of the key achievements and
applications of correlated disordered media in optics and photonics.

Figure 1:  Early achievements and applications of correlated disordered media
in optics and photonics. (a) Male mandrill (Mandrillus sphinx) blue facial
skin and cross-section of its dermis in the structurally colored area (blue)
that reveals parallel collagen fibres organised in a correlated array. Adapted
with permission from Prum and Torres (2004). (b) Modified light transport
(described by the inverse transport mean free path $\ell_{\text{t}}$) by
engineering short-range structural correlations thanks to Coulomb repulsion
between charged particles (filled symbols and solid line). Hard sphere systems
(open symbols and dotted line) exhibit weaker correlations. The dashed line is
a model neglecting completely structural correlations. Adapted with permission
from Rojas-Ochoa  _et al._ (2004). (c) Anomalous light transport in Lévy
glasses. A fractal heterogeneity is engineered by adding transparent spheres
with sizes varying over orders of magnitude in a host matrix. Transport can be
modelled by a truncated Lévy walk. On finite size samples, this leads to an
anomalous scaling of the total transmittance $T\sim L^{\alpha/2}$ (experiments
shown as symbols, lines are fits). Adapted with permission from Barthelemy
_et al._ (2008). (d) Light localization in randomly perturbed inverse opal
photonic crystals (upper left inset). Spatially-localized modes are observed
near the photonic band edge (lower left inset). Their typical spatial extent
(the localization length $\xi$) depends strongly on the degree of disorder.
Adapted with permission from Conti and Fratalocchi (2008). (e) Existence of
photonic gaps in amorphous photonic materials. Simulations of the spectral
density - a quantity proportional to the density of states - in a connected
amorphous diamond structure exhibiting short-range order shows a photonic gap
near $d/\lambda\simeq 0.23$, where $d$ is the average bond length. Adapted
with permission from Edagawa  _et al._ (2008). (f) Random lasing in two-
dimensional photonic structures with correlated disorder. Short-range
correlations are shown to increase the lasing efficiency at certain
frequencies due to enhanced optical confinement. Adapted with permission from
Noh  _et al._ (2011).

Important efforts have been made in recent years to elucidate the role of
structural correlations on the emergence of photonic gaps and on Anderson
localization in two-dimensional Florescu  _et al._ (2009); Conley  _et al._
(2014); Froufe-Pérez  _et al._ (2016, 2017) and three-dimensional disordered
systems Ricouvier  _et al._ (2019); Klatt  _et al._ (2019); Aubry  _et al._
(2020); Haberko  _et al._ (2020). Near-field interaction and light
polarization considerably complexify theoretical modelling Escalante and
Skipetrov (2017); Vynck  _et al._ (2016), explaining the widespread use of
full-wave numerical methods to address this issue, alongside phenomenological
models Naraghi  _et al._ (2015). The so-called hyperuniform disordered
structures Torquato and Stillinger (2003) have received considerable attention
in this context, leading to a wider exploration of their optical properties
Leseur  _et al._ (2016); Bigourdan  _et al._ (2019); Sheremet  _et al._ (2020)
and advances on top-down and bottom-up fabrication techniques Weijs  _et al._
(2015); Muller  _et al._ (2017); Ricouvier  _et al._ (2017); Maimouni  _et
al._ (2020).

In a different context, the interplay of order and disorder appeared quite
early as an essential ingredient to explain the colored appearance of certain
plants and animals Kinoshita and Yoshioka (2005). Research on natural photonic
structures continued at a fast pace with important findings, such as the
ubiquity of correlated disorder in animals exhibiting vivid diffuse blue
colors Noh  _et al._ (2010); Yin  _et al._ (2012); Magkiriadou  _et al._
(2012); Johansen  _et al._ (2017); Moyroud  _et al._ (2017), the use of short-
range correlations to reduce light reflectance Deparis  _et al._ (2009);
Siddique  _et al._ (2015) or structural anisotropy to enhance whiteness
Burresi  _et al._ (2014). Efforts are nowadays made to realize artificial
materials exhibiting correlated disorder to create materials with versatile
visual appearances Forster  _et al._ (2010); Takeoka (2012); Park  _et al._
(2014); Shang  _et al._ (2018); Schertel  _et al._ (2019a); Chan  _et al._
(2019); Goerlitzer  _et al._ (2018); Jacucci  _et al._ (2020a).

In this review, we will introduce the key concepts and techniques in the study
of light in correlated disordered media, assess the current state of knowledge
on the topic, and define the main challenges that lie ahead of us. Compared to
existing reviews on correlated disorder and disorder engineering in optics and
photonics Wiersma (2013); Shi  _et al._ (2013a); Yu  _et al._ (2020); Wang and
Zhao (2020), we provide here a broader view on the field and sufficient
technical details for the readers who wish dive into it, be it from the
theoretical or experimental side. This article is also an attempt to bridge
the gap between different research fields for which excellent textbooks
already exist, namely on random heterogeneous materials Torquato (2013),
multiple light scattering in complex media Tsang and Kong (2001); Sheng
(2006); Akkermans and Montambaux (2007); Carminati and Schotland (2021) and
periodic photonic crystals Joannopoulos  _et al._ (2011), and which may serve
as complementary literature. We focus on two and three-dimensional dielectric
materials, intentionally leaving aside one-dimensional dielectric structures
(i.e., layered media) Izrailev  _et al._ (2012) and metallic nanostructures
Shalaev (2002). We also do not discuss the very fertile fields of
metamaterials and metasurfaces, which show some apparent similarities with the
present topic in terms of theoretical models and concepts Mackay and Lakhtakia
(2015), yet with different scopes of application.

The remainder of the article is structured as follows. Section II introduces
the basic concepts and important quantities for light scattering and transport
in correlated disordered media, namely the extinction, scattering and
transport mean free paths. We derive mathematically explicit results, as a
function of the degree of structural correlations, from rigorous multiple
scattering theories for both continuous permittivity media and discrete
particulate media, emphasizing conceptual similarities between these two
viewpoints. Section III addresses the statistical description of the
structural properties of correlated disordered media. Different classes of
correlated systems are discussed, together with numerical and experimental
techniques to realize and characterize them. Section IV reviews experimental
and theoretical studies wherein structural correlations yield dramatic
variations of light transport parameters, including enhanced scattering in
colloidal suspensions of particles, optical transparency in hyperuniform
systems and anomalous diffusion in media with large-scale (fractal)
heterogeneities. Section V is concerned with emergent mesoscopic phenomena
relying on an interplay of order and disorder, most of which are not yet fully
understood. This includes the formation of photonic gaps and localized states
in disordered systems, and the statistical properties of near-field speckles
and local density of states. Section VI describes various applications of
correlated disordered media in optics and photonics, namely light management
in thin films, random lasing and visual appearance design. Section VII
concludes the review with a discussion on the open challenges in the field.

##  II Theory of multiple light scattering by correlated disordered media

The theoretical study of light propagation in disordered media is a
notoriously difficult problem that has experienced many developments for more
than a century. In this section, we introduce the basic concepts of multiple
light scattering by heterogeneities with the aim to give a solid theoretical
ground to the role of structural correlations on light scattering and
transport.

In Sec. II.1, we first focus on the “constitutive” linear relation between the
average electric field and the average polarization density in disordered
media, which allows us to introduce the notions of effective permittivity
tensor and extinction mean free path. Many of the derived results have been
used in the study of the effective optical response and homogeneization
processes of periodic and amorphous dielectrics Van Kranendonk and Sipe
(1977); Mackay and Lakhtakia (2015). We then derive the main equations that
govern the propagation of the average intensity and introduce the scattering
and transport mean free paths, two experimentally measurable quantities that
form the backbone of radiative transfer theory Chandrasekhar (1960).

Within this unique theoretical framework, we then address the light scattering
problem for systems described by either a continuous permittivity that
fluctuates in space [Sec. II.2] or discrete particles correlated in their
position [Sec. II.3], see Fig. 2. We derive analytical expressions for the
characteristic lengths, allowing us to show, on rigorous grounds, how
structural correlations impact scattering and transport.

We summarize the main outcomes of the theoretical analysis in Sec. II.4, for
the readers who wish to skip mathematical details.

Figure 2:  Disordered media may be described by a continuous permittivity
model (left) in the most general case, or by a particulate model (right) in
the case where the permittivity variation is compact.

###  II.1 Theoretical framework

####  II.1.1 Average field and self-energy

We consider a region of space filled with a non-magnetic, isotropic material
(relative permeability, $\mu(\mathbf{r})=1$), described by a scalar spatially-
varying relative permittivity $\epsilon({\bf r})$ in a uniform host medium
with relative permittivity $\epsilon_{\text{h}}$. Throughout this review, we
consider harmonic fields at frequency $\omega$ with the $e^{-{\rm i}\omega t}$
convention and drop the explicit dependence in the permittivities, fields,
etc. for clarity. In absence of charge and currents, the electric field ${\bf
E}({\bf r})$ at frequency $\omega$ satisfies a vector wave propagation
equation

$\displaystyle\bm{\mathbf{\nabla}}\times\bm{\mathbf{\nabla}}\times{\bf E}({\bf
r})-k_{0}^{2}\epsilon_{\text{h}}{\bf E}({\bf r})=k_{0}^{2}{\bf P}({\bf
r})/\epsilon_{0},$ (1)

where $k_{0}=\omega/c$ is the vacuum wave number and ${\bf P}({\bf
r})=\epsilon_{0}(\epsilon({\bf r})-\epsilon_{\text{h}}){\bf E}({\bf r})$ is
the polarization density (electric dipole moment per unit volume). The
permittivity variation $\Delta\epsilon({\bf r})=\epsilon({\bf
r})-\epsilon_{\text{h}}$ readily appears as the source of light scattering in
the system.

Upon statistical average over disorder realizations, denoted here as
$\langle\dots\rangle$, Eq. (1) describes the propagation of the average field
$\langle{\bf E}({\bf r})\rangle$ with the average polarization density
$\langle{\bf P}({\bf r})\rangle$ as a source term. The difficulty to solve
this general problem essentially comes from the fact that the fluctuations of
the electric field and of the permittivity are not statistically-independent,
i.e., $\langle\Delta\epsilon({\bf r}){\bf E}({\bf
r})\rangle\neq\langle\Delta\epsilon({\bf r})\rangle\langle{\bf E}({\bf
r})\rangle$. The standard approach to solve this problem is to make a
preliminary ansatz about the effective permittivity of the medium, thereby
establishing a constitutive linear relation between the average field and the
average polarization density, and then perturbatively calculate correlations
in orders the statistical cumulants of the permittivity spatial fluctuations
around it.

Let us then rewrite Eq. (1) as Ryzhov  _et al._ (1965)

$\displaystyle\bm{\mathbf{\nabla}}\times\bm{\mathbf{\nabla}}\times{\bf E}({\bf
r})-k_{0}^{2}\epsilon_{\text{b}}{\bf E}({\bf
r})={\bm{\mathbf{\mathcal{X}}}}({\bf r}){\bf E}({\bf r}),$ (2)

where $\epsilon_{\text{b}}$ is a constant, auxiliary, background permittivity
that can differ from the permittivity of the host medium
$\epsilon_{\text{h}}$, and ${\bm{\mathbf{\mathcal{X}}}}$ is an effective
scattering potential (or electric susceptibility) defined as

$\displaystyle{\bm{\mathbf{\mathcal{X}}}}({\bf r})=k_{0}^{2}(\epsilon({\bf
r})-\epsilon_{\text{b}})\bm{\mathbf{1}},$ (3)

with $\bm{\mathbf{1}}$ the unit tensor. The average field $\langle{\bf
E}\rangle$ then fulfills the vector wave equation

$\displaystyle\bm{\mathbf{\nabla}}\times\bm{\mathbf{\nabla}}\times\langle{\bf
E}\rangle({\bf r})-k_{\text{b}}^{2}\langle{\bf E}\rangle({\bf
r})=\langle{\bm{\mathbf{\mathcal{X}}}}({\bf r}){\bf E}({\bf r})\rangle,$ (4)

where $k_{\text{b}}^{2}=k_{0}^{2}\epsilon_{\text{b}}$. We now introduce the
susceptibility tensor ${\bm{\mathbf{\Sigma}}}$, usually called “self-energy”
or “mass operator”, following the language of many-body scattering theory
Dyson (1949a), as

$\displaystyle\langle{\bm{\mathbf{\mathcal{X}}}}({\bf r}){\bf E}({\bf
r})\rangle\equiv\int\bm{\mathbf{\Sigma}}({\bf r},{\bf r}^{\prime})\langle{\bf
E}\rangle({\bf r}^{\prime})d{\bf r}^{\prime},$ (5)

with

$\displaystyle\bm{\mathbf{\Sigma}}({\bf r},{\bf
r}^{\prime})=k_{0}^{2}\left({\bm{\mathbf{\epsilon}}}_{\text{eff}}({\bf r},{\bf
r}^{\prime})-\epsilon_{\text{b}}\bm{\mathbf{1}}\delta({\bf r}-{\bf
r}^{\prime})\right).$ (6)

The self-energy depends on the non-local effective permittivity tensor that
results from multiple scattering in the disordered medium. Hereafter, we
assume that the system has a proper thermodynamic limit in which it becomes
spatially homogeneous and translationally invariant on average, i.e.
${\bm{\mathbf{\epsilon}}}_{\text{eff}}({\bf r},{\bf
r}^{\prime})={\bm{\mathbf{\epsilon}}}_{\text{eff}}({\bf r}-{\bf r}^{\prime})$.
Aspects related to the effective medium description in large but finite
systems and the deep connection with the Ewald-Oseen extinction theorem Hynne
and Bullough (1987); Van Kranendonk and Sipe (1977) will thus not be discussed
here. In Fourier space, the self-energy is given by

$\displaystyle\bm{\mathbf{\Sigma}}({\bf k},{\bf k}^{\prime})$ $\displaystyle=$
$\displaystyle\iint e^{-{\rm i}{\bf k}\cdot{\bf r}}\bm{\mathbf{\Sigma}}({\bf
r}-{\bf r}^{\prime})e^{{\rm i}{\bf k}^{\prime}\cdot{\bf r}^{\prime}}d{\bf
r}d{\bf r}^{\prime}$ (7) $\displaystyle=$
$\displaystyle(2\pi)^{3}{\bm{\mathbf{\Sigma}}}({\bf k})\delta({\bf k}-{\bf
k}^{\prime}).$ (8)

It is further convenient to decompose ${\bm{\mathbf{\Sigma}}}({\bf k})$ into
its transverse ($\perp$) and longitudinal ($\parallel$) components as

$\displaystyle{\bm{\mathbf{\Sigma}}}({\bf k})={\Sigma}_{\perp}({\bf
k})\left(\bm{\mathbf{1}}-{\bf u}\otimes{\bf
u}\right)+{\Sigma}_{\parallel}({\bf k}){\bf u}\otimes{\bf u},$ (9)

where ${\bf u}={\bf k}/|{\bf k}|$ and ${\bf u}\otimes{\bf u}$ is the outer
tensor product between ${\bf u}$ and itself. Defining ${\bf e}$ as the unit
polarization vector with ${\bf e}\cdot{\bf u}=0$, the transverse component of
the self-energy reads

$\displaystyle{\Sigma}_{\perp}({\bf k})={\bf
e}\cdot{\bm{\mathbf{\Sigma}}}({\bf k}){\bf e}.$ (10)

####  II.1.2 Refractive index and extinction mean free path

To understand the role played by the self-energy in wave propagation and
scattering, we can seek for transversal solutions of the vector wave
propagation equation of the form

$\displaystyle\langle{\bf E}\rangle({\bf r})$ $\displaystyle=$ $\displaystyle
E_{0}{\bf e}e^{{\rm i}{\bf k}_{\text{eff}}\cdot{\bf r}},$ (11)

with ${\bf k}_{\text{eff}}=k_{0}n_{\text{eff}}{\bf u}$ the wavevector
describing propagation in a homogeneous medium with effective refractive index
$n_{\text{eff}}$. Assuming a statistically isotropic system, substituting Eq.
(11) in Eq. (4) with Eq. (5) and following the decomposition in Eq. (9) with
Eq. (10) leads to a transcendental equation for the effective wave number

$\displaystyle k_{\text{eff}}=k_{0}n_{\text{eff}}$ $\displaystyle=$
$\displaystyle\sqrt{k_{\text{b}}^{2}+{\Sigma}_{\perp}(k_{\text{eff}})},$ (12)
$\displaystyle\equiv$ $\displaystyle k_{\text{r}}+{\rm
i}\frac{1}{2\ell_{\text{e}}}.$

The real part of the effective index
$\text{Re}[n_{\text{eff}}]=k_{\text{r}}/k_{0}$ describes the phase velocity of
the average field (often called “coherent” or “ballistic” component) in the
material, while the imaginary part
$\text{Im}[n_{\text{eff}}]=(2k_{0}\ell_{\text{e}})^{-1}$ describes its
exponential decay with propagation due to absorption and/or scattering on a
characteristic length scale that is the extinction mean free path,
$\ell_{\text{e}}$. In the weak extinction regime, i.e.
$\text{Im}{\Sigma}_{\perp}(k_{\text{eff}})\ll
k_{\text{b}}^{2}+\text{Re}{\Sigma}_{\perp}(k_{\text{eff}})$ and
$k_{\text{r}}\ell_{\text{e}}\gg 1$, Eq. (12) leads to

$\displaystyle\frac{1}{\ell_{\text{e}}}\simeq\frac{\text{Im}\Sigma_{\perp}(k_{\text{r}})}{k_{\text{r}}}.$
(13)

In non-absorbing dielectric materials, extinction is purely driven by
scattering ($\ell_{\text{e}}=\ell_{\text{s}}$ with $\ell_{\text{s}}$ the
scattering mean free path). The problem of light scattering by correlated
disordered media can therefore be apprehended by determining the self-energy
of the system.

####  II.1.3 Multiple-scattering expansion

A key ingredient in solving multiple light scattering problems is the
electromagnetic dyadic Green tensor ${\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime})$, that is the solution of Eq. (2) describing the electric field
produced at a point ${\bf r}$ by a radiating point electric dipole at ${\bf
r}^{\prime}$ in a homogenous medium with permittivity $\epsilon_{\text{b}}$,
and is given by

${\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime})=-\frac{\bm{\mathbf{1}}}{3k_{\text{b}}^{2}}\delta({\bf r}-{\bf
r}^{\prime})\\\ +\lim_{a\rightarrow 0}\Theta(|{\bf r}-{\bf
r}^{\prime}|-a)\left\\{\left(\bm{\mathbf{1}}+\frac{\bm{\mathbf{\nabla}}\otimes\bm{\mathbf{\nabla}}}{k_{\text{b}}^{2}}\right)\frac{e^{{\rm
i}k_{\text{b}}|{\bf r}-{\bf r}^{\prime}|}}{4\pi|{\bf r}-{\bf
r}^{\prime}|}\right\\},$ (14)

where $\Theta$ is the Heaviside step function. The Dirac delta function in the
right hand side of Eq. (14) gives the well-known singularity in the source
region while the second term corresponds to the non-singular, principal value
of the Green function, $\text{PV}\left\\{{\bf G}_{0}\right\\}$ Yaghjian
(1980); Van Bladel and Van Bladel (1991). The Green function enables writing a
general solution of the wave equation as a Lippmann-Schwinger equation

$\displaystyle{\bf E}({\bf r})={\bf E}_{\text{b}}({\bf r})+\int{{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bm{\mathbf{\mathcal{X}}}}({\bf
r}^{\prime}){\bf E}({\bf r}^{\prime})d{\bf r}^{\prime}},$ (15)

where ${\bf E}_{\text{b}}({\bf r})$ is the solution of the homogeneous
problem, which can be seen as a background (incident) field with wave number
$k_{\text{b}}$. Equation (15) can conveniently be written in operator form as

$\displaystyle{\bf E}={\bf
E}_{\text{b}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}{\bf
E},$ (16)

where ${\bm{\mathbf{\mathcal{G}}}}_{\text{b}}$ is an integral operator.
Equation (16) can be formally solved by successive iterations, leading to a
multiple-scattering expansion on orders of ${\bm{\mathbf{\mathcal{X}}}}$
(i.e., single scattering, double scattering, etc.). Eventually, all multiple-
scattering orders are taken into account by defining a transition operator
${\bm{\mathbf{\mathcal{T}}}}$ relating the polarization induced in the medium
to the background field, as Messiah (1999)

$\displaystyle{\bf E}={\bf
E}_{\text{b}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{T}}}}{\bf
E}_{\text{b}},$ (17)

with

$\displaystyle{\bm{\mathbf{\mathcal{T}}}}$ $\displaystyle=$
$\displaystyle{\bm{\mathbf{\mathcal{X}}}}+{\bm{\mathbf{\mathcal{X}}}}{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}+\cdots$
(18) $\displaystyle=$
$\displaystyle\left[\bm{\mathbf{1}}-{\bm{\mathbf{\mathcal{X}}}}{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\right]^{-1}{\bm{\mathbf{\mathcal{X}}}}.$

Keeping only the lowest order of the expansion,
${\bm{\mathbf{\mathcal{T}}}}={\bm{\mathbf{\mathcal{X}}}}$, is known as the
Born approximation, wherein one neglects the contribution of the field
scattered by neighboring permittivity fluctuations, resulting in a local,
single-scattering process. In the general case, multiple scattering occurs
leading to a transition operator ${\bm{\mathbf{\mathcal{T}}}}$ that is
spatially non-local, ${\bm{\mathbf{\mathcal{T}}}}{\bf
E}_{\text{b}}\equiv\int{{\bf T}({\bf r},{\bf r}^{\prime}){\bf
E}_{\text{b}}({\bf r}^{\prime})d{\bf r}^{\prime}}$.

Upon statistical average of Eqs. (16) and (17), and having Eq. (5), we finally
reach a general expression for the average field as a function of the self-
energy operator $\bm{\mathbf{\mathit{\Sigma}}}$, known as the Dyson equation
Yvon (1937); Dyson (1949a, b); Rytov  _et al._ (1989)

$\displaystyle\langle{\bf E}\rangle$ $\displaystyle=$ $\displaystyle{\bf
E}_{\text{b}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\bm{\mathbf{\mathit{\Sigma}}}\langle{\bf
E}\rangle,$ (19)

with

$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}=\langle{\bm{\mathbf{\mathcal{T}}}}\rangle\left[\bm{\mathbf{1}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\langle{\bm{\mathbf{\mathcal{T}}}}\rangle\right]^{-1}.$
(20)

In summary, the disordered medium is described as a permittivity that
fluctuates around an auxiliary background permittivity $\epsilon_{\text{b}}$
via the scattering potential ${\bm{\mathbf{\mathcal{X}}}}$ [Eq. (3)]. The
field propagates from fluctuation to fluctuation via the dyadic Green tensor
${\bm{\mathbf{\mathcal{G}}}}_{\text{b}}$ in the homogeneous background with
wave number $k_{\text{b}}$ [Eq. (14)]. The multiple scattering process on the
scattering potential ${\bm{\mathbf{\mathcal{X}}}}$ is described (to infinite
order) via the transition operator ${\bm{\mathbf{\mathcal{T}}}}$ of the entire
system [Eq. (18)]. The average transition operator
$\langle{\bm{\mathbf{\mathcal{T}}}}\rangle$ finally defines the self-energy
$\bm{\mathbf{\mathit{\Sigma}}}$ [Eq. (20)], which describes the propagation of
the average field $\langle{\bf E}\rangle$ in the disordered medium, and
therefore the extinction (i.e., scattering and absorption) efficiency of the
system [Eq. (12)].

####  II.1.4 Average intensity and four-point irreducible vertex

The theory above describes how an incident wave gets attenuated by scattering
and absorption in a heterogeneous medium, but not how the light intensity (or
energy density) propagates in the medium. Indeed, the energy that is removed
by scattering from the average field intensity $|\left\langle{\bf E}({\bf
r})\right\rangle|^{2}$, also known as coherent intensity, is redistributed
along different directions to form, upon sequences of scattering events, a
spatial distribution for the average intensity $\left\langle|{\bf E}({\bf
r})|^{2}\right\rangle$ that can strongly differ from a simple exponentially-
decaying profile. The formalism used for the average field will now be applied
to the average intensity. This will lead to the definition of the scattering
and transport mean free paths, two additional typical lengths at the root of
light propagation in disordered media Apresyan and Kravtsov (1996); van Rossum
and Nieuwenhuizen (1999); Rytov  _et al._ (1989).

Let us then consider the spatial correlation function of the electric field,
or “coherence matrix” Mandel and Wolf (1995), $\bm{\mathbf{C}}({\bf r},{\bf
r}^{\prime})\equiv\left\langle{\bf E}({\bf r})\otimes{\bf E}^{*}({\bf
r}^{\prime})\right\rangle$, with ∗ denoting complex conjugate. Starting from
the Lippmann-Schwinger equation [Eq. (16)], one can easily show that
$\bm{\mathbf{C}}$ depends on the spatial correlation of the polarization
density in the effective scattering potential
$\left\langle({\bm{\mathbf{\mathcal{X}}}}({\bf r}){\bf E}({\bf
r}))\otimes({\bm{\mathbf{\mathcal{X}}}}^{*}({\bf r}^{\prime}){\bf E}^{*}({\bf
r}^{\prime}))\right\rangle$. Similarly to the self-energy that allowed us to
relate the average polarization to the average field [Eq. (5)] eventually
leading to the Dyson equation [Eq. (19)], we can introduce here an operator
$\bm{\mathbf{\Gamma}}$ known as the “four-point irreducible vertex” – or
intensity vertex – that relates the effective polarization density correlation
to the electric field correlation. This eventually leads to a closed-form and
exact equation, known as the Bethe-Salpeter equation Salpeter and Bethe
(1951), which reads

$\bm{\mathbf{C}}({\bf r},{\bf
r}^{\prime})=\left\langle\bm{\mathbf{E}}(\bm{\mathbf{r}})\right\rangle\otimes\left\langle\bm{\mathbf{E}}^{*}({\bf
r}^{\prime})\right\rangle+\int\left\langle\bm{\mathbf{G}}(\bm{\mathbf{r}},\bm{\mathbf{r}}_{1})\right\rangle\otimes\left\langle\bm{\mathbf{G}}^{*}({\bf
r}^{\prime},{\bf r}_{1}^{\prime})\right\rangle\\\
\cdot\bm{\mathbf{\Gamma}}({\bf r}_{1},{\bf r}_{2},{\bf r}_{1}^{\prime},{\bf
r}_{2}^{\prime})\cdot\bm{\mathbf{C}}({\bf r}_{2},{\bf
r}_{2}^{\prime})\mathrm{d}{\bf r}_{1}\mathrm{d}{\bf
r}_{1}^{\prime}\mathrm{d}{\bf r}_{2}\mathrm{d}{\bf r}_{2}^{\prime},$ (21)

where the average Green function $\left\langle\bm{\mathbf{G}}\right\rangle$ is
given by the Dyson equation [Eq. (19)],

$\displaystyle\left\langle\bm{\mathbf{\mathcal{G}}}\right\rangle={\bm{\mathbf{\mathcal{G}}}}_{\text{b}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\bm{\mathbf{\mathit{\Sigma}}}\left\langle\bm{\mathbf{\mathcal{G}}}\right\rangle.$
(22)

The symbol $\cdot$ denotes here a tensor contraction, defined such that
$(\bm{\mathbf{A}}\otimes\bm{\mathbf{B}})\cdot(\bm{\mathbf{e}}\otimes\bm{\mathbf{f}})=(\bm{\mathbf{A}}\bm{\mathbf{e}})\otimes(\bm{\mathbf{B}}\bm{\mathbf{f}})$
and
$(\bm{\mathbf{A}}\otimes\bm{\mathbf{B}})\cdot(\bm{\mathbf{C}}\otimes\bm{\mathbf{D}})=(\bm{\mathbf{A}}\bm{\mathbf{C}})\otimes(\bm{\mathbf{B}}\bm{\mathbf{D}})$,
where $\bm{\mathbf{e}},\bm{\mathbf{f}}$ are vectors and
$\bm{\mathbf{A}},\bm{\mathbf{B}},\bm{\mathbf{C}},\bm{\mathbf{D}}$ matrices.

The first term in Eq. (21) is the source term, given by the coherent intensity
that decays exponentially within the scattering medium. The second term
expresses the field correlation as a multiple scattering process, wherein the
propagation is described by the average Green tensors and scattering by the
vertex $\bm{\mathbf{\Gamma}}$ that connects two pairs of points (for the field
and the complex conjugate). Following similar steps as those leading to Eq.
(20), we obtain a completely general and exact expression for
$\bm{\mathbf{\mathit{\Gamma}}}$ (here, in operator notation) as Carminati and
Schotland (2021)

$\displaystyle\bm{\mathbf{\mathit{\Gamma}}}=\left[{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}^{*}\right]^{-1}$
$\displaystyle\big{[}$
$\displaystyle\bm{\mathbf{1}}+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle+{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}^{*}\left\langle{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle$
(23) $\displaystyle+$
$\displaystyle\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}^{*}\left\langle{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle\big{]}^{-1}.$

We proceed by establishing a link between $\bm{\mathbf{\mathit{\Gamma}}}$ and
the light scattering and transport properties of the medium. Towards this aim,
it is convenient to rewrite the Bethe-Salpeter equation in Fourier space. A
closed-form analytical formula can be obtained by assuming that the scattering
events take place on distances larger than the wavelength, such that the
average Green tensor can be approximated as its transverse component only,

$\displaystyle\left\langle\bm{\mathbf{G}}({\bf k})\right\rangle$
$\displaystyle=$ $\displaystyle\left[k^{2}\bm{\mathbf{P}}({\bf
u})-k_{\text{b}}^{2}\bm{\mathbf{1}}-\bm{\mathbf{\Sigma}}({\bf
k})\right]^{-1},$ (24) $\displaystyle\simeq$ $\displaystyle\left\langle
G_{\perp}({\bf k})\right\rangle\bm{\mathbf{P}}({\bf u}),$ (25)

where $\bm{\mathbf{P}}({\bf u})=\bm{\mathbf{1}}-{\bf u}\otimes{\bf u}$ is the
transverse projection operator and $\left\langle G_{\perp}({\bf
k})\right\rangle=\left[k^{2}-k_{\text{b}}^{2}-\Sigma_{\perp}({\bf
k})\right]^{-1}$ is the (scalar) transverse component. After some algebra
provided in details in Appendix B, we reach

$\left[\left({\bf k}-\frac{{\bf q}}{2}\right)^{2}-\left({\bf k}+\frac{{\bf
q}}{2}\right)^{2}-\Sigma_{\perp}^{*}\left({\bf k}-\frac{{\bf
q}}{2}\right)+\Sigma_{\perp}\left({\bf k}+\frac{{\bf
q}}{2}\right)\right]\bm{\mathbf{L}}_{\perp}({\bf k},{\bf q})\\\
=\left[\left\langle G_{\perp}\left({\bf k}+\frac{{\bf
q}}{2}\right)\right\rangle-\left\langle G_{\perp}^{*}\left({\bf k}-\frac{{\bf
q}}{2}\right)\right\rangle\right]\int\bar{\bm{\mathbf{\Gamma}}}_{\perp}\left({\bf
k}+\frac{{\bf q}}{2},{\bf k}^{\prime}+\frac{{\bf q}}{2},{\bf k}-\frac{{\bf
q}}{2},{\bf k}^{\prime}-\frac{{\bf
q}}{2}\right)\cdot\bm{\mathbf{L}}_{\perp}\left({\bf k}^{\prime},{\bf
q}\right)\frac{\mathrm{d}\bm{\mathbf{k}}^{\prime}}{(2\pi)^{3}},$ (26)

where we relied on the statistical homogeneity and translational invariance of
the medium and neglected the source term. The field correlation, described by
a new function

$\displaystyle\bm{\mathbf{L}}_{\perp}({\bf k},{\bf
q})\equiv\bm{\mathbf{C}}_{\perp}\left({\bf k}+\frac{{\bf q}}{2},{\bf
k}-\frac{{\bf q}}{2}\right),$ (27)

with

$\displaystyle\bm{\mathbf{C}}_{\perp}\left(\bm{\mathbf{k}},\bm{\mathbf{k}}^{\prime}\right)=\bm{\mathbf{P}}({\bf
u})\otimes\bm{\mathbf{P}}({\bf
u}^{\prime})\cdot\bm{\mathbf{C}}\left(\bm{\mathbf{k}},\bm{\mathbf{k}}^{\prime}\right),$
(28)

depends only on the transverse part of the intensity vertex, which is given by

$\bm{\mathbf{\Gamma}}_{\perp}(\bm{\mathbf{k}},\bm{\mathbf{\upkappa}},\bm{\mathbf{k}}^{\prime},\bm{\mathbf{\upkappa}}^{\prime})=\bm{\mathbf{P}}({\bf
u})\otimes\bm{\mathbf{P}}({\bf
u}^{\prime})\cdot\bm{\mathbf{\Gamma}}(\bm{\mathbf{k}},\bm{\mathbf{\upkappa}},\bm{\mathbf{k}}^{\prime},\bm{\mathbf{\upkappa}}^{\prime})\\\
=(2\pi)^{3}\delta({\bf k}-\bm{\mathbf{\upkappa}}-{\bf
k}^{\prime}-\bm{\mathbf{\upkappa}}^{\prime})\bar{\bm{\mathbf{\Gamma}}}_{\perp}(\bm{\mathbf{k}},\bm{\mathbf{\upkappa}},\bm{\mathbf{k}}^{\prime},\bm{\mathbf{\upkappa}}^{\prime}).$
(29)

Equation (26) is very general, as it considers all multiple scattering events
within the medium and does not make any explicit assumption on the kind of
disorder. Note however that neglecting the longitudinal component of the Green
tensor implicitly excludes situations wherein near-field contributions are
important (e.g., in dense packings of high-index resonant particles).

####  II.1.5 Radiative transfer limit and scattering mean free path

Further approximations are required to extract analytical expressions for the
scattering and transport properties of the medium. First, we take the large-
scale approximation $|{\bf q}|\ll\\{|{\bf k}|,|{\bf k}^{\prime}|\\}$, also
known as the radiative transfer limit Barabanenkov and Finkel’berg (1968);
Ryzhik  _et al._ (1996); Henkel (1997), which assumes that the average
intensity varies on length scales $2\pi/|{\bf q}|$ much larger than the
typical wavelength in the medium $2\pi/k_{\text{r}}$. Since the average field
and intensity vary on the scale of $\ell_{\text{e}}$, this approximation holds
in the weak extinction limit, $k_{\text{r}}\ell_{\text{e}}\gg 1$. Equation
(26) thus becomes

$\left[-2{\bf k}\cdot{\bf q}+2i\operatorname{Im}\Sigma_{\perp}({\bf
k})\right]\bm{\mathbf{L}}_{\perp}({\bf k},{\bf
q})=2i\operatorname{Im}\left\langle G_{\perp}({\bf k})\right\rangle\\\
\times\int\bar{\bm{\mathbf{\Gamma}}}_{\perp}\left({\bf k},{\bf
k}^{\prime},{\bf k},{\bf k}^{\prime}\right)\cdot\bm{\mathbf{L}}_{\perp}({\bf
k}^{\prime},{\bf q})\frac{\mathrm{d}{\bf k}^{\prime}}{(2\pi)^{3}}.$ (30)

In a dilute material, the scattering correction in the average Green function
is considered to be weak, i.e. $|\Sigma_{\perp}|\ll k_{\text{b}}^{2}$. Using
the relation

$\displaystyle\lim_{\epsilon\to
0^{+}}\frac{1}{x-x_{0}-i\epsilon}=\operatorname{PV}\left[\frac{1}{x-x_{0}}\right]+i\pi\delta(x-x_{0}),$
(31)

where $\operatorname{PV}$ stands for the Cauchy principal value operator, the
imaginary part of the average Green function reduces to

$\displaystyle\operatorname{Im}\left\langle G_{\perp}({\bf
k})\right\rangle=\pi\delta\left[k^{2}-k_{\text{b}}^{2}-\operatorname{Re}\Sigma_{\perp}({\bf
k})\right].$ (32)

This relation fixes the real part of the effective wavevector
$k_{\text{r}}=\operatorname{Re}k_{\text{eff}}$ to

$\displaystyle
k_{\text{r}}=\sqrt{k_{\text{b}}^{2}+\operatorname{Re}\Sigma_{\perp}(k_{\text{r}})},$
(33)

which is the so-called “on-shell approximation”. Second, we assume that the
field is fully depolarized, which is a perfectly acceptable hypothesis after
many scattering events (i.e., far from a coherent, polarized source) Bicout
and Brosseau (1992); Gorodnichev  _et al._ (2014); Vynck  _et al._ (2016).
This gives

$\displaystyle\bm{\mathbf{C}}\left(\bm{\mathbf{k}},\bm{\mathbf{k}}^{\prime}\right)=C\left(\bm{\mathbf{k}},\bm{\mathbf{k}}^{\prime}\right)\bm{\mathbf{1}}$
(34)

leading to

$\displaystyle\bm{\mathbf{L}}_{\perp}({\bf k},{\bf q})=L({\bf k},{\bf
q})\bm{\mathbf{P}}({\bf u})\otimes\bm{\mathbf{P}}({\bf
u}^{\prime})\cdot\bm{\mathbf{1}}.$ (35)

An inverse Fourier transform of the trace of Eq. (30) together with Eqs. (32)
and (35) eventually leads to the famous Radiative Transfer Equation (RTE)
Chandrasekhar (1960)

$\displaystyle\left[{\bf u}\cdot\bm{\mathbf{\nabla}}_{{\bf
r}}+\frac{1}{\ell_{\text{e}}}\right]I({\bf r},{\bf
u})=\frac{1}{\ell_{\text{s}}}\int p({\bf u},{\bf u}^{\prime})I\left({\bf
r},{\bf u}^{\prime}\right)\mathrm{d}{\bf u}^{\prime},$ (36)

$I$ is the specific intensity, which can be interpreted as a local (at
position $\bm{\mathbf{r}}$) and directional (on direction ${\bf u}$) radiative
flux, defined as

$\displaystyle\delta(k-k_{\text{r}})I({\bf r},{\bf u})=L({\bf r},{\bf k}).$
(37)

$\ell_{\text{s}}$ and $p({\bf u},{\bf u}^{\prime})$ are the scattering mean
free path and the phase function, describing respectively the average distance
between two scattering events and the angular diagram for an incident
planewave along ${\bf u}^{\prime}$ scattered along direction ${\bf u}$. Both
quantities are related to the intensity vertex via the relation

$\frac{1}{\ell_{\text{s}}}p({\bf u},{\bf
u}^{\prime})=\frac{1}{32\pi^{2}}\operatorname{Tr}\left[\bm{\mathbf{P}}({\bf
u})\otimes\bm{\mathbf{P}}({\bf u})\right.\\\
\left.\cdot\bar{\bm{\mathbf{\Gamma}}}(k_{\text{r}}{\bf u},k_{\text{r}}{\bf
u}^{\prime},k_{\text{r}}{\bf u},k_{\text{r}}{\bf
u}^{\prime})\cdot\bm{\mathbf{P}}({\bf u}^{\prime})\otimes\bm{\mathbf{P}}({\bf
u}^{\prime})\cdot\bm{\mathbf{1}}\right],$ (38)

and the phase function is normalized as

$\displaystyle\int p({\bf u},{\bf u}^{\prime})\mathrm{d}{\bf u}^{\prime}=1.$
(39)

The trace corresponds to an average over all possible polarization states for
the incident and scattered fields. The RTE [Eq. (36)] can be seen as an energy
balance Chandrasekhar (1960). The spatial variation of the specific intensity
(term involving the derivative) is due to the loss induced by extinction along
the direction ${\bf u}$ [term involving $\ell_{\text{e}}$] and the gain from
scattering from direction ${\bf u}^{\prime}$ to direction ${\bf u}$ [term
involving $\ell_{\text{s}}$ and $p({\bf u},{\bf u}^{\prime})$]. Equations (38)
and (39) show that $\ell_{\text{s}}$, the key quantity to describe the
scattering strength of a medium, is obtained in the radiative transfer limit
from the angular integral of the intensity vertex.

Previously, we showed that the extinction mean free path $\ell_{\text{e}}$
could be obtained from the self-energy $\bm{\mathbf{\Sigma}}$ and noted that,
in absence of absorption, we should have $\ell_{\text{e}}=\ell_{\text{s}}$,
the latter being defined from the intensity vertex $\bm{\mathbf{\Gamma}}$. It
is worth emphasizing at this point that the two operators are indeed formally
linked by the Ward identity, which may be seen as a generalization of the
extiction (optical) theorem and ensures energy conservation Lagendijk and Van
Tiggelen (1996); Apresyan and Kravtsov (1996); Tsang and Kong (2001); Sheng
(2006); Cherroret  _et al._ (2016).

####  II.1.6 Transport mean free path and diffusion approximation

Many experiments on light in disordered media are performed in situations
where light experiences not just a few but many scattering events on average.
It turns out that many observables accessible experimentally depend on yet
another length scale, that is the transport mean free path $\ell_{\text{t}}$,
which we will introduce here.

We start by taking the first moment of Eq. (36) (i.e., multiplying both sides
by ${\bf u}$ and integrating over ${\bf u}$) which directly leads to

$\displaystyle\int\left[{\bf u}\cdot\bm{\mathbf{\nabla}}_{{\bf r}}I({\bf
r},{\bf u})\right]{\bf u}\mathrm{d}{\bf u}+\frac{1}{\ell_{\text{t}}}{\bf
j}({\bf r})=0,$ (40)

where ${\bf j}({\bf r})=\int I({\bf r},{\bf u}){\bf u}\mathrm{d}{\bf u}$ is
the radiative flux vector and

$\displaystyle\ell_{\text{t}}\equiv\frac{\ell_{\text{s}}}{1-g},$ (41)

is the transport mean free path.

$\displaystyle g=\int p({\bf u},{\bf u}^{\prime}){\bf u}\cdot{\bf
u}^{\prime}\mathrm{d}{\bf u}$ (42)

is the average cosine of the phase function, or scattering anisotropy factor.
It can be interpreted as a renormalization factor of the scattering mean free
path due to predominantly forward ($0<g<1$) or backward ($-1<g<0$) single
scattering. Physically, $\ell_{\text{t}}$ describes the average distance
required for light to be fully randomized (e.g., forward scattering increases
this distance at constant scattering strength $\ell_{\text{s}}$). Structural
correlations impact the transport mean free path via both the scattering mean
free path $\ell_{\text{s}}$ and the scattering anisotropy described by $g$.

After propagation on distances much larger than $\ell_{\text{t}}$, the
specific intensity becomes quasi-isotropic. Expanding the specific intensity
in the RTE [Eq. (36)] on Legendre polynomials to first order in ${\bf u}$,
that is the so-called $P_{1}$-approximation Ishimaru (1978), leads to the
diffusion equation. In the steady-state regime, it reads

$\displaystyle-\mathcal{D}\Delta u({\bf r})=s({\bf r})$ (43)

where $u({\bf r})=v_{\text{E}}^{-1}\int I({\bf r},{\bf u})\mathrm{d}{\bf u}$
is the energy density with $v_{\text{E}}$ the energy velocity Lagendijk and
Van Tiggelen (1996), $\mathcal{D}=v_{\text{E}}\ell_{\text{t}}/3$ is the
diffusion constant and $s$ is a source term given by the ballistic intensity.
Resolving this equation in a slab geometry of thickness $L$ under planewave
illumination at normal incidence gives the following asymptotic behavior for
the total transmittance

$\displaystyle T\sim\frac{5}{3}\frac{\ell_{\text{t}}}{L},$ (44)

also known as Ohm’s law. Many transport observations in the diffusive limit
depend directly on the transport mean free path, including the linewidth of
the coherent backscattering cone Akkermans  _et al._ (1986, 1988), the time-
resolved transmittance and reflectance Contini  _et al._ (1997), and long-
range speckle correlations Shapiro (1999); Scheffold and Maret (1998); Fayard
_et al._ (2015).

In summary, the average transition operator of the medium
$\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle$ defines the four-point
irreducible vertex $\bm{\mathbf{\mathit{\Gamma}}}$ [Eq. (23)]. Neglecting
near-field interaction between scattering elements, taking the radiative
transfer limit, using the on-shell approximation and assuming fully
depolarized light allow us to relate the transverse component of
$\bm{\mathbf{\mathit{\Gamma}}}$ to the scattering mean free path
$\ell_{\text{s}}$ and phase function $p({\bf u},{\bf u}^{\prime})$ [Eq. (38)].
Looking at the radiative flux vector leads to the definition of the transport
mean free path $\ell_{\text{t}}$ [Eq. (41)], related to the intensity vertex
via Eqs. (38) and (39).

The theoretical framework described here will now be used to get closed-form
expressions for the different optical length scales in the cases of random
media described by a continuous permittivity [Sec. II.2] and assemblies of
identical particles [Sec. II.3].

###  II.2 Media with fluctuating continuous permittivity

####  II.2.1 First insights on structural correlations

We consider a statistically homogeneous and isotropic disordered medium,
described by a spatially-dependent permittivity $\epsilon({\bf
r})=\langle\epsilon\rangle+\Delta\epsilon({\bf r})$, where $\Delta\epsilon$ is
the fluctuating part with statistics

$\displaystyle\langle\Delta\epsilon({\bf r})\rangle$ $\displaystyle=$
$\displaystyle 0,$ (45) $\displaystyle\langle\Delta\epsilon({\bf
r})\Delta\epsilon({\bf r}^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\langle\epsilon\rangle^{2}\delta_{\epsilon}^{2}h_{\epsilon}(|{\bf
r}-{\bf r}^{\prime}|).$ (46)

Here,
$\delta_{\epsilon}^{2}=\langle\Delta\epsilon^{2}\rangle/\langle\epsilon\rangle^{2}$
is the normalized variance of $\epsilon$ and $h_{\epsilon}(|{\bf r}-{\bf
r}^{\prime}|)=\langle\Delta\epsilon({\bf r})\Delta\epsilon({\bf
r}^{\prime})\rangle/\langle\Delta\epsilon^{2}\rangle$ is the normalized
permittivity-permittivity correlation function ($h_{\epsilon}(0)=1$).
Hereafter, we assume ergodicity, such that the ensemble average is equivalent
to a volume average over in the infinite-volume limit, and isotropic
permittivity fluctuations, keeping in mind, however, that anisotropic
fluctuations may take place even in isotropic materials Landau  _et al._
(2013).

The statistical properties on $\epsilon({\bf r})$ straightforwardly translate
into statistical properties on ${\bm{\mathbf{\mathcal{X}}}}({\bf r})$ via Eq.
(3). The self-energy $\bm{\mathbf{\mathit{\Sigma}}}$ can be expressed in terms
of ${\bm{\mathbf{\mathcal{X}}}}$ by inserting the expression for the
transition operator ${\bm{\mathbf{\mathcal{T}}}}$ given by Eq. (18) into Eq.
(20), leading to

$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}=\left\langle{\bm{\mathbf{\mathcal{X}}}}\left[\bm{\mathbf{1}}-{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}\right]^{-1}\right\rangle\left\langle\left[\bm{\mathbf{1}}-{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}\right]^{-1}\right\rangle^{-1}.$
(47)

In the simplest approach, we proceed by assuming that the scattering potential
${\bm{\mathbf{\mathcal{X}}}}$ weakly fluctuates around its average value
$\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle$. Expanding the last
expression near $\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle$ leads
to

$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}\sim\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle+\left\langle\left({\bm{\mathbf{\mathcal{X}}}}-\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle\right){\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\left({\bm{\mathbf{\mathcal{X}}}}-\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle\right)\right\rangle+\cdots$
(48)

At this stage, we need to define explicitly the constant auxiliary background
permittivity $\epsilon_{\text{b}}$, which describes the reference value around
which the permittivity fluctuates. A reasonable choice is to set it to the
average permittivity,
$\epsilon_{\text{b}}=\left\langle\epsilon\right\rangle\equiv\epsilon_{\text{av}}$,
which, for a two-component medium with permittivities $\epsilon_{\text{p}}$
and $\epsilon_{\text{h}}$ at filling fractions $f$ and $1-f$, respectively,
would simply be
$\epsilon_{\text{av}}=f\epsilon_{\text{p}}+(1-f)\epsilon_{\text{h}}$. Having
then $\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle=0$, the leading
term for the self-energy becomes
$\left\langle{\bm{\mathbf{\mathcal{X}}}}{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}\right\rangle$,
such that

$\displaystyle\bm{\mathbf{\Sigma}}({\bf r}-{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle k_{\text{av}}^{4}\delta_{\epsilon}^{2}h_{\epsilon}(|{\bf
r}-{\bf r}^{\prime}|){\bf G}_{\text{av}}({\bf r}-{\bf r}^{\prime}),$ (49)

where ${\bf G}_{\text{av}}$ is the Green tensor in a homogeneous medium with
permittivity $\epsilon_{\text{av}}$. Correlated permittivity fluctuations
mutually interacting via ${\bf G}_{\text{av}}$ are readily responsible for the
non-local character of the self-energy [Eq. (6)]. In Fourier space, Eq. (49)
becomes

$\displaystyle\bm{\mathbf{\Sigma}}({\bf k})$ $\displaystyle=$ $\displaystyle
k_{\text{av}}^{4}\delta_{\epsilon}^{2}\int h_{\epsilon}(|{\bf k}-{\bf
k}^{\prime}|){\bf G}_{\text{av}}({\bf k}^{\prime})\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}}.$ (50)

The extinction mean free path $\ell_{\text{e}}$ can finally be determined
using Eq. (13) with $k_{\text{r}}=k_{\text{av}}$. We consider media composed
of non-absorbing materials [$\epsilon({\bf r})$ always real], such that
$k_{\text{av}}=k_{0}\sqrt{\epsilon_{\text{av}}}$. The imaginary part of the
Green tensor is given by

$\displaystyle\text{Im}{\bf G}_{\text{av}}({\bf
k})=\pi\delta\left(k^{2}-k_{\text{av}}^{2}\right){\bf P}({\bf u}).$ (51)

Using Eq. (10) for the transverse component with ${\bf e}$ the polarization
such that ${\bf e}\cdot{\bf u}=0$, we eventually obtain

$\displaystyle\frac{1}{\ell_{\text{e}}}$ $\displaystyle\simeq$
$\displaystyle\frac{\text{Im}\Sigma_{\perp}(k_{\text{av}})}{k_{\text{av}}},$
$\displaystyle=$
$\displaystyle\frac{k_{\text{av}}^{4}}{16\pi^{2}}\delta_{\epsilon}^{2}\int
h_{\epsilon}(k_{\text{av}}|{\bf u}-{\bf u}^{\prime}|)\left({\bf e}\cdot{\bf
P}({\bf u}^{\prime}){\bf e}\right)d{\bf u}^{\prime}.$

Introducing the scattering wavenumber $q=k_{\text{av}}|{\bf u}-{\bf
u}^{\prime}|$, we eventually reach

$\displaystyle\frac{1}{\ell_{\text{e}}}=\frac{k_{0}^{4}}{8\pi
k_{\text{av}}^{2}}\epsilon_{\text{av}}^{2}\delta_{\epsilon}^{2}\int_{0}^{2k_{\text{av}}}P\left(\frac{q}{2k_{\text{av}}}\right)h_{\epsilon}(q)qdq,$
(53)

where

$\displaystyle P(k)\equiv 1-2k^{2}+2k^{4},$ (54)

is due to the vector nature of light. The integral over $q$, which corresponds
to an integral over the scattering angle, should be made up to twice the wave
number in the homogeneous effective background with permittivity
$\epsilon_{\text{av}}$.

Equation (53), obtained in non-absorbing disordered media with weak
permittivity fluctuations, constitute the first analytical expression for the
extinction mean free path in correlated media. Very importantly, it shows that
spatial correlations, described here by $h_{\epsilon}$, are as important for
light scattering as the amplitude of the permittivity fluctuations,
$\epsilon_{\text{av}}^{2}\delta_{\epsilon}^{2}=\left\langle\Delta\epsilon^{2}\right\rangle$.

####  II.2.2 Lorentz local fields

The approximation of weak fluctuations is prohibitive in many realistic cases.
Fortunately, this constraint can be eliminated by properly handling the
singularity of the dyadic Green function at the origin [Eq. (14)], which
constitutes the basis of a strong fluctuation theory Finkel’berg (1964);
Ryzhov  _et al._ (1965); Tsang and Kong (1981). Related approaches were
introduced by Bedeaux and Mazur (1973) and Felderhof (1974) in the description
of the optical response of non-polar fluids, and more recently by Torquato and
Kim (2021) in the study of two-phase composite media. Extensions to chiral and
anisotropic media have also been proposed Ryzhov and Tamoikin (1970); Michel
and Lakhtakia (1995); Mackay and Lakhtakia (2015), but will not be discussed
here.

The singularity of the Green function is handled by considering the scattering
medium as being made of infinitesimal volume elements within which the
polarization density ${\bf P}({\bf r})$ is constant. This physical viewpoint
is nothing but the basis of the renowned discrete dipole approximation Draine
and Flatau (1994); Lakhtakia (1992). Let us then write the Green function as a
sum of two contributions

$\displaystyle{\bf G}_{\text{b}}({\bf r},{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle\Theta(|{\bf r}-{\bf r}^{\prime}|-a)\tilde{{\bf
G}}_{\text{b}}({\bf r},{\bf r}^{\prime})$ (55) $\displaystyle+$
$\displaystyle\Theta(a-|{\bf r}-{\bf
r}^{\prime}|)\bm{\mathbf{g}}_{\text{b}}({\bf r},{\bf r}^{\prime}),$

where $\bm{\mathbf{g}}_{\text{b}}$ contains the singular part of the Green
function and $\tilde{{\bf G}}_{\text{b}}$ is the so-called Lorentz propagator,
which is then purely non-local. The contributions are distinguished as
belonging or not to a spherical volume of radius $a$ and volume $v=4\pi
a^{3}/3$ around ${\bf r}^{\prime}$. Choosing $a$ such that $k_{\text{b}}a\ll
1$, the singular part reads

$\displaystyle\left.\bm{\mathbf{{g}}}_{\text{b}}({\bf r},{\bf
r}^{\prime})\right|_{k_{\text{b}}a\ll 1}$ $\displaystyle=$
$\displaystyle-\frac{1}{3k_{\text{b}}^{2}}\delta({\bf r}-{\bf
r}^{\prime})\bm{\mathbf{1}}+{\rm
i}\frac{k_{\text{b}}}{6\pi}\bm{\mathbf{1}}+\cdots$ (56)

Keeping the lowest-order terms in the real and imaginary parts, the Lippmann-
Schwinger equation [Eq. (15)] can be rewritten as

$\displaystyle{\bf E}({\bf r})$ $\displaystyle=$ $\displaystyle\left\\{{\bf
E}_{\text{b}}({\bf r})+\left(-\frac{1}{3k_{\text{b}}^{2}}+{\rm
i}\frac{k_{\text{b}}v}{6\pi}\right){\bm{\mathbf{\mathcal{X}}}}({\bf r}){\bf
E}({\bf r})\right\\}$ (57) $\displaystyle+$ $\displaystyle\int\tilde{{\bf
G}}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bm{\mathbf{\mathcal{X}}}}({\bf
r}^{\prime}){\bf E}({\bf r}^{\prime})d{\bf r}^{\prime}.$

The actual field at ${\bf r}$ is then given by the sum of the external field
${\bf E}_{\text{b}}({\bf r})$ and the local contributions (terms in braces),
plus the non-local contributions coming from neighboring permittivity
fluctuations (integral term).

In this framework, the field exciting a small volume element around ${\bf r}$,
${\bf E}_{\text{exc}}({\bf r})$, is the sum of the incident (background) field
and the field scattered by other permittivity fluctuations. Going back to the
operator notation, we thus reach an important set of equalities

$\displaystyle{\bf E}_{\text{exc}}$ $\displaystyle=$ $\displaystyle{\bf
E}_{\text{b}}+\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}{\bf
E}=\left[\bm{\mathbf{1}}-\bm{\mathbf{\it
g}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}\right]{\bf E}$ (58) $\displaystyle=$
$\displaystyle{\bf
E}_{\text{b}}+\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}{\bf
E}_{\text{exc}}=\left[\bm{\mathbf{1}}+\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}{\bm{\mathbf{\mathcal{T}}}}\right]{\bf
E}_{\text{b}}.$

We have introduced here a new quantity, $\tilde{{\bm{\mathbf{\mathcal{T}}}}}$,
that plays the role of a local transition operator. From Eqs. (58), we
straightforwardly obtain

$\displaystyle{\bm{\mathbf{\mathcal{T}}}}=\tilde{{\bm{\mathbf{\mathcal{T}}}}}\left[\bm{\mathbf{1}}-\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right]^{-1}.$
(59)

The transition operator of the system is now expressed as a multiple-
scattering expansion on small independent scattering elements, connected via
the (non-local) Lorentz propagator. Consistently with this picture, from Eqs.
(58), we also obtain

$\displaystyle\tilde{{\bm{\mathbf{\mathcal{T}}}}}={\bm{\mathbf{\mathcal{X}}}}\left[\bm{\mathbf{1}}-\bm{\mathbf{\it
g}}_{\text{b}}{\bm{\mathbf{\mathcal{X}}}}\right]^{-1}.$ (60)

For $k_{\text{b}}a\ll 1$, the (local) transition operator $\tilde{{\bf T}}$ is
directly proportional to the polarizability of the volume element,

$\displaystyle\tilde{{\bf T}}({\bf r},{\bf r}^{\prime})$ $\displaystyle=$
$\displaystyle k_{\text{b}}^{2}\frac{\alpha({\bf r})}{v}\delta({\bf r}-{\bf
r}^{\prime})\bm{\mathbf{1}},$ (61)

with

$\displaystyle\alpha({\bf r})=\frac{\alpha_{0}({\bf r})}{1-{\rm
i}\frac{k_{\text{b}}^{3}}{6\pi}\alpha_{0}({\bf r})},\;\alpha_{0}({\bf
r})=3v\frac{\epsilon({\bf r})-\epsilon_{\text{b}}}{\epsilon({\bf
r})+2\epsilon_{\text{b}}},$ (62)

where $\alpha_{0}$ is the quasi-static polarizability that depends on space
via $\epsilon({\bf r})$.

####  II.2.3 Average exciting field

To determine the self-energy $\bm{\mathbf{\mathit{\Sigma}}}$ of the system, it
is convenient to introduce a self-energy
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}$ for the exciting field defined from a
Dyson equation

$\displaystyle\langle{\bf E}_{\text{exc}}\rangle={\bf
E}_{\text{b}}+{\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}}\tilde{\bm{\mathbf{\mathit{\Sigma}}}}\langle{\bf
E}_{\text{exc}}\rangle,$ (63)

thereby leading to

$\displaystyle\tilde{\bm{\mathbf{\mathit{\Sigma}}}}=\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\left[\bm{\mathbf{1}}-\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right]^{-1}\right\rangle\left\langle\left[\bm{\mathbf{1}}-\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right]^{-1}\right\rangle^{-1}.$
(64)

Note the similarity of this expression with Eq. (47), where the self-energy
$\bm{\mathbf{\mathit{\Sigma}}}$ was expressed directly in terms of the
scattering potential ${\bm{\mathbf{\mathcal{X}}}}$. From Eqs. (20), (59) and
(64), one shows that the two self-energies are related as

$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}=\tilde{\bm{\mathbf{\mathit{\Sigma}}}}\left[\bm{\mathbf{1}}+\bm{\mathbf{\it
g}}_{\text{b}}\tilde{\bm{\mathbf{\mathit{\Sigma}}}}\right]^{-1}.$ (65)

Let us then expand $\tilde{\bm{\mathbf{\mathit{\Sigma}}}}$ in Eq. (64) near
$\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\rangle=\tilde{{\bm{\mathbf{\mathcal{T}}}}}-\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}$,

$\displaystyle\tilde{\bm{\mathbf{\mathit{\Sigma}}}}\sim\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle+\left\langle\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle+\cdots\equiv\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}+\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{2}+\cdots.$
(66)

We have introduced here a new “dressed” propagator,

$\displaystyle\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}=\left[\bm{\mathbf{1}}-\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle\right]^{-1}\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}.$
(67)

Noting that $\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\rangle$ corresponds to
an average polarizability of the medium (for $k_{\text{b}}a\ll 1$, see Eqs.
(61)), one understands that $\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}$
describes the fact that the fields propagate from fluctuation to fluctuation
via a medium with a permittivity that can differ from the background
permittivity $\epsilon_{\text{b}}$ Bedeaux and Mazur (1973); Felderhof (1974).

The self-energy (for the average exciting field) in Eq. (66) now explicitly
depends on the spatial correlations of the fluctuations of
$\tilde{{\bm{\mathbf{\mathcal{T}}}}}$, i.e. of the polarizability of small
volume elements. In most practical cases, the expansion is stopped at the
second order, corresponding to the so-called “bilocal” approximation Tsang and
Kong (2001), due to the lack of information on higher-order correlation
functions in real systems.

Expanding $\bm{\mathbf{\mathit{\Sigma}}}$ in Eq. (65) near
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}$, we obtain

$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}$ $\displaystyle\sim$
$\displaystyle\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}\left[\bm{\mathbf{1}}+\bm{\mathbf{\it
g}}_{\text{b}}\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}\right]^{-1}+\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{2}\left[\bm{\mathbf{1}}+\bm{\mathbf{\it
g}}_{\text{b}}\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}\right]^{-2}+\cdots$
(68) $\displaystyle\equiv$
$\displaystyle\bm{\mathbf{\mathit{\Sigma}}}_{1}+\bm{\mathbf{\mathit{\Sigma}}}_{2}+\cdots$

The self-energy is now expressed in terms of the scattering properties of
vanishingly small, individual scattering elements.

####  II.2.4 Long-wavelength solutions for strongly fluctuating media:
Bruggeman versus Maxwell-Garnett

As in the case of weakly fluctuating media discussed above, we now need to
give an explicit definition of the constant auxiliary background permittivity
$\epsilon_{\text{b}}$, that describes the homogeneous effective medium in
which the permittivity fluctuations scatter light. We will see that this sole
parameter constitutes the essential difference between the two most famous
“mixing rules” in the literature, due to Bruggeman Bruggeman (1935) and
Maxwell Garnett Maxwell Garnett (1904), presented here in a unique theoretical
framework. Despite the arbitrariness in the choice of $\epsilon_{\text{b}}$,
it is important to realize that all models would eventually be strictly
equivalent when carried out to infinite order. The applicability of a model is
thus mainly a question of accuracy at the lowest orders and convergence.

A first possibility for $\epsilon_{\text{b}}$ is to set it such that
$\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\rangle=0$. In the limit of small
volume elements, this corresponds to having a zero average polarizability, see
Eqs. (61)-(62). Considering a two-component medium with relative
permittivities $\epsilon_{\text{p}}$ (at filling fraction $f$) and
$\epsilon_{\text{h}}$ (at filling fraction $1-f$) in the quasi-static limit
($\alpha=\alpha_{0}$ in Eq. (62)), one obtains

$\displaystyle\frac{\epsilon_{\text{p}}-\epsilon_{\text{BG}}}{\epsilon_{\text{p}}+2\epsilon_{\text{BG}}}f+\frac{\epsilon_{\text{h}}-\epsilon_{\text{BG}}}{\epsilon_{\text{h}}+2\epsilon_{\text{BG}}}(1-f)=0$
(69)

which is nothing but the Bruggeman mixing rule Bruggeman (1935) with
$\epsilon_{\text{b}}\equiv\epsilon_{\text{BG}}$. The generalization to
$N$-component media is straightforward. Having
$k_{\text{b}}^{2}=k_{\text{BG}}^{2}=k_{0}^{2}\epsilon_{\text{BG}}$,
$\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{b}}=\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{BG}}$
and $\bm{\mathbf{\Sigma}}\sim\tilde{\bm{\mathbf{\Sigma}}}$ since
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}=0$, we eventually arrive to

$\displaystyle\bm{\mathbf{\Sigma}}({\bf
k})=k_{\text{BG}}^{4}\delta_{\alpha}^{2}\int h_{\alpha}(|{\bf k}-{\bf
k}^{\prime}|)\tilde{{{\bf G}}}_{\text{BG}}({\bf k}^{\prime})\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}},$ (70)

with $\delta_{\alpha}^{2}=\langle\Delta\alpha^{2}\rangle/v^{2}$ a normalized
variance of the polarizability and $h_{\alpha}(|{\bf k}-{\bf k}^{\prime}|)$
the Fourier transform of the normalized polarizability-polarizability
correlation function $h_{\alpha}(|{\bf r}-{\bf
r}^{\prime}|)=\langle\Delta\alpha({\bf r})\Delta\alpha({\bf
r}^{\prime})\rangle/\langle\Delta\alpha^{2}\rangle$. The function $h_{\alpha}$
therefore plays the role of $h_{\epsilon}$ in the weakly fluctuating
permittivity model to describe structural correlations.

Assuming non-absorbing media and following the same steps as those leading to
Eq. (53) with $\text{Im}\tilde{{{\bf G}}}_{\text{BG}}({\bf
k})=\pi\delta\left(k^{2}-k_{\text{BG}}^{2}\right){\bf P}({\bf u})$, we obtain

$\displaystyle\frac{1}{\ell_{\text{e}}}=\frac{k_{0}^{4}}{8\pi
k_{\text{BG}}^{2}}\epsilon_{\text{BG}}^{2}\delta_{\alpha}^{2}\int_{0}^{2k_{\text{BG}}}P\left(\frac{q}{2k_{\text{BG}}}\right)h_{\alpha}(q)qdq,$
(71)

with $q=k_{\text{BG}}|{\bf u}-{\bf u}^{\prime}|$. Equation (71) is strikingly
similar to Eq. (53), the essential differences being (i) on the permittivity
of the homogeneous effective medium and (ii) on the description of the system
via a local polarizability instead of a local permittivity.

A second possibility for the choice of $\epsilon_{\text{b}}$ is to set it to
the permittivity of the host medium, i.e.
$\epsilon_{\text{b}}=\epsilon_{\text{h}}$, in which case
$\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\rangle\neq 0$. Considering again a
two-component system for volume elements in the quasi-static limit, we obtain

$\displaystyle\tilde{\bm{\mathbf{\Sigma}}}_{1}({\bf r}-{\bf
r}^{\prime})=k_{\text{h}}^{2}\rho\alpha_{0}\delta({\bf r}-{\bf
r}^{\prime})\bm{\mathbf{1}},$ (72)

with $\rho=f/v$ the average number density of the small volume elements with
permittivity $\epsilon_{\text{p}}$, and

$\displaystyle\tilde{\bm{\mathbf{\Sigma}}}_{2}({\bf r}-{\bf
r}^{\prime})=k_{\text{h}}^{4}\delta_{\alpha}^{2}h_{\alpha}(|{\bf r}-{\bf
r}^{\prime}|)\hat{{{\bf G}}}_{\text{h}}({\bf r}-{\bf r}^{\prime}).$ (73)

Using Eq. (68), we arrive to an expression for the self-energy in reciprocal
space as

$\displaystyle\bm{\mathbf{\Sigma}}({\bf k})$ $\displaystyle=$ $\displaystyle
k_{0}^{2}\left(\epsilon_{\text{MG}}-\epsilon_{\text{h}}\right)\bm{\mathbf{1}}$
(74) $\displaystyle+$ $\displaystyle
k_{\text{h}}^{4}\delta_{\alpha}^{2}f_{\text{L}}\int h_{\alpha}(|{\bf k}-{\bf
k}^{\prime}|)\hat{{{\bf G}}}_{\text{h}}({\bf k}^{\prime})\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}},$

where
$f_{\text{L}}=(\partial\epsilon_{\text{MG}}/\partial\rho)/(\epsilon_{\text{h}}\alpha_{0})$
and

$\displaystyle\epsilon_{\text{MG}}=\epsilon_{\text{h}}+\epsilon_{\text{h}}\frac{\rho\alpha_{0}}{1-\rho\alpha_{0}/3},$
(75)

is the Maxwell-Garnett mixing rule Maxwell Garnett (1904); Markel (2016). By
contrast with the previous case, the lowest-order term now provides the
renormalization of the wave number in the effective medium, structural
correlations appearing at the next order. The factor $f_{\text{L}}$ is a
local-field correction coming from the fact that the fluctuation of
polarizability ($\Delta\alpha$ in $\delta_{\alpha}^{2}$) is evaluated with
respect to the host medium. Following again the same steps as those leading to
Eq. (53), noting that $\epsilon_{\text{MG}}$ is real in the quasi-static limit
for non-absorbing media, and showing that

$\displaystyle\text{Im}\hat{{{\bf G}}}_{\text{h}}({\bf k})=\pi
f_{\text{L}}\delta\left(k^{2}-k_{\text{MG}}^{2}\right){\bf P}({\bf u}),$ (76)

see Appendix A.2, we eventually obtain

$\displaystyle\frac{1}{\ell_{\text{e}}}=\frac{k_{0}^{4}}{8\pi
k_{\text{MG}}^{2}}f_{\text{L}}^{2}\epsilon_{\text{h}}^{2}\delta_{\alpha}^{2}\int_{0}^{2k_{\text{MG}}}P\left(\frac{q}{2k_{\text{MG}}}\right)h_{\alpha}(q)qdq,$
(77)

with $q=k_{\text{MG}}|{\bf u}-{\bf u}^{\prime}|$. The expression for
$\ell_{\text{e}}$ takes the same form as previously with the permittivity of
the homogeneous medium around the fluctuations being given by the Maxwell-
Garnett mixing rule and a prefactor to account for local-field corrections.

All in all, the similarity between Eqs. (53), (71) and (77) demonstrates the
deep physical implication of structural correlations for light scattering,
whatever the effective medium approach used.

####  II.2.5 Expressions for the scattering and transport mean free paths from
the average intensity

We conclude this part on continuous permittivity media by deriving the
expressions for $\ell_{\text{s}}$ and $\ell_{\text{t}}$ from the theory for
the average intensity. An expansion of the intensity vertex
$\bm{\mathbf{\mathit{\Gamma}}}$ in Eq. (23) with ${\bm{\mathbf{\mathcal{T}}}}$
given by Eq. (18) to the first two orders leads to

$\displaystyle\bm{\mathbf{\mathit{\Gamma}}}\sim\left\langle{\bm{\mathbf{\mathcal{T}}}}{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle-\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle\left\langle{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle\sim\left\langle{\bm{\mathbf{\mathcal{X}}}}{\bm{\mathbf{\mathcal{X}}}}^{*}\right\rangle-\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle\left\langle{\bm{\mathbf{\mathcal{X}}}}^{*}\right\rangle.$
(78)

This expansion is valid in the weak extinction limit
$k_{\text{r}}\ell_{\text{e}}\gg 1$. Similarly to the case of weakly
fluctuating media, we set the auxiliary background permittivity as
$\epsilon_{\text{b}}=\left\langle\epsilon\right\rangle\equiv\epsilon_{\text{av}}$,
such that $\left\langle{\bm{\mathbf{\mathcal{X}}}}\right\rangle=0$, and assume
a non-absorbing material. This leads to

$\bar{\bm{\mathbf{\Gamma}}}(k_{\text{av}}{\bf u},k_{\text{av}}{\bf
u}^{\prime},k_{\text{av}}{\bf u},k_{\text{av}}{\bf
u}^{\prime})=k_{0}^{4}\epsilon_{\text{av}}^{2}\delta_{\epsilon}^{2}h_{\epsilon}(k_{\text{av}}|{\bf
u}-{\bf u}^{\prime}|)\\\ \times\bm{\mathbf{1}}\otimes\bm{\mathbf{1}}.$ (79)

The scattering mean free path is then obtained by integrating Eq. (38) over
${\bf u}^{\prime}$, leading to

$\displaystyle\frac{1}{\ell_{\text{s}}}=\frac{k_{0}^{4}}{8\pi
k_{\text{av}}^{2}}\epsilon_{\text{av}}^{2}\delta_{\epsilon}^{2}\int_{0}^{2k_{\text{av}}}P\left(\frac{q}{2k_{\text{av}}}\right)h_{\epsilon}(q)q\mathrm{d}q,$
(80)

with $q=k_{\text{av}}|{\bf u}-{\bf u}^{\prime}|$. Similarly, the transport
mean free is obtained by calculating the average cosine of Eq. (38) and using
Eqs. (41)-(42),

$\displaystyle\frac{1}{\ell_{\text{t}}}=\frac{k_{0}^{4}}{16\pi
k_{\text{av}}^{4}}\epsilon_{\text{av}}^{2}\delta_{\epsilon}^{2}\int_{0}^{2k_{\text{av}}}P\left(\frac{q}{2k_{\text{av}}}\right)h_{\epsilon}(q)q^{3}\mathrm{d}q.$
(81)

In absence of absorption, one should have $\ell_{\text{s}}=\ell_{\text{e}}$
which is exactly what we obtain here comparing Eqs. (80) and (53).

###  II.3 Particulate media

####  II.3.1 Expansion for identical scatterers

Let us now consider in detail a system whose morphology consists in localized
(i.e., compact) permittivity variations in a uniform background. We take the
most natural choice for the background permittivity
$\epsilon_{\text{b}}=\epsilon_{\text{h}}$ from the start but the theory can
also be developed with an arbitrary $\epsilon_{\text{b}}$. We also restrict
the discussion to composite media made of identical inclusions with relative
permittivity $\epsilon_{\text{p}}$ circumscribed in a volume $v$, centered at
positions ${\bf R}=\left[{\bf R}_{1},{\bf R}_{2},\dots\bf R_{N}\right]$. The
medium permittivity then reads

$\displaystyle\epsilon({\bf r})=\sum_{j}\epsilon_{\text{p}}({\bf r}-{\bf
R}_{j})\Theta(a-|{\bf r}-{\bf R}_{j}|).$ (82)

The particle configuration is described statistically by the probability
distribution function $p({\bf R})$. Implicitly, we neglect here the
possibility to have orientational correlations between particles (otherwise,
the distribution should include orientational variables). Under the ergodic
hypothesis, defining the statistical average as an average over all possible
particle positions as $\langle\mathbf{f}({\bf R})\rangle=\int\mathbf{f}({\bf
R})p({\bf R})d{\bf R}$, where $\mathbf{f}({\bf R})$ is an arbitrary tensor,
the statistical properties of the medium can be described by $n$-particle
probability density functions Tsang  _et al._ (2004); Lebowitz and Percus
(1963)

$\displaystyle\rho_{n}({\bf r}_{1},\cdots{\bf r}_{n})=\big{\langle}$
$\displaystyle\sum_{j_{1}\neq j_{2}\cdots\neq j_{n}}$
$\displaystyle\delta({\bf r}_{1}-{\bf R}_{j_{1}})$ (83) $\displaystyle\cdots$
$\displaystyle\delta({\bf r}_{n}-{\bf R}_{j_{n}})\big{\rangle},$

or equivalently by $n$-particle correlation functions

$\displaystyle g_{n}({\bf r}_{1},\cdots{\bf
r}_{n})=\frac{1}{\rho^{n}}\rho_{n}({\bf r}_{1},\cdots{\bf r}_{n}),$ (84)

where $\rho$ is the constant particle number density reached in the limit of
infinite system size ($\rho=\lim_{N,V\rightarrow\infty}N/V$).

Similar to the case of random media described by a continuous permittivity,
the first step is to derive an expression for the transition operator
${\bm{\mathbf{\mathcal{T}}}}$ of the medium. The fact to have a discrete set
of identical particles allows us to express the multiple-scattering problem in
such a way as to separate the effects associated to (local) particulate
resonances and (non-local) structural correlations on light scattering. We
start by rewritting the integral equation for the total field, Eq. (16) using
integral operators, for particulate media,

$\displaystyle{\bf E}={\bf
E}_{\text{h}}+\sum_{j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E},$ (85)

with the effective scattering potential ${\bm{\mathbf{\mathcal{X}}}}_{j}({\bf
r}-{\bf R}_{j})\equiv k_{0}^{2}[\epsilon_{\text{p}}({\bf r}-{\bf
R}_{j})\Theta(a-|{\bf r}-{\bf R}_{j}|)-\epsilon_{\text{h}}]\bm{\mathbf{1}}$.
We can then express the polarization induced in particle $j$ in terms of the
polarization induced in all particles as

$\displaystyle{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf E}$ $\displaystyle=$
$\displaystyle{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E}_{\text{h}}+{\bm{\mathbf{\mathcal{X}}}}_{j}\sum_{k}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{k}{\bf
E},$ (86) $\displaystyle=$ $\displaystyle{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E}_{\text{h}}+{\bm{\mathbf{\mathcal{X}}}}_{j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E}+{\bm{\mathbf{\mathcal{X}}}}_{j}\sum_{k\neq
j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{k}{\bf
E},$ (87) $\displaystyle=$ $\displaystyle{\bm{\mathbf{\mathcal{T}}}}_{j}{\bf
E}_{\text{h}}+{\bm{\mathbf{\mathcal{T}}}}_{j}\sum_{k\neq
j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{k}{\bf
E}.$ (88)

In the second expression, we separated the contribution from the particle $j$
on itself and the contribution from all other particles. The last expression
was obtained by introducing the transition operator
${\bm{\mathbf{\mathcal{T}}}}_{j}$ of an individual particle centered at ${\bf
R}_{j}$ as

$\displaystyle{\bm{\mathbf{\mathcal{T}}}}_{j}={\bm{\mathbf{\mathcal{X}}}}_{j}\left[\bm{\mathbf{1}}-{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{j}\right]^{-1}.$
(89)

Very importantly, ${\bm{\mathbf{\mathcal{T}}}}_{j}$ can be determined for
particles of virtually any size, shape and composition, either analytically
using Mie theory for simple geometries like spherical particles Bohren and
Huffman (2008) or numerically by any method for solving Maxwell’s equations
otherwise Mishchenko  _et al._ (2000). This allows considering strongly
resonant particles exhibiting high-order multipolar resonances as the building
blocks of the disordered medium.

Inserting Eq. (88) into Eq. (85) and iterating over scattering sequences, we
reach an expression for the transition operator ${\bm{\mathbf{\mathcal{T}}}}$
of the entire system [Eq. (17)] in terms of the transition operator of the
individual particle, as

$\displaystyle{\bm{\mathbf{\mathcal{T}}}}$ $\displaystyle=$
$\displaystyle\sum_{j}{\bm{\mathbf{\mathcal{T}}}}_{j}+\sum_{j}{\bm{\mathbf{\mathcal{T}}}}_{j}\sum_{k\neq
j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{T}}}}_{k}$ (90)
$\displaystyle+$
$\displaystyle\sum_{j}{\bm{\mathbf{\mathcal{T}}}}_{j}\sum_{k\neq
j}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{T}}}}_{k}\sum_{l\neq
k}{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}{\bm{\mathbf{\mathcal{T}}}}_{l}+\cdots$

This series expansion is the root of multiple scattering theory for
particulate media and was introduced in the pioneering work of Kirkwood and
Yvon as early as 1936 Kirkwood (1936); Yvon (1937) to determine the
permittivity of molecular liquids. Similar multiple scattering equations were
later discussed by Foldy and Lax Foldy (1945); Lax (1951). It is further
interesting to note the occurrence of so-called “recurrent scattering”, that
is, scattering sequences that involve the same particle multiple times (for
instance, $l$ can be equal to $j$ in the last displayed term). The series in
Eq. (90) can in fact be rewritten as a “cluster” expansion according to the
number of particles involved in each term Finkel’berg (1964); Felderhof  _et
al._ (1982).

It is interesting to remark that the Green function in Eq. (89) for the
transition operator ${\bm{\mathbf{\mathcal{T}}}}_{j}$ of a specific particle
$j$ always connects two points that belong to the same particle, whereas the
Green function in Eq. (90) for the transition operator
${\bm{\mathbf{\mathcal{T}}}}$ of the entire system always connects two points
that belong to different particles. This is conceptually analogous to Eq. (55)
for continuous permittivity media where the Green function was split into
local and non-local terms. To determine the self-energy
$\bm{\mathbf{\mathit{\Sigma}}}$ in Dyson equation [Eq. (19)], we may then
follow the same strategy as done for continuous permittivity media and
consider the exciting field ${\bf E}_{\text{exc}}$. Let us then write the
Green function ${\bm{\mathbf{\mathcal{G}}}}_{\text{h}}$ as $\bm{\mathbf{\it
g}}_{\text{h}}$ when connecting two points in the same particle and
$\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}$ otherwise. Removing the local
contribution on the induced polarization in Eq. (87) and rewriting
${\bm{\mathbf{\mathcal{X}}}}_{k}{\bf E}$ in terms of the exciting field leads
to

$\displaystyle{\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E}_{\text{exc}}={\bm{\mathbf{\mathcal{X}}}}_{j}{\bf
E}_{\text{h}}+{\bm{\mathbf{\mathcal{X}}}}_{j}\sum_{k}\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}_{k}{\bf
E}_{\text{exc}},$ (91)

with
$\tilde{{\bm{\mathbf{\mathcal{T}}}}}_{j}={\bm{\mathbf{\mathcal{X}}}}_{j}\left[\bm{\mathbf{1}}-\bm{\mathbf{\it
g}}_{\text{h}}{\bm{\mathbf{\mathcal{X}}}}_{j}\right]^{-1}={\bm{\mathbf{\mathcal{T}}}}_{j}$,
see Eq. (89). Note that the sum now runs over all particles $k$. Further
defining
$\tilde{{\bm{\mathbf{\mathcal{T}}}}}\equiv\sum_{j}\tilde{{\bm{\mathbf{\mathcal{T}}}}}_{j}$,
we naturally recover Eq. (59) for continuous media,

$\displaystyle{\bm{\mathbf{\mathcal{T}}}}=\tilde{{\bm{\mathbf{\mathcal{T}}}}}\left[\bm{\mathbf{1}}-\tilde{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right]^{-1}.$
(92)

Expanding the self-energy $\tilde{\bm{\mathbf{\mathit{\Sigma}}}}$ for the
average exciting field near
$\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\rangle=\tilde{{\bm{\mathbf{\mathcal{T}}}}}-\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}$
then leads to Eq. (66), and expanding $\bm{\mathbf{\mathit{\Sigma}}}$ near
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}$ to Eq. (68).

The problem of scattering by particulate media is thus described in a strictly
similar manner to that of scattering by strongly fluctuating continuous
permittivity media (i.e., including local-field corrections). The essential
difference is the fact that the volume elements composing the continuous
medium are no longer vanishingly small but actual scattering particles of a
finite size.

Recall that Eq. (66) is reached by neglecting particle correlations beyond the
second order. Higher-order particle correlations may yet be taken into account
by treating them as sequences of two-particle correlations (which is formally
exact for crystalline media). This is the so-called quasicrystalline
approximation (QCA) originally due to Lax (1952), further developed by
Fikioris and Waterman (1964) and Tsang and Kong (1980, 1982) among others, and
used nowadays in various contexts Tsang  _et al._ (2000); Kristensson (2015);
Wang and Zhao (2018); Gower  _et al._ (2018).

####  II.3.2 Extinction mean free path and effective medium theories

Writing the transition operator of an individual particle as
${\bm{\mathbf{\mathcal{T}}}}_{j}\equiv{\bf T}_{0}({\bf r}-{\bf R}_{j},{\bf
r}^{\prime}-{\bf R}_{j})$, we obtain for the lowest order in Eq. (66),
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{1}\equiv\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle=\left\langle\sum_{j}\tilde{{\bm{\mathbf{\mathcal{T}}}}}_{j}\right\rangle$,

$\displaystyle\tilde{\bm{\mathbf{\Sigma}}}_{1}({\bf r}-{\bf
r}^{\prime})=\rho\int{\bf T}_{0}({\bf r}-{\bf r}_{\text{p}},{\bf
r}^{\prime}-{\bf r}_{\text{p}})d{\bf r}_{\text{p}}.$ (93)

This term leads to the extinction (optical) theorem for uncorrelated particle
assemblies at vanishingly small particle densities, which relates the
extinction mean free path to the forward scattering characteristics of the
individual particle, $\ell_{\text{e}}^{-1}\sim\rho\left({\bf
e}\cdot\text{Im}[{\bf T}_{0}(k_{\text{h}}{\bf u},k_{\text{h}}{\bf u})]{\bf
e}\right)/k_{\text{h}}$.

For the second order,
$\tilde{\bm{\mathbf{\mathit{\Sigma}}}}_{2}\equiv\left\langle\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}\Delta\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle=\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle-\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle\hat{\bm{\mathbf{\mathcal{G}}}}_{\text{h}}\left\langle\tilde{{\bm{\mathbf{\mathcal{T}}}}}\right\rangle$,
we get

$\displaystyle\tilde{\bm{\mathbf{\Sigma}}}_{2}({\bf r}-{\bf r}^{\prime})$
$\displaystyle=$ $\displaystyle\rho\int\left[\delta(|{\bf r}_{\text{p}}-{\bf
r}_{\text{q}}|)+\rho h_{2}(|{\bf r}_{\text{p}}-{\bf r}_{\text{q}}|)\right]$
(94) $\displaystyle\times$ $\displaystyle{\bf T}_{0}({\bf r}-{\bf
r}_{\text{p}},{\bf r}^{\prime\prime}-{\bf r}_{\text{p}})\hat{{{\bf
G}}}_{\text{h}}({\bf r}^{\prime\prime}-{\bf r}^{\prime\prime\prime})$
$\displaystyle\times$ $\displaystyle{\bf T}_{0}({\bf
r}^{\prime\prime\prime}-{\bf r}_{\text{q}},{\bf r}^{\prime}-{\bf
r}_{\text{q}})d{\bf r}^{\prime\prime}d{\bf r}^{\prime\prime\prime}d{\bf
r}_{\text{p}}d{\bf r}_{\text{q}},$

where

$\displaystyle h_{2}({\bf r})\equiv g_{2}({\bf r})-1$ (95)

is the total pair correlation function and $g_{2}$ is defined from Eq. (84).

To reach a general expression for the self-energy
$\bm{\mathbf{\bm{\mathbf{\mathit{\Sigma}}}}}$ via Eq. (68), we need to define
the operator $\bm{\mathbf{\it g}}_{\text{h}}$ that describes the local
propagation of radiation within each particle. As first pointed out by
Sullivan and Deutch Sullivan and Deutch (1976), $\bm{\mathbf{\it
g}}_{\text{h}}$ may be chosen to obtain different lowest order results for the
effective permittivities and refractive index, such as Maxwell-Garnett Bedeaux
and Mazur (1973); Felderhof (1974), Onsager-Bütcher Onsager (1936); Böttcher
_et al._ (1978); Hynne and Bullough (1987) or Wertheim Wertheim (1973).
Differences vanish when all orders of the expansion are taken into account,
but the rate of convergence of the different formulations is influenced by the
particular choice of the Green function Geigenmüller and Mazur (1986);
Sullivan and Deutch (1976); Bedeaux and Mazur (1973); Bedeaux  _et al._
(1987).

For pedagogical reasons, we restrict ourselves here to the Maxwell-Garnett
result obtained in the long-wavelength limit. For some applications dealing
with particles with complex shapes and relatively low scattering contrast, one
can simplify the problem using the well-known Rayleigh-Gans approximation
Bohren and Huffman (2008) or subsequent generalizations Acquista (1976). In
this framework, the transition operator can be written as

$\displaystyle{\bf T}_{0}({\bf r}-{\bf R}_{j},{\bf r}^{\prime}-{\bf
R}_{j})=k_{h}^{2}\frac{\alpha_{\text{p}}}{v}\Theta\left[a-|{\bf r}-{\bf
R}_{j}|\right]\delta({\bf r}-{\bf r}^{\prime})\bm{\mathbf{1}}$

with $\alpha_{\text{p}}$ the particle polarizability. Inserting this
expression in Eqs. (93) and (94) straightforwardly leads to an expression for
the self-energy that reads in reciprocal space as

$\displaystyle\bm{\mathbf{\Sigma}}({\bf k})$ $\displaystyle=$ $\displaystyle
k_{0}^{2}\left(\epsilon_{\text{MG}}-\epsilon_{\text{h}}\right)\bm{\mathbf{1}}+\rho\frac{k_{\text{h}}^{4}\alpha_{\text{p}}^{2}}{\epsilon_{\text{h}}\alpha_{\text{p}}}\frac{\partial\epsilon_{\text{MG}}}{\partial\rho}$
$\displaystyle\times$ $\displaystyle\int S(|{\bf k}-{\bf
k}^{\prime}|)\left|\frac{3j_{1}(|{\bf k}-{\bf k}^{\prime}|a)}{|{\bf k}-{\bf
k}^{\prime}|a}\right|^{2}\hat{{{\bf G}}}_{\text{h}}({\bf
k}^{\prime})\frac{d{\bf k}^{\prime}}{(2\pi)^{3}}.$

$j_{n}$ is the spherical Bessel function of the first kind and order $n$,
$\epsilon_{\text{MG}}$ is the Maxwell-Garnett permittivity given by Eq. (75)
with $\alpha_{0}\equiv\alpha_{\text{p}}$, and

$\displaystyle S({\bf k})\equiv 1+\rho h_{2}({\bf k}),$ (98)

is the static structure factor, probably the most famous quantity for light
scattering studies in correlated disordered media, as we will see throughout
this review.

Following the same steps as those for continuous permittivity random media,
assuming non-absorbing materials, we eventually reach a simple expression for
the extinction mean free path,

$\displaystyle\frac{1}{\ell_{\text{e}}}=\frac{2\pi\rho}{k_{\text{MG}}^{4}}\int_{0}^{2k_{\text{MG}}}F(q)S(q)q\mathrm{d}q,$
(99)

with $q=k_{\text{MG}}|{\bf u}-{\bf u}^{\prime}|$. We have defined here the
form factor

$\displaystyle
F(q)=k_{\text{MG}}^{2}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}\left(\frac{q}{2k_{\text{MG}}}\right)f_{\text{L}}^{2}\left|\frac{3j_{1}(qa)}{qa}\right|^{2},$
(100)

and the Rayleigh differential scattering cross-section

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}(k)=k_{\text{h}}^{4}\frac{\alpha_{\text{p}}^{2}}{(4\pi)^{2}}P(k),$
(101)

where $P(k)$ is given by Eq. (54).

Remarkably, the respective contributions of the individual scattering
elements, via the form factor $F(q)$, and of their spatial arrangement, via
the structure factor $S(q)$, on the extinction (or scattering) strength of the
medium are treated independently. Structural correlations act as a weighting
function to the optical response of a random (uncorrelated) assembly of
identical scatterers, physically due to far-field interference between pairs
of particles.

Equation (99) was obtained here in the long-wavelength limit, for small (non-
resonant) particles. The more general situation of resonant particles is
significantly more difficult to address within the theory for the average
field. A heuristic extension of the Maxwell-Garnett approximation to resonant
dipolar particles has been proposed by Doyle (1989) and further analyzed in
Grimes and Grimes (1991); Ruppin (2000). The approach provides some
understanding to the spectral resonances observed in certain light scattering
experiments. A more rigorous framework is given by the so-called Coherent
Potential Approximation (CPA) Tsang and Kong (1980, 2001), which may be seen
as a generalization of the approach leading to the Bruggeman mixing rule for
continuous permittivity media, as presented above. Considering scattering
elements that are not only the particles but also the host medium, one looks
for an auxiliary background permittivity $\epsilon_{\text{b}}$ such that the
average transition operator vanishes
$\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle=0$, or equivalently that
the background Green function ${\bm{\mathbf{\mathcal{G}}}}_{\text{b}}$ equals
the actual averaged Green function $\langle\bm{\mathbf{\mathcal{G}}}\rangle$.
Different CPA-like models have been developed based on different self-
consistent conditions Soukoulis  _et al._ (1994); Busch and Soukoulis (1995).

####  II.3.3 Scattering and transport mean free paths for resonant scatterers

In this last part, we show that the case of resonant particles can be more
easily solved within the theory for the average intensity. For this, we assume
that a solution for the real part of the effective index, $k_{\text{r}}$,
exists. In the weak extinction limit (i.e. $k_{\text{r}}\ell_{e}\gg 1$) and
using the expansion of the transition operator ${\bm{\mathbf{\mathcal{T}}}}$
in Eq. (90), the intensity vertex given by Eq. (23) can be reduced to its
lowest order terms as

$\displaystyle\bm{\mathbf{\mathit{\Gamma}}}$ $\displaystyle\sim$
$\displaystyle\left\langle{\bm{\mathbf{\mathcal{T}}}}{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle-\left\langle{\bm{\mathbf{\mathcal{T}}}}\right\rangle\left\langle{\bm{\mathbf{\mathcal{T}}}}^{*}\right\rangle$
(102) $\displaystyle\sim$
$\displaystyle\left\langle\sum_{i,j}{\bm{\mathbf{\mathcal{T}}}}_{i}{\bm{\mathbf{\mathcal{T}}}}_{j}^{*}\right\rangle-\left\langle\sum_{i}{\bm{\mathbf{\mathcal{T}}}}_{i}\right\rangle\left\langle\sum_{j}{\bm{\mathbf{\mathcal{T}}}}_{j}^{*}\right\rangle.$

In Fourier space, this leads to

$\bar{\bm{\mathbf{\Gamma}}}(k_{\text{r}}{\bf u},k_{\text{r}}{\bf
u}^{\prime},k_{\text{r}}{\bf u},k_{\text{r}}{\bf u}^{\prime})=\rho{\bf
T}_{0}(k_{\text{r}}{\bf u},k_{\text{r}}{\bf u}^{\prime})\otimes{\bf
T}_{0}^{*}(k_{\text{r}}{\bf u},k_{\text{r}}{\bf u}^{\prime})\\\ \times
S(k_{\text{r}}({\bf u}-{\bf u}^{\prime})).$ (103)

In the radiative transfer limit, taking the on-shell approximation and using
Eqs. (38)-(39) leads to

$\displaystyle\frac{1}{\ell_{\text{s}}}=\frac{\rho}{16\pi^{2}}\int\left|{\bf
P}({\bf u}){\bf T}_{0}(k_{\text{r}}{\bf u},k_{\text{r}}{\bf u}^{\prime}){\bf
e}^{\prime}\right|^{2}S(k_{\text{r}}({\bf u}-{\bf u}^{\prime}))\mathrm{d}{\bf
u}$

where ${\bf e}^{\prime}$ is the polarization vector perpendicular to ${\bf
u}^{\prime}$. This expression is valid for a spherical particle of arbitrary
size. Notice also that ${\bf T}_{0}$ is evaluated for the incident and
scattered wavevectors not in the host medium but in the effective medium, i.e.
with wavenumber $k_{\text{r}}$. Equation (II.3.3) can eventually be
reformulated using the form factor

$\displaystyle
F(q)=k_{\text{r}}^{2}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}(q),$ (105)

where $q=k_{\text{r}}|{\bf u}-{\bf u}^{\prime}|$ and the differential
scattering cross-section is now defined as

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}(q)=\frac{1}{16\pi^{2}}\left|{\bf
P}({\bf u}){\bf T}_{0}(k_{\text{r}}|{\bf u}-{\bf u}^{\prime}|){\bf
e}^{\prime}\right|^{2}.$ (106)

This leads to

$\displaystyle\frac{1}{\ell_{\text{s}}}=\frac{2\pi\rho}{k_{\text{r}}^{4}}\int_{0}^{2k_{\text{r}}}F(q)S(q)q\mathrm{d}q,$
(107)

The approaches based on the average field and the average intensity lead to
the same final expressions [Eqs. (99) and (107), respectively] for non-
absorbing media.

A similar derivation using Eqs. (38), (41) and (42) finally leads to a closed-
form expression of the transport mean free path,

$\displaystyle\frac{1}{\ell_{\text{t}}}=\frac{\pi\rho}{k_{\text{r}}^{6}}\int_{0}^{2k_{\text{r}}}F(q)S(q)q^{3}\mathrm{d}q.$
(108)

Equations (107) and (108) are surely the most widely-used expressions in the
literature on light scattering and transport in correlated disordered media
Fraden and Maret (1990); Rojas-Ochoa  _et al._ (2004); Reufer  _et al._
(2007). These two expressions, originally derived from phenomenological
arguments, have been obtained here within a rigorous theoretical framework.

###  II.4 Summary and further remarks

The literature on multiple light scattering theory is vast and many approaches
have been developed throughout the years. We adopted here a unique theoretical
framework that can handle both families of systems described by a continuous
permittivity that fluctuates in space or by a set of identical particles that
are randomly arranged in space. This has the great benefit of highlighting the
main physical principles behind the role of spatial correlations on light
scattering and transport, as well as the underlying approximations. Analytical
expressions for the characteristic lengths were obtained in the weak
extinction limit ($k_{\text{r}}\ell_{\text{e}}\gg 1$, with $k_{\text{r}}$ the
effective wave number and $\ell_{\text{e}}$ the extinction mean free path),
and considering non-absorbing media such that the extinction is only due to
scattering ($\ell_{\text{e}}=\ell_{\text{s}}$).

For media described by a continuous permittivity, the scattering mean free
path is given by Eq. (71) or Eq. (77). The expression to use depends on the
choice of the “reference” homogeneous medium around which the permittivity
fluctuates: the former is thus associated to the Bruggeman mixing rule and the
latter to the Maxwell-Garnett mixing rule. The common ground that links these
two approaches, as discussed for instance by Mackay and Lakhtakia (2015), is
often overlooked.

For assemblies of small identical and non-absorbing particles, the scattering
mean free path is instead given by Eq. (99). A comparison with Eq. (77)
unveils a fundamental concept, namely that a fluid of tiny identical particles
with a fluctuating density (e.g., a fluid of molecules) can be assimilated to
a continuous medium with a fluctuating permittivity (or polarizability).
Indeed, taking the limit $qa\rightarrow 0$ in Eq. (99), one finds the
equivalence

$\displaystyle\left\langle\Delta\rho({\bf r})\Delta\rho({\bf
r}^{\prime})\right\rangle\alpha_{\text{p}}^{2}\equiv\left\langle\Delta\alpha({\bf
r})\Delta\alpha({\bf r}^{\prime})\right\rangle/v^{2}.$ (109)

This is the reason why Rayleigh’s and Einstein’s quantitative results on the
mean free paths (explaining the blue color of the sky) were essentially
identical Rayleigh (1899); Einstein (1910).

Finally, for assemblies of large, possibly resonant particles, the scattering
and transport mean free paths are respectively given by Eqs. (107) and (108).
These expressions were derived from a theory for the average intensity in the
large-scale approximation. Their strength lies in the fact that the
contributions from the individual particles and from structural correlations
are treated separately, thereby providing considerable physical insight.

##  III Structural properties of correlated disordered media

Correlated disordered media can exhibit a rich variety of complex
morphologies, which will impact light scattering and transport in many
different ways. This section is concerned with the statistical description of
the structural properties of these materials. For a more complete and thorough
description, we recommend the textbook by Torquato (2013). After comparing the
quantities describing spatial correlations and derived in the previous section
on realistic systems [Sec. III.1] and introducing a fundamental relation
between fluctuations of particle number and spatial correlations in point
patterns [Sec. III.2], we define the main classes of correlated disordered
media according to their pair correlation function $g_{2}(r)$ and structure
factor $S(q)$ [Sec. III.3]. We then review the main techniques to construct
numerically [Sec. III.4] and fabricate experimentally [Sec. III.5] such
complex materials in practice. The section is concluded with a brief summary
of the experimental techniques to characterize structural correlations in real
materials [Sec. III.6].

###  III.1 Continuous permittivity versus particulate models in practice

As shown in the previous section, light scattering in correlated disordered
media can be described either by a continuous permittivity model that relies
on a spatial permittivity correlation function or by a particulate model that
relies on a two-point correlation function in the specific case of localized
permittivity variations. To start this section, it is interesting to compare
these two pictures in practical cases.

We consider a classical and very relevant example for photonics, that is a 2D
assembly of impenetrable disks (diameter $a$) at two different packing
fractions $p=0.10$ and $0.50$, see Fig. 3. The disk packings were generated
with a compression algorithm Skoge  _et al._ (2006) that will be briefly
described in Sec. III.4. The top-left panel of Fig. 3 shows the permittivity
correlation function
$g_{\epsilon}(|\mathbf{r}-\mathbf{r}^{\prime}|)=h_{\epsilon}(|\mathbf{r}-\mathbf{r}^{\prime}|)+1$
for the disk packings. The function first displays a rapid decrease of
correlation followed by regular, vanishing oscillations. The very short-range
correlation here is mostly associated to the finite size of the disks and the
oscillations, which are stronger for higher packing fractions, are a signature
of correlations between neighboring particles.

Figure 3:  Description of structural correlations for hard disk (diameter $a$)
packings at packing fractions $p=10\%$ (black curves) and $p=50\%$ (blue
curves). (Left) Correlation function $g_{\epsilon}(r)=h_{\epsilon}(r)+1$
describing correlations due to the finite-size of the disks and to their
positional correlation, for two different topologies: (top) a packing of disks
and (bottom) a continuous network. (Right) Pair correlation function
$g_{2}(r)=h_{2}(r)+1$, which describe only the positional correlation between
the disk centers.

The generated disk positions can also serve as a basis to generate more
complex structures. An example often encountered in photonics is based on
Delaunay tesselations (described in Sec. III.4). In the bottom-left panel of
Fig. 3, we show the permittivity correlation function for these “inverted”
structures, as generated from the disk packing above. Strikingly, the behavior
of the correlation functions are very similar to those of the disk packing.
The major difference is observed at small distances due to a very different
morphology, but the curves become indistinguishable for
$|\mathbf{r}-\mathbf{r}^{\prime}|\gtrsim a$.

Finally, we can consider the pair-correlation function
$g_{2}(|\mathbf{r}-\mathbf{r}^{\prime}|)$. The results are shown in the right
panel of Fig. 3. The pair-correlation function for the disk packing is zero
for $|\mathbf{r}_{a}-\mathbf{r}_{b}|$ between $0$ to $2R$ due to the
impenetrability of the disks, and exhibits strong oscillations indicating
structural correlations in the relative position between particles. Note that
the amplitude of the oscillations is much larger than for
$g_{\epsilon}(|\mathbf{r}-\mathbf{r}^{\prime}|)$.

This simple comparison quite importantly shows that the use of the two-point
permittivity correlation is hardly sufficient to distinguish very different
disordered media. In fact, the permittivity correlation of the inverted disk
structure at $p=0.50$ would be strictly identical to that of the direct
structure by definition, while light scattering would evidently be markedly
different. This is a strong indication that scattering is both dramatically
affected by the local morphology of the system, which yields optical
resonances, and by structural correlations in the relative position between
scattering elements. Although the definition of a “scattering element” is
questionable for inverted structures such as connected networks, we will focus
on the particulate description of scattering media in the remainder of this
section.

###  III.2 Fluctuation-correlation relation

The description of point patterns underlying the structure of correlated
disordered media is central, and many descriptors may be used in general.

An important attribute of point patterns is the variance of the number $N$ of
points contained within windows $\Omega_{i}$ with volume $V$ and shifted by
$\mathbf{r}$. This quantity has a long history, several derivations are found
for both continuous and discrete disorder models Ornstein and Zernike (1914);
de Boer (1949); Van Kranendonk and Sipe (1977); Landau and Lifshitz (1980);
Martin and Yalcin (1980); Torquato and Stillinger (2003). Probabilistic
calculations eventually lead to a closed-form expression for the variance of
$N$

$\displaystyle\frac{\langle N^{2}(\mathbf{r},\Omega)\rangle-\langle
N(\mathbf{r},\Omega)\rangle^{2}}{\langle
N(\mathbf{r},\Omega)\rangle}=1+\rho\int_{\Omega_{1}\cap\Omega_{2}}h_{2}(\mathbf{r})d\mathbf{r},$
(110)

where $h_{2}(\mathbf{r})=g_{2}(\mathbf{r})-1$ is again the total correlation
function and the integral is performed over the intersection volume between
the two shifted windows Torquato and Stillinger (2003). Equation (110)
requires $h_{2}(\mathbf{r})$ to be integrable on the intersection volume and
go to 0 at large distances $r$. Interestingly, in the limit of large windows,
one finds that Torquato and Stillinger (2003)

$\displaystyle\lim_{V\rightarrow\infty}\frac{\langle
N^{2}(\mathbf{r},\Omega)\rangle-\langle
N(\mathbf{r},\Omega)\rangle^{2}}{\langle
N(\mathbf{r},\Omega)\rangle}=\lim_{|\mathbf{q}|\rightarrow 0}S(\mathbf{q}),$
(111)

which corresponds to the simplified definition given in Appendix C.2.

Equations (110) and (111) are remarkable in that they describe the spatial
fluctuations in the number of points in the pattern from its pair correlation
between points – or equivalently, its structure factor near 0, which is a
measurable quantity (e.g., by small-angle scattering, see Sec. III.6). A
Poisson point pattern $p_{N}=\langle N\rangle^{N}\exp\left[-\langle
N\rangle\right]/N!$ yields $\langle N^{2}\rangle=\langle N\rangle+\langle
N\rangle^{2}$, which, as expected, corresponds to a fully uncorrelated system
with $h_{2}(\mathbf{r})=0$ or $\lim_{|\mathbf{q}|\rightarrow
0}S(\mathbf{q})=1$. Implementing structural correlations at constant density
$\rho$ therefore results into weaker or stronger point density fluctuations.
Negative correlations are obtained when
$\rho\int_{\Omega_{1}\cap\Omega_{2}}h_{2}(\mathbf{r})d\mathbf{r}<0$, leading
to a sub-Poissonian fluctuations, while positive correlations are obtained
when $\rho\int_{\Omega_{1}\cap\Omega_{2}}h_{2}(\mathbf{r})d\mathbf{r}>0$,
leading to a super-Poissonian fluctuations. In the literature, such structures
are sometimes denoted as negatively- and positively-correlated, respectively
Davis and Mineev (2008). As illustrated in Fig. 4, they correspond to
situations in which the points either repel or attract themselves. As we will
see in Sec. IV, the impact of negative and positive correlations on optical
transport are markedly different.

Figure 4:  Illustration of the fluctuation-correlation relation with three
point patterns: (left) A negatively-correlated disordered medium, where points
tend to repel themselves; (middle) A Poisson point pattern, that is
uncorrelated; (right) A positively-correlated disordered medium, wherein
clustering is present. Two windows denoted as $\Omega_{1}$ and $\Omega_{2}$ in
which a number of points is counted are shifted by $\mathbf{r}$.

###  III.3 Classes of correlated disordered media

Figure 5 summarizes the most important classes of correlated disordered media
and their properties, that we will now describe specifically.

Figure 5:  Classes of correlated disordered media. (From top to bottom)
Illustration of a correlated disordered medium; Pair correlation function;
Structure factor; SEM image of a fabricated correlated disordered structure.
(From left to right) Disordered short-range correlated structures, SEM image
Adapted with permission from García  _et al._ (2008); Polycrystalline
structure, SEM image Adapted with permission from Salvarezza  _et al._ (1996);
Imperfect ordered structures, SEM image Adapted with permission from Garcia
_et al._ (2012); Disordered hyperuniform structures, SEM image Adapted with
permission from Haberko and Scheffold (2013); Disordered hierarchical
structures, SEM image Adapted with permission from Burresi  _et al._ (2012).

####  III.3.1 Short-range correlated disordered structures

Consider a volume containing a disordered ensemble of mobile, impenetrable
particles – i.e. a fluid of hard particles – at a low density. With increasing
particle density, the particles tend to organize themselves to fill space. In
this regime of low to moderate densities, the system exhibits no structural
correlation in the long range yet the impenetrability of the particles impose
a short-range correlation that increases with the packing fraction Hansen and
McDonald (1990). As shown in Fig. 5 (1st column), short-range structural
correlations give rise to decaying oscillations in the pair correlation
function $g_{2}$. The most likely distance to find a neighboring particle is
given by the position of the first peak and the decay of the higher-order
peaks, which is generally rapid, allows defining a correlation length. In
reciprocal space, such oscillations are also observed. When increasing short-
range correlations, the structure factor goes from a flat response around 1 to
sharper peaks whose amplitudes decrease with increasing $q$. Short-range
structural correlations can be described, for instance, by analytical or semi-
analytical solutions of the Ornstein-Zernike equation using the so-called
Percus-Yevick approximation Percus and Yevick (1958); Wertheim (1963) in three
dimensions and the Baus-Colot approximation in two dimensions Baus and Colot
(1987), respectively. At higher densities, these models become less accurate,
although, for slightly polydisperse systems, errors appear to cancel, and
predictions by Percus-Yevick approximation can describe experimental data up
to random close packing or jamming Frenkel  _et al._ (1986); Scheffold and
Mason (2009).

Short-range correlated disordered systems constitute the primary class found
in colloidal systems with isotropic interactions, since short-range
correlations stem from the impenetrability of particles in suspension.
Interestingly, densely packed sphere assemblies are also hyperuniform, but not
stealthy Zachary  _et al._ (2011), for details see section III.3.4. The
behavior of other repulsive particles, such as charge-stabilized particles,
can often be mapped onto the isotropic hard-sphere case Pusey and Van Megen
(1986); Gast and Russel (1998). The recent advent of colloids interacting via
sticky patches could open a pathway towards more complex structures through
self-assembly He  _et al._ (2020).

In general, short-range structural correlations are not limited to sphere
assemblies but can also be encoded in connected networks Florescu  _et al._
(2009); Liew  _et al._ (2011); Muller  _et al._ (2013), in which case the
individual scattering centers are more difficult to identify. This form of
correlated structure is widespread in natural photonic structures such as bird
feathers Saranathan  _et al._ (2012) and very popular in artificial photonic
structures fabricated by top-down techniques Liew  _et al._ (2011); Riboli
_et al._ (2014). Dry foams are promising candidates for correlated network
structures that can be made by self-assembly Ricouvier  _et al._ (2019).

####  III.3.2 Polycrystalline structures

For a disordered ensemble of identical hard particles, one reaches a liquid-
crystal coexistence at about 49% and a purely crystalline phase for
concentrations above 54.5% Pusey and Van Megen (1986); Pusey (1991); Gast and
Russel (1998); Zhu  _et al._ (1997). The equilibrium structure appears to be
face-centered-cubic but hexagonal close packed structures are also observed
and are found to be at least meta-stable Pusey  _et al._ (1989). This liquid
to crystal phase transition, also known as the Kirkwood-Alder transition Gast
and Russel (1998), is purely driven by the higher entropy of the crystalline
phase compared to the liquid phase. The densest packing of monodisperse
spheres in three dimensions is approximately 74%, also referred to as the
close-packing of equal spheres. Monodisperse particles usually assemble in
finite-sized crystal clusters. These clusters are randomly arranged and form a
polycrystalline materials Astratov  _et al._ (2002); Yang  _et al._ (2010b),
see Fig. 5 (2nd column). In the bulk, crystallites are formed by homogeneous
nucleation throughout the sample Pusey  _et al._ (1989). The size of the
crystal clusters is then typically several tens of $\mu$m, much larger than
the particle diameter $\sim\lambda$ but smaller than the usual sample size,
which is typically in the millimeter to centimeter range. The radially-
averaged pair correlation function exhibits peaks indicating the position of
the $n$th-order neighboring particles, as well as minima approaching zero.
Positional correlations vanish for distances exceeding the size of the crystal
clusters. Similarly, the structure factor shows well-defined Debye-Scherrer
rings due to Bragg scattering from randomly oriented crystal planes, that can
be identified in light scattering Pusey  _et al._ (1989), similarly to powder
diffraction in X-ray crystallography.

The formation of clusters of regular arrays in fluids of hard particles is
strongly influenced by the polydispersity of the particles, since particles of
very different sizes do not naturally arrange in a crystal. Indeed, for hard
sphere fluids with a polydispersity larger than 6-12%, crystallization is
avoided in three dimensions Pusey (1987). The spheres remain disordered and
particles enter a solid glass phase at about 58%. The glass can be further
compressed until the spheres ’jam’ forming what is known as a “randomly closed
packed” or “maximally jammed structure” in the literature Torquato  _et al._
(2000). The presence of some hidden structural order, crystalline precursors
or locally favoured structures in the glass and jammed phase is still being
discussed Zhang  _et al._ (2016).

Due to the unavoidable finite polydispersity, experimental realizations of
crystalline photonic structures based on colloidal suspensions are the
exception rather than the rule even at high packing fractions. By careful
synthesis of colloidal particles made from polystyrene or silica (SiO2), it is
however possible to induce crystallization rather easily Salvarezza  _et al._
(1996). These materials are usually polycrystalline and display some defects
and stacking faults to a varying degree. Polycrystalline structures are also
observed in natural photonic structures, such as opals. Interestingly,
relatively little is known about the comparison of scattering and light
transport between random-close-packed assemblies of spheres and
polycrystalline materials Yang  _et al._ (2010b), in particular when the size
of crystallites is gradually reduced to smaller length scales.

####  III.3.3 Imperfect ordered structures

The two previous classes of correlated disorder were obtained by “adding
order” into a fully-disordered (uncorrelated) system. Materials with
correlated disorder can also be obtained starting from the other limit, that
is a periodic system with random perturbations, see Fig. 5 (3rd column). In
systems of infinite size, both the pair-correlation function $g_{2}$ and the
structure factor $S$ are characterized by a series of Dirac peaks located at
$\mathbf{r}-\mathbf{r}^{\prime}=u_{1}\mathbf{a}_{1}+u_{2}\mathbf{a}_{2}+u_{3}\mathbf{a}_{3}$
with $u_{i}\in\mathbb{Z}$ and $\mathbf{a}_{i}$ the lattice vectors, and
$\mathbf{G}=v_{1}\mathbf{b}_{1}+v_{2}\mathbf{b}_{2}+v_{1}\mathbf{b}_{2}$ with
$v_{i}\in\mathbb{Z}$ and $\mathbf{b}_{i}$ the reciprocal lattice vectors,
respectively Kittel (1976); Joannopoulos  _et al._ (2011). If the position of
a point of the lattice is randomly shifted (e.g., with normal distribution)
around its nominal position, this results in a broadening of the Dirac peaks
with a width that depends on the disorder amplitude. By contrast with the
previous classes of disordered systems, disorder in such a periodic-on-average
structure does not impact the correlation length, which remains infinite. The
striking consequence of this is that the structure factor remains
characterized by Dirac peaks of vanishing width (for systems of infinite size)
yet with an amplitude that decays with $q$. This reasoning does not only hold
for periodic media but also for quasiperiodic systems (e.g. Fibonacci
structures, Penrose tiles, etc.), which also exhibit a discrete pattern in
reciprocal space - meaning structural order and long-range correlation. Real
systems are however never exactly periodic, due to fabrication imperfections
and in practice this also leads to a finite correlation length Meseguer  _et
al._ (2002); López (2003); Koenderink  _et al._ (2005); Nelson  _et al._
(2011).

A plethora of studies of ordered photonic crystal structures with
imperfections, both numerical and experimental, can be found in the literature
Soukoulis (2012). Defects were also added intentionally, either at random or
selected positions, to study the interplay between defect states, density of
states, wave tunneling and percolation, random diffuse scattering, and
directed Bragg scattering of light Florescu  _et al._ (2010); García  _et al._
(2009); Fernandes  _et al._ (2013); Aeby  _et al._ (2020). Moreover, the
interaction between the band structures and defect scattering is fascinating,
since it might lead to other critical coherent transport phenomena such as
Anderson localization of light John (1987). Defect states can also be
introduced in a photonic crystal to deliberately implement a particular
function, such as optical sensing applications, lasing, or optical circuitry
Joannopoulos  _et al._ (2011); Soukoulis (2012); Nelson  _et al._ (2011).

####  III.3.4 Disordered hyperuniform structures

One of the important characteristics of point patterns is how the number of
points contained in a given volume fluctuate with various disorder
realizations Torquato (2013). This quantity is related to the notion of
spatial uniformity. For a Poisson point process, one shows from Eq. (110) that
the variance in the number of points $N$ contained in a $d$-dimensional sphere
of radius $R$ grows as the sphere volume, i.e. $\langle N^{2}\rangle-\langle
N\rangle^{2}=\langle N\rangle\sim R^{d}$. This result holds for many
disordered point patterns. By contrast, the same analysis performed on a
periodic pattern shows that the variance grows with the surface of the sphere,
$\langle N^{2}\rangle-\langle N\rangle^{2}\sim R^{d-1}$. In a founding work,
Torquato and Stillinger (2003) proposed to define a general class of point
patterns, dubbed “hyperuniform”, the property of which is to exhibit point
number fluctuations scaling as the surface of the window, that is slower than
expected for usual disordered media. Hyperuniformity encompasses periodic,
quasi-periodic but also – very interestingly in the framework of this review –
a subclass of disordered systems, see Fig. 5 (4th column) for an illustration
and Torquato (2018) for a recent review. It was observed numerically that
maximally jammed packings of spheres and platonic solids tend to a
hyperuniform structure Jiao and Torquato (2011). While such long-range
fluctuations can hardly be observed on the pair-correlation function,
hyperuniform point patterns can be recognized from the behavior of the
structure factor at low values

$\lim_{\mathbf{q}\rightarrow 0}S(\mathbf{q})=0.$ (112)

Of particular interest in photonics are so-called “stealthy” hyperuniform
structures, for which $S(\mathbf{q})=0$ for $0<q\leqslant q_{\text{max}}$,
where $q_{\text{max}}$ may be set to an arbitrary value. The region of zero
structure factor is often followed by oscillations similar to those found in
short-range disordered correlated media Froufe-Pérez  _et al._ (2016).

The concept of hyperuniformity in photonics has first been introduced in a
numerical study by Florescu  _et al._ (2009). Important efforts have been put
since then on the fabrication of hyperuniform disordered systems, which could
be achieved so far by lithography in 2D Man  _et al._ (2013) and 3D Muller
_et al._ (2013), block copolymer assembly Zito  _et al._ (2015), emulsion
routes Weijs  _et al._ (2015); Ricouvier  _et al._ (2017); Piechulla  _et al._
(2018), and spinodal solid-state dewetting Salvalaglio  _et al._ (2020).

####  III.3.5 Disordered fractal structures

In all classes of disordered point patterns discussed above, the average
number of points $N$ contained in a $d$-dimensional sphere of radius $R$ is
expected to grow as $\langle N\rangle\propto R^{d}$ \- by doubling the
observation radius for a 3D point pattern, the number of points increases by a
factor $2^{3}=8$. This scaling is however not a general rule. Introduced by
Mandelbrot (1967), the concept of fractals encompasses systems for which the
power-law scaling of the mass with the system size does not have the Euclidean
dimension as an exponent. More specifically, for fractal point patterns, we
have $\langle N\rangle\propto R^{d_{\text{f}}}$, where $d_{\text{f}}$ is a
non-integer fractal dimension. Fractality has a dramatic impact on the
structure, as illustrated in Fig. 5 (5th column). First, it is statistically
self-similar, meaning that the structure is statistically identical whatever
the scale on which it is looked at (though lower and upper bounds are always
met in practice). Second, it exhibits enormous local density fluctuations and
high lacunarity Allain and Cloitre (1991), meaning that both very dense and
very empty regions are found. As a result, the pair correlation function can
be shown to decay as a power-law as $g_{2}(R)\sim R^{-\alpha}$ and similarly
for the structure factor, $S(q)-1\sim|q|^{-(d-\alpha)}$ assuming that
$0<\alpha<d$. Depending on the process of structure formation, one can
directly relate the exponent $\alpha$ with the fractal dimension
$d_{\text{f}}$. For instance, clustering described by the Soneira-Peebles
model gives $\alpha=d-d_{\text{f}}$, leading to
$S(q)-1\sim|q|^{-d_{\text{f}}}$ Soneira and Peebles (1977). It is important to
note that real systems generally exhibit lower and upper bounds in their
fractal nature. This implies that the power-law decays are observed on finite
range ($g_{2}(R)$ eventually goes to 1 at large $R$).

Fractal disordered optical materials are encountered in a wide variety of
colloidal aggregates that form naturally for certain charged particles Meakin
(1987) as well as in certain emulsions Bibette  _et al._ (1993). They can also
be designed in a laboratory by inserting spacing particles with a size
distribution that covers several orders of magnitude in a statistically-
homogeneous disordered medium Barthelemy  _et al._ (2008); Bertolotti  _et
al._ (2010a).

###  III.4 Numerical construction of correlated disordered media

Numerical simulations of the complex heterogeneous morphologies play a key
role in colloidal chemistry and soft matter physics. For light scattering
studies, modelled structured materials are taken as input data for solving
Maxwell’s equations. Here, we present some standard numerical approaches that
have been used in the literature to generate correlated disordered structures.

Random packings of hard spheres in different dimensions is of great interest
due to their structural and thermodynamic properties. The numerical generation
of such ensembles plays a key role in research, especially in the case of
random close packings (RCP) Song  _et al._ (2008). When the packing fraction
of the system is kept below a few tens of percent, a random sequential
absorption model (RSA) model Widom (1966) is suitable to generate large (non-
equilibrium) ensembles in any dimensionality. In the RSA model, new points are
randomly added to the system following a uniform distribution. The new point
is rejected if closer than a given distance to any of the previous points in
the pattern. A careful management of the coordinates storage in appropriate
structures lead to very efficient algorithms. However, the former algorithms
become dramatically inefficient when the rejection rate is high, the maximum
filling fraction being $\phi\simeq 54\%$ in 2D Wang (2000) and $\phi\simeq
38\%$ in 3D Meakin and Jullien (1992).

In order to achieve larger packing fractions up to the jamming packing (at a
filling fraction $\phi\simeq 64\%$ in 3D) several approaches leading to
efficient algorithms have been developed. The Lubachevsky-Stillinger (or
compression) algorithm Lubachevsky and Stillinger (1990) generates random
packings of any physically-realistic $\phi$ by placing a set of $N$ particles
of vanishingly small size in a closed or periodic domain at random, and then
letting the particles grow at a given rate, move and collide (elastically),
until the desired $\phi$ is reached. This algorithm has been successfully used
in the study of sphere packings in any dimensionality Skoge  _et al._ (2006)
and is largely used in photonics Conley  _et al._ (2014); Froufe-Pérez  _et
al._ (2016). An approach named “ideal amorphous solids” was also proposed by
Lee  _et al._ (2010) to realize maximally random jammed packings of
polydisperse particles. The method relies on building aggregates of touching
spheres by placing spheres one by one around a center of mass.

Enlarging the kind of correlations encountered in sphere packings requires the
use of interaction potentials beyond the hard sphere model. Simple two-body
interaction potentials such as Lennard-Jones together with standard Monte-
Carlo techniques have been used to generate assemblies of scatterers in
different phases De Sousa  _et al._ (2016). Two and three-body interaction
models like the Stillinger-Weber model Stillinger and Weber (1985) can be used
to generate fully connected dielectric networks showing the same statistical
structural properties (coordination and angle statistics) as amorphous silicon
or diamond. Stealthy hyperuniform point patterns have been generated
numerically using a suitable pairwise, long-range potential in the real space
Froufe-Pérez  _et al._ (2016). Quite often, the structures are generated using
molecular dynamics. Being a very vast and mature field, various softwares are
nowadays available, including NAMD Phillips  _et al._ (2005) and CHARMMS
Brooks  _et al._ (2009), which are both extensively used in the field of
chemistry and biochemistry. Other software packages include HOOMD-blue
Anderson  _et al._ (2008), which is implemented for GPU, and LAMMPS LAMMPS
(2019), which exploits massive parallelization Plimpton (1995).

Instead of using constraints in real space as in the case of hard spheres,
targeted interaction potentials can be obtained by imposing constraints on the
structure factor in reciprocal space Uche  _et al._ (2004), which allowed
realizing, e.g., stealthy hyperuniform structures Batten  _et al._ (2008);