11institutetext: Kevin Vynck 22institutetext: Université Claude Bernard Lyon
1, CNRS, iLM, F-69622 Villeurbanne, France, 22email<EMAIL_ADDRESS>33institutetext: Armel Pitelet 44institutetext: Université de Bordeaux,
Institut d’Optique Graduate School, CNRS, LP2N, F-33400 Talence, France,
44email<EMAIL_ADDRESS>55institutetext: Louis Bellando
66institutetext: Université de Bordeaux, CNRS, LOMA, F-33400 Talence, France,
66email<EMAIL_ADDRESS>77institutetext: Philippe Lalanne
88institutetext: Université de Bordeaux, Institut d’Optique Graduate School,
CNRS, LP2N, F-33400 Talence, France, 88email:
<EMAIL_ADDRESS>
# Specular reflection and transmission of electromagnetic waves by disordered
metasurfaces

Kevin Vynck  Armel Pitelet  Louis Bellando and Philippe Lalanne

###### Abstract

Planar, disordered assemblies of small particles incorporated in layered media
– sometimes called “disordered metasurfaces” in the recent literature – are
becoming widespread in optics and photonics. Their ability to scatter light
with exotic angular and spectral features in reflection and transmission, as
well as their suitability to scalable fabrication techniques, makes them
promising candidates for certain applications, ranging from thin-film
photovoltaics to visual appearance design. This chapter introduces the basic
concepts and theoretical models for the specular (a.k.a. coherent) reflectance
and transmittance of electromagnetic waves by disordered metasurfaces. After
describing the classical scattering formalism for discrete media, we establish
known analytical expressions for the reflection and transmission coefficients
of disordered particle monolayers on layered substrates. Two classical models,
based on the independent scattering approximation (ISA) and the effective
field approximation (EFA), are presented. Their accuracy is examined by
comparing predictions with those obtained from rigorous full-wave computations
using an in-house multiple-scattering code. This chapter may serve as a
starting point to students and researchers who wish to dive into the topic and
explore the potential of disordered metasurfaces for applications.

##  1 Introduction

The phenomenon of wave scattering by random rough surfaces is encountered in
many physical settings and techniques. These include radar remote sensing and
imaging of planetary surfaces, optical characterization of surfaces and LIDAR,
underwater acoustics, radio communications, and seismology, to cite only a
few. A typical problem of interest is a coherent wave impinging on a surface
with random heterogeneities (e.g., random height variations, supported
particles at random positions, etc.). The intensity of the reflected signal
generally takes the form of a speckle, showing bright and dark spots due to
constructive and destructive interferences between waves coming from different
positions of the sample. As might be expected, the statistical properties of
the measured signal contain statistical information on the surface topology.
Considerable advances have been made in the past decades on establishing this
relationship, as discussed in several excellent textbooks
beckmann1987scattering ; stover1995optical ; bedeaux2004optical ;
vesperinas2006scattering ; voronovich2013wave .

This period has also been marked by the exploration of novel strategies to
control the emission, propagation and confinement of light with nanostructures
– a discipline nowadays known as “nanophotonics” benisty2022introduction .
Research has been greatly stimulated by advances in nanofabrication
techniques, which enabled the realization of a large variety of dielectric and
metallic subwavelength structures with nanometer-scale resolution, in one, two
and three dimensions lopez2003materials ; mascolo2011production . One major
outcome has been to unveil the potential of finely-engineered high-index
resonant nano-objects to exhibit exotic properties, such as a light bending
effect mirin2009light , a strong scattering anisotropy gomez2011electric  or a
side-dependent coupling to waveguides wu2020intrinsic . By optimizing the
composition, size and shape of the nano-objects and placing them on a
substrate, one creates a so-called “metasurface” that can deviate, focus, or
alter the polarization state of an incident beam. The topic has become
extremely popular and is covered by many review and expert opinion articles
kuznetsov2016optically ; lalanne2017metalenses ; kamali2018review ;
shaltout2019spatiotemporal .

While in most situations, the nano-objects are ordered, i.e., arranged
deterministically on a lattice, and the surfaces are realized by top-down
fabrication techniques (e.g., lithography), increasing efforts are being made
to realize metasurfaces by bottom-up techniques, like colloidal chemistry and
self-assembly. A major motivation is the expected benefit in terms of
scalability and manufacturing cost. The resulting structures are usually
disordered, even though some degree of control over structural correlations is
possible roach2022controlling . Light interaction with disordered metasurfaces
is receiving growing attention since a few years, with a focus on
functionalities that do not require a precise control over the position of the
inclusions, such as light absorption for photodetection, sensing and
photovoltaics moreau2012controlled ; chevalier2015absorbing ;
stewart2017toward ; piechulla2021antireflective , light extraction for organic
LEDs jouanin2016designer ; donie2021planarized , and light scattering for
augmented reality displays bertin2018correlated  and visual appearance design
vynck2022visual .

Disordered metasurfaces may thus be seen as a special type of random
scattering surface and, naturally, one is led to the question of whether (and
how well) previously-established theoretical models can be used for the
purpose. The present chapter aims at answering, at least in part, to this
question.

In modelling studies, the intensity scattered by a surface is generally
decomposed into two components: a so-called coherent component, which
corresponds to the specularly reflected and transmitted waves111The term
“specular” is used here independently for both reflected and transmitted
waves, but one may also use the term “ballistic” for the transmitted waves.,
and a so-called incoherent component, which corresponds to the diffuse (non-
specular) intensity. Formally, as will be shown below, the former is given by
the ensemble-averaged electric field produced by the system and leads to the
definition of complex reflection and transmission coefficients, whereas the
latter is given by the fluctuations of the electric field around its average
value and leads to the definition of an angle-resolved scattering diagram. The
theoretical modelling of the diffuse intensity created by planar, disordered
assemblies of resonant particles in layered media remains, to our knowledge,
an open problem, as models are in general limited to particle monolayers
embedded in a uniform background loiko2018incoherent  or to particles much
smaller than the wavelength sasihithlu2016surface . We will therefore focus
here on the specular response of disordered metasurfaces, for which a more
exhaustive literature exists.

The problem of determining the reflection and transmission coefficients of
particle monolayers is closely related to that of the electromagnetic
homogenization. For particles that are very small compared to the wavelength
in the embedding medium (i.e., in the quasi-static regime), the system may be
seen as a homogeneous layer with an effective permittivity tensor222The
effective permittivity of a monolayer of particles is not a scalar in
general., which may be obtained from the Maxwell-Garnett mixing rule
sasihithlu2016surface , the Yamaguchi model yamaguchi1974optical  or its
extensions fedotov2003optical ; toudert2012advanced . An alternative approach,
still applicable to monolayers of very small particles only, proposed by
Bedeaux and Vlieger bedeaux1973phenomenological ; vlieger1980statistical  and
implemented in the software GranFilm lazzari2002granfilm , relies on the
notions of excess fields and surface susceptibilities to solve the
electromagnetic problem. However, for larger particles possibly exhibiting
high-order resonances, one needs to resort to the electromagnetic scattering
theory by discrete media tsang2004scattering , an approach that has notably
been followed by García-Valenzuela and colleagues pena2006coherent ;
garcia2012multiple .

In light of the increasing attention given to disordered metasurfaces in
recent years, we propose in this chapter to give some keys to understand the
physical origin and underlying approximations, and test the validity of known
analytical expressions for the reflection and transmission coefficients of
particle monolayers in layered media. The chapter is decomposed as follows. In
Section 2, starting from the simple case of an individual, finite-size
heterogeneity in a uniform medium [Fig. 1(a)], we introduce the basic
theoretical concepts in electromagnetic scattering, namely the wave
propagation equation and the transition operator, the dyadic Green function
and its angular spectrum representation, and the scattering amplitude.
Eventually, we arrive at an expression of the field scattered by a particle in
a planewave basis. In Section 3, we treat the problem of the coherent
intensity scattered by a monolayer of particles in a layered geometry. After
introducing the multiple-scattering equations and formally defining the
coherent and incoherent intensities, we present two models for the specular
reflection and transmission by (infinite) monolayers of particles in a uniform
medium [Fig. 1(b)]. We conclude this section with a brief explanation of how
the results can be generalized to handle particle monolayers above a layered
substrate [Fig. 1(c)]. In Section 4, finally, we test the validity of these
two analytical models on different systems, made of either metallic or
dielectric particles, on a bare semi-infinite or layered substrate, by
comparing the analytical predictions with those from rigorous, full-wave
computations using an in-house multiple-scattering code.

Figure 1: Illustration of the different problems of interest in this chapter
and definition of certain variables. A discrete medium composed of particles
is illuminated by a planewave ${\bf E}_{\text{b}}$ with amplitude
$E_{\text{b}}$, wavevector ${\bf k}_{\text{i}}$ and polarization
$\hat{\bm{e}}_{\text{i}}$. (a) Electromagnetic scattering by an isolated
particle in a uniform medium. In Section 2, we will derive an expression of
the scattered field in a planewave basis. (b) Specular reflection and
transmission by a planar, disordered assembly of identical particles in a
uniform medium. Section 3 aims mainly at presenting two models for the complex
reflection and transmission coefficients of particle monolayers in uniform
media, $r_{\text{coh}}$ and $t_{\text{coh}}$, respectively. (c) Specular
reflection and transmission by a disordered metasurface made of identical
particles on a layered substrate. We will explain, at the end of Section 3,
how the models can be used to predict the specular reflection, described by
the coefficient $r_{\text{st}}$, of particle monolayers on layered substrates.
In Section 4, we will test those models using full-wave multiple-scattering
computations.

##  2 Basics of electromagnetic scattering by particles

###  2.1 Wave equations

We consider a finite region of space filled with a non-magnetic material
(relative permeability, $\mu(\mathbf{r})=1$), described by a relative
permittivity ${\bm{\epsilon}}({\bf r})$ in a uniform host medium with relative
scalar permittivity ${\epsilon_{\text{b}}}$. The permittivity variation
$\delta{\bm{\epsilon}}({\bf r})\equiv({\bm{\epsilon}}({\bf
r})-{\epsilon_{\text{b}}}{\bf I})$, with ${\bf I}$ the unit tensor, defines a
compact heterogeneity, which will be our particle lateron, though no
assumption on the particle composition, size or shape is in fact necessary. We
consider harmonic fields at frequency $\omega$ with the $e^{-i\omega t}$
convention and drop hereafter the explicit dependence in the permittivities,
fields, etc. for simplicity.

We start from the macroscopic Maxwell’s (curl) equations for the electric and
magnetic fields, ${\bf E}$ and ${\bf H}$, at frequency $\omega$ with a current
density source $\mathbf{J}$,

$\displaystyle\bm{\nabla}\times{\bf E}({\bf r})$ $\displaystyle=$
$\displaystyle i\omega\mu_{0}{\bf H}({\bf r}),$ (1)
$\displaystyle\bm{\nabla}\times{\bf H}({\bf r})$ $\displaystyle=$
$\displaystyle-i\omega\epsilon_{0}\bm{\epsilon}({\bf r}){\bf E}({\bf
r})+\mathbf{J}({\bf r}).$ (2)

Taking the curl of Eq. (1) and inserting Eq. (2) in the resulting expression
leads to a vector wave propagation equation for the electric field

$\displaystyle\bm{\nabla}\times\bm{\nabla}\times{\bf E}({\bf
r})-\frac{\omega^{2}}{c^{2}}\bm{\epsilon}({\bf r}){\bf E}({\bf
r})=i\omega\mu_{0}\mathbf{J}({\bf r}),$ (3)

where $c^{2}=(\epsilon_{0}\mu_{0})^{-1}$, $c$ being the speed of
electromagnetic waves in vacuum.

In a scattering problem, it is convenient to decompose the field as the sum of
a background field and a scattered field, ${\bf E}={\bf E}_{\text{b}}+{\bf
E}_{\text{s}}$. The background field ${\bf E}_{\text{b}}$ is the solution of
the wave propagation equation with the source term but without the
heterogeneity,

$\displaystyle\bm{\nabla}\times\bm{\nabla}\times{\bf E}_{\text{b}}({\bf
r})-k_{\text{b}}^{2}{\bf E}_{\text{b}}({\bf r})=i\omega\mu_{0}\mathbf{J}({\bf
r}),$ (4)

with $k_{\text{b}}^{2}=k_{0}^{2}{\epsilon_{\text{b}}}$ and $k_{0}=\omega/c$.
From Eqs. (3) and (4), one therefore reaches a wave equation for the scattered
field only,

$\displaystyle\bm{\nabla}\times\bm{\nabla}\times{\bf E}_{\text{s}}({\bf
r})-k_{\text{b}}^{2}{\bf E}_{\text{s}}({\bf
r})=k_{0}^{2}\delta{\bm{\epsilon}}({\bf r}){\bf E}({\bf r}).$ (5)

The scattered field in the background medium is thus generated by the total
field in the heterogeneity volume, which is nothing but the polarization
density ${\bf P}({\bf r})/\epsilon_{0}=\delta{\bm{\epsilon}}({\bf r}){\bf
E}({\bf r})$.

###  2.2 Lippmann-Schwinger equation

The electromagnetic problem described by Eq. (5) can conveniently be
reformulated by introducing the notion of dyadic Green function tai1994dyadic
. The dyadic Green function ${\bf G}({\bf r},{\bf r}^{\prime})$ of a medium
describes the electric field produced at point ${\bf r}$ by a radiating point
electric dipole at point ${\bf r}^{\prime}$. While it can be defined for an
arbitrary environment, we are interested here in the dyadic Green function in
the background medium, ${\bf G}_{\text{b}}$, which is then the solution of

$\displaystyle\bm{\nabla}\times\bm{\nabla}\times{\bf G}_{\text{b}}({\bf
r},{\bf r}^{\prime})-k_{\text{b}}^{2}{\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime})={\bf I}\delta({\bf r}-{\bf r}^{\prime}).$ (6)

Multiplying both sides by $k_{0}^{2}\delta{\bm{\epsilon}}({\bf
r}^{\prime}){\bf E}({\bf r}^{\prime})$, integrating over ${\bf r}^{\prime}$
and using Eq. (5) leads to

$\displaystyle{\bf E}_{\text{s}}({\bf r})=k_{0}^{2}\int{\bf G}_{\text{b}}({\bf
r},{\bf r}^{\prime})\delta{\bm{\epsilon}}({\bf r}^{\prime}){\bf E}({\bf
r}^{\prime})d{\bf r}^{\prime}.$ (7)

Only the field in the heterogeneity volume contributes to the scattered field.
Thus, ones reaches the so-called Lippmann-Schwinger equation for the total
field,

$\displaystyle{\bf E}({\bf r})={\bf E}_{\text{b}}({\bf r})+k_{0}^{2}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime})\delta{\bm{\epsilon}}({\bf
r}^{\prime}){\bf E}({\bf r}^{\prime})d{\bf r}^{\prime}.$ (8)

The equation is physically very insightful, but not very practical as such
because the field at ${\bf r}$ depends explicitly on the field at ${\bf
r}^{\prime}$.

###  2.3 Transition operator

To reach more practical expressions, we can express the scattered field in Eq.
(7) by successive iterations as a series of scattering events,

$\displaystyle{\bf E}_{\text{s}}({\bf r})$ $\displaystyle=$ $\displaystyle
k_{0}^{2}\int{\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime})\delta{\bm{\epsilon}}({\bf r}^{\prime}){\bf E}_{\text{b}}({\bf
r}^{\prime})d{\bf r}^{\prime}$ (9) $\displaystyle+$ $\displaystyle
k_{0}^{4}\int{\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime})\delta{\bm{\epsilon}}({\bf r}^{\prime}){\bf G}_{\text{b}}({\bf
r}^{\prime},{\bf r}^{\prime\prime})\delta{\bm{\epsilon}}({\bf
r}^{\prime\prime}){\bf E}_{\text{b}}({\bf r}^{\prime\prime})d{\bf
r}^{\prime}d{\bf r}^{\prime\prime}$ $\displaystyle+$ $\displaystyle...$

The first integral term gives the scattered field due to one interaction
within the heterogeneity, the second to two interactions within the
heterogeneity, etc. The scattered field is now expressed in terms of the
background field only. Eventually, all multiple-scattering orders can be
incorporated into a transition operator ${\bf T}$, hereafter called
T-operator333${\bf T}$ is called here the T-operator to avoid the confusion
with the T-matrix, which, in the numerical method of the same name
mishchenko1996t , is generally defined in a basis of vector spherical wave
functions., such that

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\int{\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime}){\bf T}({\bf r}^{\prime},{\bf r}^{\prime\prime}){\bf
E}_{\text{b}}({\bf r}^{\prime\prime})d{\bf r}^{\prime}d{\bf
r}^{\prime\prime}.$ (10)

with

$\displaystyle{\bf T}({\bf r},{\bf
r}^{\prime})=k_{0}^{2}\delta{\bm{\epsilon}}({\bf r})\left[\delta({\bf r}-{\bf
r}^{\prime})+\int{\bf G}_{\text{b}}({\bf r},{\bf r}^{\prime\prime}){\bf
T}({\bf r}^{\prime\prime},{\bf r}^{\prime})d{\bf r}^{\prime\prime}\right],$
(11)

Equation (10) can be physically interpreted as follows: a background field
${\bf E}_{\text{b}}$ arriving a point ${\bf r}^{\prime\prime}$ in the
heterogeneity induces, via the T-operator ${\bf T}$, a dipole moment at
position ${\bf r}^{\prime}$, which then radiates in the background medium
towards point ${\bf r}$ via the dyadic Green function ${\bf G}_{\text{b}}$.
Note that ${\bf T}$ is spatially non-local in general, as it contains all
multiple-scattering events relating two points within the heterogeneity. It is
a complicated, yet intrinsic quantity: once known, it allows predicting the
scattered field for any background field.

Because the microscopic details of the heterogeneity are “hidden” in the
T-operator, we will hereafter talk about particles explicitly when referring
to an individual, finite-size heterogeneity.

The T-operator can be obtained analytically using Mie theory for simple
objects (e.g., particles with spherical symmetry) and numerically by solving
Maxwell’s equations otherwise bohren2008absorption . Interestingly, for
assemblies of particles, one can define a T-operator of the whole system,
which can be expressed in terms of the T-operators of the individual
particles. The computation can be achieved by solving the multiple-scattering
problem using, for instance, so-called the T-matrix method mishchenko1996t .

###  2.4 Dyadic Green function in a homogeneous medium

To proceed further, we need to define more properly the dyadic Green function
above. We will skip the derivation, which can be found elsewhere tai1994dyadic
; jackson1999classical ; novotny2012principles ; carminati2021principles . In
translationally-invariant and isotropic media, the dyadic Green function ${\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime})\equiv{\bf G}_{\text{b}}({\bf r}-{\bf
r}^{\prime})$ in Cartesian coordinates is given by

$\displaystyle{\bf G}_{\text{b}}({\bf r}-{\bf r}^{\prime})=-\frac{\delta({\bf
R})}{3k_{\text{b}}^{2}}{\bf
I}+\text{PV}\left\\{\frac{\exp[ik_{\text{b}}R]}{4\pi R}\left[{\bf I}-{\bf
u}\otimes{\bf u}+\frac{ik_{\text{b}}R-1}{(k_{\text{b}}R)^{2}}\left({\bf
I}-3{\bf u}\otimes{\bf u}\right)\right]\right\\},$ (12)

where $\delta(.)$ is the Dirac delta function, PV stands for Principal Value,
$\otimes$ denotes the tensor product, and we defined ${\bf R}={\bf r}-{\bf
r}^{\prime}=R{\bf u}$ for conciseness. The first term corresponds to the
singularity of the dyadic Green function at the origin and the second term to
the non-local part, which itself contains evanescent and propagating
components. In the far-field, $k_{\text{b}}R\gg 1$, only propagating waves
remain and the expression reduces to

$\displaystyle{\bf G}_{\text{b}}({\bf r}-{\bf
r}^{\prime})\sim\frac{\exp[ik_{\text{b}}R]}{4\pi R}\left[{\bf I}-{\bf
u}\otimes{\bf u}\right].$ (13)

Equation (13) shows that the field radiated by a compact heterogeneity behaves
at large distances as a transverse, outgoing spherical wave, where the
transverse nature is due to the term between brackets.

The representation in outgoing spherical waves is well suited to scattering
problems in volumes but when the problem involves planar geometries (e.g.,
particle monolayers, interfaces), it is more convenient to express the waves
as a linear combination of planewaves. This is the so-called angular spectrum
representation, a.k.a. Weyl expansion or Weyl identity weyl1919ausbreitung ,
and that we will now apply to the dyadic Green function.

For this, let us start by expressing the dyadic Green function in real space
in terms of its Fourier transform,

$\displaystyle{\bf G}_{\text{b}}({\bf r}-{\bf
r}^{\prime})=\frac{1}{(2\pi)^{3}}\int\frac{1}{k^{2}-k_{\text{b}}^{2}}\left[{\bf
I}-\frac{{\bf k}\otimes{\bf k}}{k_{\text{b}}^{2}}\right]\exp\left[i{\bf
k}\cdot({\bf r}-{\bf r}^{\prime})\right]d{\bf k}.$ (14)

Defining ${\bf r}=\left[{\bf r}_{\parallel},z\right]$, ${\bf k}=\left[{\bf
k}_{\parallel},k_{z}\right]$, $k_{\parallel}=|{\bf k}_{\parallel}|$, and
$\gamma=\sqrt{k_{\text{b}}^{2}-k_{\parallel}^{2}}$ with $\text{Re}[\gamma]>0$
and $\text{Im}[\gamma]>0$, Eq. (14) can then be rewritten as

$\displaystyle{\bf G}_{\text{b}}({\bf r}-{\bf
r}^{\prime})=\frac{1}{(2\pi)^{2}}\int{\bf g}_{\text{b}}({\bf
k}_{\parallel},z,z^{\prime})\exp\left[i{\bf k}_{\parallel}\cdot({\bf
r}_{\parallel}-{\bf r}^{\prime}_{\parallel})\right]d{\bf k}_{\parallel},$ (15)

with

$\displaystyle{\bf g}_{\text{b}}({\bf
k}_{\parallel},z,z^{\prime})=\frac{1}{2\pi}\int\frac{1}{k_{z}^{2}-\gamma^{2}}\left[{\bf
I}-\frac{{\bf k}\otimes{\bf
k}}{k_{\text{b}}^{2}}\right]\exp\left[ik_{z}(z-z^{\prime})\right]dk_{z}.$ (16)

The solution of the integral over $k_{z}$ in Eq. (16) is evaluated taking care
of the poles by using Cauchy’s residue theorem to the contour integrals,
leading to

$\displaystyle{\bf g}_{\text{b}}({\bf
k}_{\parallel},z,z^{\prime})=-\frac{\delta(z-z^{\prime})}{k_{\text{b}}^{2}}\bm{\delta}_{zz}+\text{PV}_{zz}\left\\{\frac{i}{2}\frac{1}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]\exp\left[i\gamma|z-z^{\prime}|\right]\right\\}$
(17)

where $\bm{\delta}_{ij}$ is the Kronecker delta tensor and
${\bm{\kappa}}^{\pm}=[{\bf k}_{\parallel},\pm\gamma]$ with
$\pm=\text{sign}(z-z^{\prime})$. Inserting this expression into Eq. (15) and
keeping the non-singular part only (i.e., $z\neq z^{\prime}$), we obtain the
following angular spectrum represenation of the dyadic Green function,

$\displaystyle{\bf G}_{\text{b}}({\bf r}-{\bf
r}^{\prime})=\frac{i}{2(2\pi)^{2}}\int\frac{1}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]\exp\left[i{\bm{\kappa}}^{\pm}\cdot({\bf
r}-{\bf r}^{\prime})\right]d{\bf k}_{\parallel}.$ (18)

Note that the integral is made over the entire (infinite) reciprocal plane
described by ${\bf k}_{\parallel}$. Because
$\gamma=\sqrt{k_{\text{b}}^{2}-k_{\parallel}^{2}}$, wavevectors fulfilling
$k_{\parallel}^{2}\leq k_{\text{b}}^{2}$ or
$k_{\parallel}^{2}>k_{\text{b}}^{2}$ correspond to propagating or evanescent
waves, respectively.

###  2.5 Scattering of a planewave by a particle

To complete this section, we will now express the scattered field ${\bf
E}_{\text{s}}$ produced by an individual particle in the angular spectrum
representation. Let us first consider that the particle is centered at
position ${\bf r}_{j}$ and is described by a T-operator ${\bf T}_{j}$. The
relative position of an individual particle in space does not affect its
scattering properties, of course, but the relative position between particles
in an assembly does. This step is therefore important for our purpose, as will
be shown lateron. Equation (10) then reads

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\int{\bf G}_{\text{b}}({\bf r},{\bf
r}^{\prime}){\bf T}_{j}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{b}}({\bf r}^{\prime\prime})d{\bf
r}^{\prime}d{\bf r}^{\prime\prime}.$ (19)

The T-operator can be expressed in terms of its Fourier transform as

$\displaystyle{\bf T}_{j}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{6}}\int\exp\left[i{\bf p}^{\prime}\cdot({\bf
r}^{\prime}-{\bf r}_{j})\right]$ (20) $\displaystyle\times$ $\displaystyle{\bf
T}_{j}({\bf p}^{\prime},{\bf p}^{\prime\prime})\exp\left[-i{\bf
p}^{\prime\prime}\cdot({\bf r}^{\prime\prime}-{\bf r}_{j})\right]d{\bf
p}^{\prime}d{\bf p}^{\prime\prime}.$

Let us then consider that the background field ${\bf E}_{\text{b}}$ is a
planewave with amplitude $E_{\text{b}}$, wavevector ${\bf
k}_{\text{i}}=k_{\text{b}}\hat{\bf k}_{\text{i}}$ and polarization
$\hat{\bm{e}}_{\text{i}}$,

$\displaystyle{\bf E}_{\text{b}}({\bf
r})=E_{\text{b}}\hat{\bm{e}}_{\text{i}}\exp\left[i{\bf k}_{\text{i}}\cdot{\bf
r}\right].$ (21)

Using Eqs. (20) and (18), we can rewrite Eq. (19) as

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\frac{iE_{\text{b}}}{2(2\pi)^{8}}$
$\displaystyle\int$ $\displaystyle\frac{1}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]{\bf
T}_{j}({\bf p}^{\prime},{\bf
p}^{\prime\prime})\hat{\bm{e}}_{\text{i}}\exp\left[i{\bm{\kappa}}^{\pm}\cdot({\bf
r}-{\bf r}^{\prime})\right]$ (22) $\displaystyle\times$
$\displaystyle\exp\left[i{\bf p}^{\prime}\cdot({\bf r}^{\prime}-{\bf
r}_{j})\right]\exp\left[-i{\bf p}^{\prime\prime}\cdot({\bf
r}^{\prime\prime}-{\bf r}_{j})\right]\exp\left[i{\bf k}_{\text{i}}\cdot{\bf
r}^{\prime\prime}\right]$ $\displaystyle\times$ $\displaystyle d{\bf
r}^{\prime}d{\bf r}^{\prime\prime}d{\bf p}^{\prime}d{\bf
p}^{\prime\prime}d{\bf k}_{\parallel}.$

Having that $\int\exp[i\mathbf{q}{\bf r}]d{\bf
r}=(2\pi)^{3}\delta(\mathbf{q})$, the integrals over ${\bf r}^{\prime}$ and
${\bf r}^{\prime\prime}$, and then ${\bf p}^{\prime}$ and ${\bf
p}^{\prime\prime}$, lead to ${\bf p}^{\prime}={\bm{\kappa}}^{\pm}$ and ${\bf
p}^{\prime\prime}={\bf k}_{\text{i}}$. After simplification, we thus obtain

$\displaystyle{\bf E}_{\text{s}}({\bf
r})=\frac{iE_{\text{b}}}{2(2\pi)^{2}}\exp\left[i{\bf k}_{\text{i}}\cdot{\bf
r}_{j}\right]$ $\displaystyle\int$ $\displaystyle\frac{1}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]{\bf
T}_{j}({\bm{\kappa}}^{\pm},{\bf k}_{\text{i}})\hat{\bm{e}}_{\text{i}}$ (23)
$\displaystyle\times$ $\displaystyle\exp\left[i{\bm{\kappa}}^{\pm}\cdot({\bf
r}-{\bf r}_{j})\right]d{\bf k}_{\parallel}.$

Similarly to Eq. (15), we can define ${\bf E}_{\text{s}}({\bf
k}_{\parallel},z,z_{j})$ such that

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\frac{1}{(2\pi)^{2}}\int{\bf
E}_{\text{s}}({\bf k}_{\parallel},z,z_{j})\exp\left[i{\bf
k}_{\parallel}\cdot({\bf r}_{\parallel}-{\bf r}_{j,\parallel})\right]d{\bf
k}_{\parallel},$ (24)

leading to

$\displaystyle{\bf E}_{\text{s}}({\bf k}_{\parallel},z,z_{j})$
$\displaystyle=$
$\displaystyle\frac{iE_{\text{b}}}{2}\frac{\exp[i\gamma|z-z_{j}|]}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]{\bf
T}_{j}({\bm{\kappa}}^{\pm},{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\exp[i{\bf k}_{\text{i}}\cdot{\bf
r}_{j}].$

Note that the wavevector ${\bf k}_{\parallel}$ still covers the entire
reciprocal space, so that both evanescent and propagating waves are
considered. As we are more interested here into waves propagating in the far
field, let us select a propagating wave with wavevector ${\bf
k}_{\text{s}}=k_{\text{b}}\hat{\bf k}_{\text{s}}$ and consider a plane above
the particle, $z>z_{j}$. We therefore use ${\bm{\kappa}}^{+}=[{\bf
k}_{\text{s},\parallel},+k_{\text{s},z}]$. For waves propagating in the far
field, we can further introduce the vector scattering amplitude, a classical
quantity in scattering theory carminati2021principles 444Important quantities,
notably the extinction cross-section and differential scattering cross-
section, can be calculated directly from the scattering amplitude.,

$\displaystyle{\bf f}_{j}(\hat{\bf k}_{\text{s}},\hat{\bf
k}_{\text{i}})=\frac{1}{4\pi}\left[{\bf I}-\hat{\bf
k}_{\text{s}}\otimes\hat{\bf k}_{\text{s}}\right]{\bf
T}_{j}(k_{\text{b}}\hat{\bf k}_{\text{s}},k_{\text{b}}\hat{\bf
k}_{\text{i}}).$ (26)

This leads to

$\displaystyle{\bf E}_{\text{s}}({\bf k}_{\text{s},\parallel},z,z_{j})=2i\pi
E_{\text{b}}\frac{\exp[ik_{\text{s},z}(z-z_{j})]}{k_{\text{s},z}}{\bf
f}_{j}(\hat{\bf k}_{\text{s}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\exp[i{\bf k}_{\text{i}}\cdot{\bf
r}_{j}],$ (27)

where $k_{\text{s},z}=k_{\text{b}}\cos(\theta_{\text{i}})$, with
$\theta_{\text{i}}$ the angle of incidence with respect to the the $z$-axis.
Equation (27) formally expresses the planewave decomposition of the field
scattered by a particle at ${\bf r}_{j}$ in the far field in terms of the
particle scattering amplitude. Similar steps will later be made to derive
expressions for the field scattered by a monolayer of particles and the
resulting reflection and transmission coefficients.

##  3 Specular reflection and transmission by particle monolayers

###  3.1 Multiple scattering by discrete media

In the previous section, we considered the electromagnetic scattering problem
for an individual particle, centered at point ${\bf r}_{j}$, described by a
T-operator ${\bf T}_{j}$, and illuminated by a background field ${\bf
E}_{\text{b}}$. In the present section, we will consider scattering by an
ensemble of such particles. In the general case, each particle may have a
different composition, size and shape, and therefore a different T-operator.
Following Eq. (19), the field scattered by a set of $N$ particles is now
simply the sum of the field scattered by each of them,

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\sum_{j=1}^{N}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bf T}_{j}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{exc}}^{j}({\bf
r}^{\prime\prime})d{\bf r}^{\prime}d{\bf r}^{\prime\prime}.$ (28)

We defined here ${\bf E}_{\text{exc}}^{j}$ as the field exciting the $j$-th
particle, which, very importantly, can differ from the background field ${\bf
E}_{\text{b}}$. Indeed, quite intuitively, this exciting field for the $j$-th
particle should be the sum of the background field and the field scattered by
all other particles, $l\neq j$, as

$\displaystyle{\bf E}_{\text{exc}}^{j}({\bf r})={\bf E}_{\text{b}}({\bf
r})+\sum_{l\neq j}^{N}\int{\bf G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bf
T}_{l}({\bf r}^{\prime}-{\bf r}_{l},{\bf r}^{\prime\prime}-{\bf r}_{l}){\bf
E}_{\text{exc}}^{l}({\bf r}^{\prime\prime})d{\bf r}^{\prime}d{\bf
r}^{\prime\prime}.$ (29)

Equations (28) and (29) define the multiple scattering problem. The exciting
field ${\bf E}_{\text{exc}}^{p}$ depends on the positions and properties of
all other particles and is generally unknown. This set of equations can be
solved numerically on finite ensembles of particles using T-matrix methods
mishchenko1996t , as mentioned above, however, the scattering response of a
specific configuration of particles is generally of limited interest. Instead,
one looks for statistical properties of the scattering response, such as the
average field or intensity, intensity fluctuations, field-field or intensity-
intensity correlations, etc carminati2021principles . Theoretically, the
problem is tackled in general by taking the configurational average of the
relevant quantity (the field, intensity, etc.) from Eqs. (28) and (29), and
making certain approximations on the interaction between particles.

###  3.2 Coherent and incoherent intensity

The intensity scattered by a medium, which is the quantity that is usually
measured experimentally, can formally be decomposed into two components,
denoted as coherent and incoherent. To show this, let us write the scattered
field as the sum of its average value and a fluctuating part,

$\displaystyle{\bf E}_{\text{s}}({\bf r})=\langle{\bf E}_{\text{s}}({\bf
r})\rangle+\delta{\bf E}_{\text{s}}({\bf r}),\qquad\text{with
}\left\langle\delta{\bf E}_{\text{s}}({\bf r})\right\rangle=0.$ (30)

The configurational average (formally defined below) is written here with
angle brackets, $\langle\cdot\rangle$. Calculating the configurational average
of the electric field norm squared, $\left\langle|{\bf E}({\bf
r})|^{2}\right\rangle$, which is directly proportional to the average
intensity, shows that it can be decomposed into two terms,

$\displaystyle\langle|{\bf E}|^{2}\rangle$ $\displaystyle=$
$\displaystyle|{\bf E}_{\text{b}}|^{2}+2\text{Re}[{\bf
E}_{\text{b}}\cdot\langle{\bf E}_{\text{s}}^{*}\rangle]+|\langle{\bf
E}_{\text{s}}\rangle|^{2}+\langle|\delta{\bf E}_{\text{s}}|^{2}\rangle,$ (31)
$\displaystyle=$ $\displaystyle|\langle{\bf E}\rangle|^{2}+\langle|\delta{\bf
E}_{\text{s}}|^{2}\rangle.$

The first term, $|\langle{\bf E}\rangle|^{2}=|{\bf
E}_{\text{b}}|^{2}+2\text{Re}[{\bf E}_{\text{b}}\cdot\langle{\bf
E}_{\text{s}}^{*}\rangle]+|\langle{\bf E}_{\text{s}}\rangle|^{2}$, is known as
the coherent intensity555$|\langle{\bf E}\rangle|^{2}$ contains an
interference term between the background and scattered fields, hence the term
“coherent”.. In volume scattering, this component defines the extinction
coefficient of the medium, which describes the attenuation rate of an incident
wave due to scattering and absorption. In surface scattering, this leads to
specularly reflected and transmitted waves defined by their reflection and
transmission coefficients.

The second term, $\langle|\delta{\bf E}_{\text{s}}|^{2}\rangle=\langle|{\bf
E}|^{2}\rangle-|\langle{\bf E}\rangle|^{2}=\langle|{\bf
E}_{\text{s}}|^{2}\rangle-|\langle{\bf E}_{\text{s}}\rangle|^{2}$, is known as
the incoherent intensity666The term “incoherent” may be misleading since the
diffuse intensity can be impacted by interference between scattered waves, as
in the case of correlated disorder vynck2021light . It corresponds to the
diffuse intensity, created by field fluctuations from realization to
realization and characterized by a differential scattering cross-section.

As stated previously, we are interested here in the specular reflection and
transmission of a monolayer of particles, and will therefore focus on theory
for the average field. For further details on the theory for the intensity, we
recommend the textbooks from Tsang and Kong tsang2004scattering  and from
Carminati and Schotland carminati2021principles .

###  3.3 Average scattered field

We start by defining the configurational average for discrete media. Our main
ingredient here is the $N$-dimensional probability density function $p({\bf
R})$ of finding $N$ particles in the configuration ${\bf R}=[{\bf r}_{1},{\bf
r}_{2},\cdots,{\bf r}_{N}]$. The configurational average of a variable
$\bm{f}({\bf r})$ is defined as

$\displaystyle\langle\bm{f}({\bf r})\rangle$ $\displaystyle=$
$\displaystyle\int\bm{f}({\bf r},{\bf R})p({\bf R})d{\bf R},$ (32)

with $\bm{f}({\bf r},{\bf R})$ the variable evaluated at ${\bf r}$ for a
specific configuration ${\bf R}$.

Let us then proceed by calculating the configurational average of the
scattered field, Eq. (28), considering a set of $N$ identical particles
described by a unique T-operator, ${\bf T}_{j}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j})\equiv{\bf T}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j})$. We thus have

$\displaystyle\langle{\bf E}_{\text{s}}({\bf r})\rangle=\sum_{j=1}^{N}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime})\left\langle{\bf T}({\bf
r}^{\prime}-{\bf r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j}){\bf
E}_{\text{exc}}^{j}({\bf r}^{\prime\prime})\right\rangle d{\bf
r}^{\prime}d{\bf r}^{\prime\prime},$ (33)

with

$\displaystyle\left\langle{\bf T}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{exc}}^{j}({\bf
r}^{\prime\prime})\right\rangle=\int{\bf T}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{exc}}^{j}({\bf
r}^{\prime\prime},{\bf R})p({\bf R})d{\bf R},$

following Eq. (32). Using $p({\bf R})=p({\bf R}|{\bf r}_{j})p({\bf r}_{j})$
with $p({\bf R}|{\bf r}_{j})$ the conditional probability density function for
${\bf R}$ having ${\bf r}_{j}$ fixed, and $d{\bf R}=d{\bf r}_{1}d{\bf
r}_{2}\cdots d{\bf r}_{N}$, we obtain

$\displaystyle\left\langle{\bf T}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{exc}}^{j}({\bf
r}^{\prime\prime})\right\rangle=\int{\bf T}({\bf r}^{\prime}-{\bf r}_{j},{\bf
r}^{\prime\prime}-{\bf r}_{j})\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}({\bf
r}^{\prime\prime},{\bf r}_{j})p({\bf r}_{j})d{\bf r}_{j},$ (35)

where $\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}$ is the average exciting
field ${\bf E}_{\text{exc}}^{j}$ having the particle $j$ fixed at ${\bf
r}_{j}$, given by

$\displaystyle\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}({\bf
r}^{\prime\prime},{\bf r}_{j})=\int{\bf E}_{\text{exc}}^{j}({\bf
r}^{\prime\prime},{\bf R})p({\bf R}|{\bf r}_{j})d{\bf r}_{1}\cdots d{\bf
r}_{j-1}d{\bf r}_{j+1}\cdots d{\bf r}_{N}.$ (36)

Inserting Eq. (35) into Eq. (33), we finally obtain

$\displaystyle\langle{\bf E}_{\text{s}}({\bf r})\rangle=\sum_{j=1}^{N}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bf T}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j})\langle{\bf
E}_{\text{exc}}^{j}\rangle_{j}({\bf r}^{\prime\prime},{\bf r}_{j})p({\bf
r}_{j})d{\bf r}_{j}d{\bf r}^{\prime}d{\bf r}^{\prime\prime}.$ (37)

The difficulty at this stage is to express $\langle{\bf
E}_{\text{exc}}^{j}\rangle_{j}$ in a suitable manner as to solve Eq. (37).
Without any approximation, $\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}$ can be
expressed, in a similar manner as done above using Eq. (29), in terms of the
average exciting field with two particles $j$ and $l$ fixed, and so on.
Different levels of approximations are obtained depending on the order at
which the truncation is made tsang2004scattering . Hence, the lowest-order
approximation, known as independent scattering approximation (ISA), completely
neglects the interaction between particles, assuming then $\langle{\bf
E}_{\text{exc}}^{j}\rangle_{j}({\bf r}^{\prime\prime},{\bf r}_{j})\simeq{\bf
E}_{\text{b}}({\bf r}^{\prime\prime})$. Instead, the so-called effective field
(or Foldy’s) approximation (EFA) considers that $\langle{\bf
E}_{\text{exc}}^{j}\rangle_{j}({\bf r}^{\prime\prime},{\bf
r}_{j})\simeq\langle{\bf E}\rangle({\bf r}^{\prime\prime})$, whereas the so-
called quasi-crystalline approximation (QCA) assumes that $\langle{\bf
E}_{\text{exc}}^{j}\rangle_{jl}({\bf r}^{\prime\prime},{\bf r}_{j},{\bf
r}_{l})\simeq\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}({\bf
r}^{\prime\prime},{\bf r}_{j})$. We will treat below the ISA and EFA models.
The QCA was used by García-Valenzuela et al. garcia2012multiple  and
constitutes, to our knowledge, the state-of-the-art on multiple-scattering
models for the coherent intensity from particle monolayers.

###  3.4 Independent scattering approximation (ISA)

The ISA assumes, as the name indicates, that the particles behave
independently from each other; an electromagnetic wave interacts only once
with each particle, the mutual interactions between particles are therefore
neglected. The exciting field [Eq. (29)] is then given by the background field
only. Evidently, this approximation should only be valid for very dilute
systems. Assuming then

$\displaystyle\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}({\bf
r}^{\prime\prime},{\bf r}_{j})\simeq{\bf E}_{\text{b}}({\bf
r}^{\prime\prime}),$ (38)

Eq. (37) becomes

$\displaystyle\langle{\bf E}_{\text{s}}({\bf r})\rangle=\sum_{j=1}^{N}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bf T}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j}){\bf E}_{\text{b}}({\bf
r}^{\prime\prime})p({\bf r}_{j})d{\bf r}_{j}d{\bf r}^{\prime}d{\bf
r}^{\prime\prime}.$ (39)

Using Eqs. (18), (20) and (21), and having $\int\exp[i({\bf
p}^{\prime}-{\bm{\kappa}}^{\pm})\cdot{\bf r}^{\prime}]d{\bf
r}^{\prime}=(2\pi)^{3}\delta({\bf p}^{\prime}-{\bm{\kappa}}^{\pm})$ and
$\int\exp[i({\bf k}_{\text{i}}-{\bf p}^{\prime\prime})\cdot{\bf
r}^{\prime\prime}]d{\bf r}^{\prime\prime}=(2\pi)^{3}\delta({\bf
k}_{\text{i}}-{\bf p}^{\prime\prime})$ leads to

$\displaystyle\langle{\bf E}_{\text{s}}({\bf
r})\rangle=\frac{iE_{\text{b}}}{2(2\pi)^{2}}$ $\displaystyle\int$
$\displaystyle\frac{1}{\gamma}\left[{\bf
I}-\frac{{\bm{\kappa}}^{\pm}\otimes{\bm{\kappa}}^{\pm}}{k_{\text{b}}^{2}}\right]{\bf
T}({\bm{\kappa}}^{\pm},{\bf k}_{\text{i}})\exp[i{\bm{\kappa}}^{\pm}\cdot{\bf
r}]\hat{\bm{e}}_{\text{i}}$ (40) $\displaystyle\times$
$\displaystyle\left\\{\sum_{j=1}^{N}\int\exp[i({\bf
k}_{\text{i}}-{\bm{\kappa}}^{\pm})\cdot{\bf r}_{j}]p({\bf r}_{j})d{\bf
r}_{j}\right\\}d{\bf k}_{\parallel}$

Considering that the particles are randomly distributed on a surface of area
$S$ at $z=0$ (for simplicity) leads to $p({\bf r}_{j})=\delta(z_{j})/S$. In
addition, since all particles obey the same statistics, the sum reduces to a
prefactor $N$. Then, we consider the limit of an infinite surface,
$\lim_{N,S\rightarrow\infty}(N/S)=\rho$, with $\rho$ the surface density. In
this limit, the term between curled brackets in Eq. (40) becomes

$\displaystyle\lim_{N,S\rightarrow\infty}\left\\{\frac{N}{S}\int\exp[i({\bf
k}_{\text{i}}-{\bm{\kappa}}^{\pm})\cdot{\bf r}_{j}]\delta(z_{j})d{\bf
r}_{j}\right\\}=(2\pi)^{2}\rho\delta({\bf k}_{\text{i},\parallel}-{\bf
k}_{\parallel}),$ (41)

thereby imposing the conservation of the parallel wavevector between the
incident and scattered waves. Two solutions appear, ${\bm{\kappa}}^{\pm}=[{\bf
k}_{\text{i},\parallel},\pm k_{\text{i},z}]\equiv{\bf k}_{\text{t/r}}$, which
correspond to the specularly reflected and transmitted planewaves. Having
$k_{\text{i},z}=k_{\text{b}}\cos(\theta_{\text{i}})$ and using Eq. (26) for
the scattering amplitude, we reach

$\displaystyle\lim_{N,S\rightarrow\infty}\langle{\bf E}_{\text{s}}({\bf
r})\rangle$ $\displaystyle=$ $\displaystyle
E_{\text{b}}\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}{\bf
f}(\hat{{\bf k}}_{\text{t/r}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\exp\left[i{\bf k}_{\text{t/r}}\cdot{\bf
r}\right],$ (42) $\displaystyle\equiv$ $\displaystyle\langle{\bf
E}_{\text{s}}(\hat{{\bf k}}_{\text{t/r}})\rangle\exp\left[i{\bf
k}_{\text{t/r}}\cdot{\bf r}\right].$ (43)

Note the similarity with Eq. (27) obtained for an individual particle, which
makes sense considering the independent scattering approximation made here.
Indeed, independent scattered waves emerging from random (but statistically
uniform) positions on a surface add up coherently to form a planewave whose
amplitude depends only on the scattering properties of the individual particle
(i.e., the source of the scattered field) and the particle density.

Having now the average scattered field as a planewave with parallel wavevector
${\bf k}_{\text{i},\parallel}$ either in reflection or in transmission, we may
define the complex reflection and transmission coefficients, $r_{\text{coh}}$
and $t_{\text{coh}}$, of the particle monolayer, as

$\displaystyle r_{\text{coh}}E_{\text{b}}$ $\displaystyle=$
$\displaystyle\hat{\bm{e}}_{\text{r}}\cdot\left\langle{\bf
E}_{\text{s}}(\hat{{\bf k}}_{\text{r}})\right\rangle,$ (44) $\displaystyle
t_{\text{coh}}E_{\text{b}}$ $\displaystyle=$ $\displaystyle
E_{\text{b}}+\hat{\bm{e}}_{\text{t}}\cdot\left\langle{\bf
E}_{\text{s}}(\hat{{\bf k}}_{\text{t}})\right\rangle,$ (45)

with $\hat{\bm{e}}_{\text{r}}$ and $\hat{\bm{e}}_{\text{t}}$ the polarization
vectors of the specularly reflected and transmitted planewaves, respectively.
Noting that ${\bf k}_{\text{t}}={\bf k}_{\text{i}}$, we obtain

$\displaystyle r_{\text{coh}}$ $\displaystyle=$
$\displaystyle\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}\left(\hat{\bm{e}}_{\text{r}}\cdot{\bf
f}(\hat{{\bf k}}_{\text{r}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\right),$ (46) $\displaystyle
t_{\text{coh}}$ $\displaystyle=$ $\displaystyle
1+\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}\left(\hat{\bm{e}}_{\text{t}}\cdot{\bf
f}(\hat{\bf k}_{\text{i}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\right).$ (47)

For spherical particles, for which the scattering response depends on the
scattering angle rather than the incident and scattered angles, and producing
no polarization conversion
($\hat{\bm{e}}_{\text{t}}=\hat{\bm{e}}_{\text{i}}$), these expressions can
further be simplified using the following definitions

$\displaystyle{\bf f}(\hat{\bf k}_{\text{i}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}$ $\displaystyle\equiv$ $\displaystyle
f(0)\hat{\bm{e}}_{\text{i}},$ (48) $\displaystyle{\bf f}(\hat{{\bf
k}}_{\text{r}},\hat{{\bf k}}_{\text{r}})\hat{\bm{e}}_{\text{r}}$
$\displaystyle\equiv$ $\displaystyle f(0)\hat{\bm{e}}_{\text{r}},$ (49)
$\displaystyle{\bf f}(\hat{{\bf k}}_{\text{r}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}$ $\displaystyle\equiv$ $\displaystyle
f_{p}(\pi-2\theta_{\text{i}})\hat{\bm{e}}_{\text{r}},$ (50) $\displaystyle{\bf
f}(\hat{\bf k}_{\text{i}},\hat{{\bf k}}_{\text{r}})\hat{\bm{e}}_{\text{r}}$
$\displaystyle\equiv$ $\displaystyle
f_{p}(\pi-2\theta_{\text{i}})\hat{\bm{e}}_{\text{i}},$ (51)

where $f(0)$ is the forward scattering amplitude and
$f_{p}(\pi-2\theta_{\text{i}})$ is the polarization-dependent scattering
amplitude in the direction of specular reflection. Equations (46) and (47)
thus become

$\displaystyle r_{\text{coh}}$ $\displaystyle=$
$\displaystyle\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}f_{p}(\pi-2\theta_{\text{i}}),$
(52) $\displaystyle t_{\text{coh}}$ $\displaystyle=$ $\displaystyle
1+\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}f(0).$ (53)

We have therefore derived first equations for the specular reflection and
transmission coefficients by particle monolayers, assuming here that the
particles do not interact with each other. These expressions can be found in
the literature on metasurfaces. Let us emphasize, however, that the
coefficients diverge at grazing angles ($\theta_{\text{i}}$ approaching
$\pi/2$). Intuitively, indeed, the interaction between particles should become
more significant at oblique incidence, as the particles start to “shadow” each
other. As we will see in Section 4, it can nevertheless be a good
approximation for dilute media (1% surface coverage).

###  3.5 Effective field approximation (EFA)

The EFA is a higher-order approximation than the ISA, in the sense that the
average exciting field on particle $j$ (with particle $j$ fixed) is
approximated not only by the background field but also the average scattered
field, which partly accounts for the interaction between particles. Assuming
thus

$\displaystyle\langle{\bf E}_{\text{exc}}^{j}\rangle_{j}({\bf
r}^{\prime\prime},{\bf r}_{j})\simeq\langle{\bf E}\rangle({\bf
r}^{\prime\prime}),$ (54)

Eq. (37) reads

$\displaystyle\langle{\bf E}_{\text{s}}({\bf r})\rangle=\sum_{j=1}^{N}\int{\bf
G}_{\text{b}}({\bf r},{\bf r}^{\prime}){\bf T}({\bf r}^{\prime}-{\bf
r}_{j},{\bf r}^{\prime\prime}-{\bf r}_{j})\langle{\bf E}\rangle({\bf
r}^{\prime\prime})p({\bf r}_{j})d{\bf r}_{j}d{\bf r}^{\prime}d{\bf
r}^{\prime\prime}.$ (55)

Having $\langle{\bf E}\rangle={\bf E}_{\text{b}}+\langle{\bf
E}_{\text{s}}\rangle$, this brings us to a similar problem as that with the
Lippmann-Schwinger equation, Eq. (8), namely that the field at a point ${\bf
r}$ depends on the field at another point ${\bf r}^{\prime\prime}$. As
previously, this can be solved by successive iterations, thereby expressing
the problem as a series of scattering events between particles (not within a
single heterogeneity as in the derivation of the T-operator).

A similar procedure as done above for the ISA can be followed and applied to
the resulting infinite sum, leading to

$\displaystyle\lim_{N,S\rightarrow\infty}\langle{\bf E}_{\text{s}}({\bf
r})\rangle$ $\displaystyle=$ $\displaystyle
E_{\text{b}}\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}{\bf
f}(\hat{{\bf k}}_{\text{t/r}},\hat{\bf
k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\exp\left[i{\bf k}_{\text{t/r}}\cdot{\bf
r}\right]$ (56) $\displaystyle+$ $\displaystyle
E_{\text{b}}\rho^{2}\left(\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}\right)^{2}{\bf
f}(\hat{{\bf k}}_{\text{t/r}},\hat{{\bf k}}_{\text{t/r}}){\bf f}(\hat{{\bf
k}}_{\text{t/r}},\hat{\bf k}_{\text{i}})\hat{\bm{e}}_{\text{i}}\exp\left[i{\bf
k}_{\text{t/r}}\cdot{\bf r}\right]$ $\displaystyle+$ $\displaystyle\cdots$

The first integral term of the sum involves one vector scattering amplitude
and is identical to Eq. (42), the second integral term involves two vector
scattering amplitudes (note the different input and output wavevector
directions), etc. The sum takes the form of a geometric series, which, using
Eqs. (44) and (45), leads to new expressions for the reflection and
transmission coefficients

$\displaystyle r_{\text{coh}}$ $\displaystyle=$
$\displaystyle\hat{\bm{e}}_{\text{r}}\cdot\left[{\bf
I}-\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}{\bf f}(\hat{{\bf
k}}_{\text{r}},\hat{{\bf
k}}_{\text{r}})\right]^{-1}\left(\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}\right){\bf
f}(\hat{{\bf k}}_{\text{r}},\hat{\bf k}_{\text{i}})\hat{\bm{e}}_{\text{i}},$
(57) $\displaystyle t_{\text{coh}}$ $\displaystyle=$
$\displaystyle\hat{\bm{e}}_{\text{t}}\cdot\left[{\bf
I}-\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}{\bf f}(\hat{\bf
k}_{\text{i}},\hat{\bf k}_{\text{i}})\right]^{-1}\hat{\bm{e}}_{\text{i}},$
(58)

where $[\cdot]^{-1}$ denotes the dyadic inverse.

For spherical particles, we can use the equalities in Eqs. (48)-(51) to reach

$\displaystyle r_{\text{coh}}$ $\displaystyle=$
$\displaystyle\frac{\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}f_{p}(\pi-2\theta_{\text{i}})}{1-\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}f(0)},$
(59) $\displaystyle t_{\text{coh}}$ $\displaystyle=$
$\displaystyle\frac{1}{1-\rho\frac{2i\pi}{k_{\text{b}}\cos(\theta_{\text{i}})}f(0)},$
(60)

We have thus reached a second set of equations for the reflection and
transmission coefficients of particle monolayers, which considers the
interaction between particles in a mean-field sense. Importantly, unlike in
the ISA, the coefficients now behave correctly at grazing angles. Taking the
limit $\theta_{\text{i}}\rightarrow\pi/2$ leads to a reflection coefficient
$r_{\text{coh}}$ that approaches -1 (and thus, a reflectance
$|r_{\text{coh}}|^{2}$ approaching 1) since
$f_{p}(\pi-2\theta_{\text{i}})\rightarrow f(0)$, and a transmission
coefficient $t_{\text{coh}}$ approaching 0. We will see in Section 4 that the
EFA leads to physically sound and quantitatively accurate predictions for
moderately dense systems (typically, 10% surface coverage) even at relatively
large incident angles ($60^{\circ}$).

As a final remark, note that these equations have been derived heuristically
in Ref. pena2006coherent  by supposing that the exciting field is the average
field transmitted through the monolayer. Instead, we show here that this
result stems from a well-handled approximation that naturally appears in the
multiple-scattering expansion.

###  3.6 Generalization to particle monolayers on layered substrates

We now move to the last step of the derivation, which is to consider the
impact of a layered environment on the specular reflection and transmission of
a particle monolayer. A first possible approach would be to start from
Maxwell’s equations and solve the electromagnetic scattering problem using the
layered medium as a background. The dyadic Green function should then be the
one of the layered geometry paulus2000accurate . While this could be done for
very small particles behaving as electric dipoles sasihithlu2016surface , in
which case scattering by an individual particle is not described by a non-
local T-operator but by a local polarizability, it turns out to be much more
challenging for large particles.

A more practical approach is to consider that the interaction with the layered
geometry is mediated by the average scattered field from the particle
monolayer in the uniform background garcia2012multiple . In other words, the
particle monolayer is treated as planar interface with reflection and
transmission coefficients given by $r_{\text{coh}}$ and $t_{\text{coh}}$, and
the reflection and transmission coefficients of the entire structure are
determined by classical formulas of electromagnetic waves in layered media
yariv2005optical . In addition to the approximations made at the level of the
particle monolayer only, the mutual interaction between particles via the
interfaces of the environment is therefore discarded in this approach, which
may be problematic in certain situations, for instance, when incident
planewaves couple efficiently to guided (photonic or plasmonic) modes in the
layered medium.

In the framework of the latter approach and using classical recursive
relations for wave propagation in thin films yariv2005optical , the reflection
coefficient $r_{\text{st}}$ of a particle monolayer above a layered substrate
(made of isotropic materials, to disregard polarization
conversion777Multilayered substrates with anisotropic materials could be
considered as well, in which case the reflection and transmission coefficients
should be written as $2\times 2$ matrices for TE and TM (or s and p)
polarizations. Several publicly available codes can be used for the purpose
hugonin2021reticolo ; bay2022pyllama .) is simply given by

$\displaystyle
r_{\text{st}}=r_{\text{coh}}+\frac{r_{\text{sub}}t_{\text{coh}}^{2}\exp[2ik_{0}n_{\text{b}}a\cos(\theta_{\text{i}})]}{1-r_{\text{coh}}r_{\text{sub}}\exp[2ik_{0}n_{\text{b}}a\cos(\theta_{\text{i}})]},$
(61)

where $r_{\text{coh}}$ and $t_{\text{coh}}$ are the reflection and
transmission coefficients of the particle monolayer – determined with any
preferred model (ISA, EFA, QCA, …), $r_{\text{sub}}$ is the reflection
coefficient of the substrate alone,
$n_{\text{b}}=\sqrt{{\epsilon_{\text{b}}}}$ is the refractive index of the
medium in which the particles are embedded, and $a$ is the height separating
the monolayer and the first substrate interface ($a$ equals the particle
radius for spherical particles). Importantly, $r_{\text{sub}}$ already takes
into account all multiple reflections in-between the interfaces if several
layers are considered. This quantity can be calculated using, for instance,
the transfer matrix mackay2020transfer  or the scattering matrix method
ko1988scattering .

In the next section, we will test the ISA and EFA models for the specular
reflection of particle monolayers on layered substrates in several
experimentally-relevant situations.

##  4 Numerical validation of theoretical predictions

To verify the validity of the ISA and EFA models for the problem of interest,
we use an in-house multiple-scattering code developed by Jean-Paul Hugonin at
Laboratoire Charles Fabry (Institut d’Optique Graduate School, CNRS,
Université Paris Saclay). The numerical method belongs to the broad family of
T-matrix methods mishchenko1996t , that is probably the most adapted to solve
wave scattering problems by discrete media. T-matrix methods benefit from the
high degree of analyticity of the functions describing the wave propagation
between particles, typically dyadic Green functions or vector spherical wave
functions, depending on the problem of interest and the specific
implementation. With the recent emergence of disordered media in photonics
vynck2021light , the modelling of multiple scattering by particles in layered
media has gained considerable attention, leading to the development of new
methods bertrand2020global  and powerful publicly available tools
egel2021smuthi .

One of the major difficulties in numerical studies of multiple scattering in
disordered media comes from the fact that the phenomena of interest may take
place on mesoscopic scales, much larger than the wavelength, implying that
care should be taken to avoid finite-size effects. This issue is often
circumvented partially by simulating larger systems illuminated by (smaller)
collimated beams, but this strategy generally looses reliability for studies
at large incident angles.

In addition to being capable of modelling large particles incorporated in an
arbitrary layered medium, made of isotropic or anisotropic materials, the in-
house code used here implements the supercell method, that is an artificial
periodization of the electromagnetic problem, wherein each supercell can
contain a large number of particles [Fig. 2] to cope with finite-size effects.
Technical details on the implementation can be found in Ref.
langlais2014cooperative . The numerical code has been used successfully in
several recent studies jouanin2016designer ; blanchard2020multipolar ;
bertrand2020global ; vynck2022visual .

Figure 2: Sketch of the supercell approach used to compute the specular
reflectance from particle monolayers on layered substrates. Non-overlapping
particles are distributed, using a Random Sequential Addition algorithm, at
random positions in a square of side $\Lambda$ with periodic boundary
conditions. The unwanted effects of this artificial periodicity on the
scattering properties are expected to vanish as $\Lambda$ increases.

To test the validity of the ISA and EFA models, we propose here to compute the
specular reflectance spectra, $R_{\text{st}}=|r_{\text{st}}|^{2}$ of several
systems under planewave illumination using the supercell approach. The
artificial periodicity translates the scattering problem as a diffraction
problem where the radiated power is distributed in several diffraction orders,
with the 0-th order corresponding to the specular component and all the others
to the diffuse component. One expects that, upon increasing the size of the
supercell at constant particle density, the power radiated in the 0-th order
converges towards a stable value.

Let us then start our numerical study by testing the convergence of the
method. We consider two systems for this convergence test, monolayers of
either 10-nm-radius spherical gold (Au) particles or 50-nm-radius spherical
silicon (Si) particles, deposited at a surface coverage of 10% ($f=0.10$) on a
semi-infinite silica (SiO2) substrate. The particle monolayers are illuminated
by a planewave at $\lambda=500$ and $470$ nm, respectively, and at either
$\theta_{\text{i}}=0^{\circ}$ or $\theta_{\text{i}}=60^{\circ}$ in either TE
or TM polarization. The simulations are repeated on 10 independent disorder
configurations for systems up to $N=100$ particles. The results given in Fig.
3 show some stabilization around a fixed value (the black dashed lines serve
as guides to the eye), although some oscillations and fluctuations remain.
Predictions appear reasonably accurate for $N\leq 50$.

Figure 3: Convergence of the supercell method for the specular reflectance
with increasing particle number. The convergence study is done in the same
conditions as for the spectra calculations (incident angles and polarizations,
number of disorder configurations). One expects to converge to a stabilized
value as the particle number (or equivalently at fixed density, the system
size) increases. The black dashed line serves as a guide to the eye. (Left)
For a monolayer of spherical Au NPs on a SiO2 substrate at a surface coverage
$f=0.10$ and a wavelength $\lambda=500$ nm. (Right) For a monolayer of
spherical Si NPs on a SiO2 substrate at a surface coverage $f=0.10$ and at a
wavelength $\lambda=470$ nm.

We thus proceed fixing $N=50$ for all systems and compute the specular
reflectance spectra in the visible range for different systems, composed of
spherical Au or Si particles at different surface coverages on different
substrates. The results are given in Figs.4-7. The numerical predictions are
compared with those from the ISA model (Eqs. (46), (47) and (61)) and the EFA
model (Eqs. (46), (47) and (61)).

For the gold particles at low surface coverage ($f=0.01$) [Fig. 4], all model
curves superimpose with the numerical data, showing clearly that the
assumption of independent scattering is fully justified here. The peak
observed in the spectrum corresponds to the plasmon resonance of individual
particles.

The situation is different for the same system at higher surface coverage
($f=0.10$) [Fig. 5], where the agreement between the models and the numerics
is only moderately satisfactory. Remarkably, the EFA does not yield
significantly better predictions compared to the ISA, suggesting that the
mutual interaction between particles is more intricate than assumed. A reason
may be the strong near-field interaction, accompanied by the formation of “hot
spots”, which are expected in dense heterogeneous metallic nanostructures.

The spectra of high-index Si particles at the same surface coverage ($f=0.10$)
[Fig. 6] exhibit sharp spectral resonances due to the Mie resonances of the
individual particles. The strength of these resonances are yet strongly
overestimated by the ISA model, especially at grazing angles, where values up
to about 4 times the numerical value are reached. This is likely being due to
the unphysical divergence of the coefficients at large angles in the ISA. By
comparison, the EFA model performs very well: despite some inaccurate
predictions near resonance wavelengths, the spectral features are overall very
well reproduced.

Even more impressive are the results for the same particle monolayer on a
layered substrate composed of a 500-nm-thick SiO2 intermediate layer of top of
a semi-infinite Si substrate. The multiple reflections between the monolayer
and the various interfaces lead to strong spectral variations, which are well
captured by the EFA, contrary to the ISA.

All in all, our full-wave computations show that, whereas the use of the ISA
should be restricted to very dilute systems, the EFA model can be
quantitatively accurate for particle monolayers with surface coverages of
about 10% even at large angles of incidence. The specular reflection and
transmission coefficients from particle monolayers have been derived here,
step by step, from Maxwell’s equations. We hope that this derivation and the
numerical tests will be useful to the scientific community.

Figure 4: Specular reflectance spectra of a monolayer of 10-nm-radius gold
(Au) nanoparticles (NPs) deposited at a surface coverage (or filling fraction)
$f=0.01$ on top of a SiO2 substrate. The three panels correspond to three
different incident angles and polarizations: $0^{\circ}$ (polarization
independent), $60^{\circ}$ in TE-polarization and in TM-polarization. The
numerical predictions obtained by full-wave multiple-scattering computations
(square markers, results averaged over 10 disorder configurations, error bar =
standard deviation) are compared to the predictions from the ISA model
(dotted-dashed line) and the ESA model (solid line). Figure 5: Same as Fig. 4
for a higher surface coverage $f=0.10$. Figure 6: Same as Fig. 5 for 50-nm-
radius silicon (Si) particles. The particles exhibit strong Mie resonances in
the blue part of the spectrum. Figure 7: Same as Fig. 6 for a layered
substrate composed of a 500-nm-thick SiO2 layer on top of a semi-infinite Si
substrate.

##  5 Acknowledgments

We are grateful to Jean-Paul Hugonin (Laboratoire Charles Fabry, Palaiseau,
France) for providing the multiple-scattering code used in this work. This
work has received financial support from the PSA group and from the french
National Agency for Research (ANR) under the projects “NanoMiX”
(ANR-16-CE30-0008) and “NANO-APPEARANCE” (ANR-19-CE09-0014).

## References

  * (1) Petr Beckmann and Andre Spizzichino.  The scattering of electromagnetic waves from rough surfaces.  Norwood, MA, Artech House, Inc., 1987, 511 p., 1987.
  * (2) John C Stover.  Optical scattering: measurement and analysis, volume 2.  SPIE optical engineering press Bellingham, 1995.
  * (3) Dick Bedeaux and Jan Vlieger.  Optical properties of surfaces.  Imperial College Press, 2004.
  * (4) Manuel Nieto Vesperinas.  Scattering and diffraction in physical optics.  World Scientific Publishing Company, 2006.
  * (5) Alexander G Voronovich.  Wave scattering from rough surfaces, volume 17.  Springer Science & Business Media, 2013.
  * (6) Henri Benisty, Jean-Jacques Greffet, and Philippe Lalanne.  Introduction to Nanophotonics.  Oxford University Press, 2022.
  * (7) Cefe Lopez.  Materials aspects of photonic crystals.  Advanced Materials, 15(20):1679–1704, 2003.
  * (8) Danilo Mascolo, Ida Petrone, Emanuele Pascale, Marco Pisco, and Antonello Cutolo.  Production technologies for manufacturing of metamaterials and photonic crystals.  In Selected Topics In Photonic Crystals And Metamaterials, pages 505–532. World Scientific, 2011.
  * (9) Nikolay A Mirin and Naomi J Halas.  Light-bending nanoparticles.  Nano Letters, 9(3):1255–1259, 2009.
  * (10) Raquel Gomez-Medina, Braulio Garcia-Camara, Irene Suárez-Lacalle, Francisco González, Fernando Moreno, Manuel Nieto-Vesperinas, and Juan José Sáenz.  Electric and magnetic dipolar response of germanium nanospheres: interference effects, scattering anisotropy, and optical forces.  Journal of Nanophotonics, 5(1):053512, 2011.
  * (11) Tong Wu, Alexandre Baron, Philippe Lalanne, and Kevin Vynck.  Intrinsic multipolar contents of nanoresonators for tailored scattering.  Physical Review A, 101(1):011803, 2020.
  * (12) Arseniy I Kuznetsov, Andrey E Miroshnichenko, Mark L Brongersma, Yuri S Kivshar, and Boris Luk’yanchuk.  Optically resonant dielectric nanostructures.  Science, 354(6314), 2016.
  * (13) Philippe Lalanne and Pierre Chavel.  Metalenses at visible wavelengths: past, present, perspectives.  Laser & Photonics Reviews, 11(3):1600295, 2017.
  * (14) Seyedeh Mahsa Kamali, Ehsan Arbabi, Amir Arbabi, and Andrei Faraon.  A review of dielectric optical metasurfaces for wavefront control.  Nanophotonics, 7(6):1041–1068, 2018.
  * (15) Amr M Shaltout, Vladimir M Shalaev, and Mark L Brongersma.  Spatiotemporal light control with active metasurfaces.  Science, 364(6441), 2019.
  * (16) Lucien Roach, Adrian Hereu, Philippe Lalanne, Etienne Duguet, Mona Tréguer-Delapierre, Kevin Vynck, and Glenna L Drisko.  Controlling disorder in self-assembled colloidal monolayers via evaporative processes.  Nanoscale, 14:3324–3345, 2022.
  * (17) Antoine Moreau, Cristian Ciracì, Jack J Mock, Ryan T Hill, Qiang Wang, Benjamin J Wiley, Ashutosh Chilkoti, and David R Smith.  Controlled-reflectance surfaces with film-coupled colloidal nanoantennas.  Nature, 492(7427):86–89, 2012.
  * (18) Paul Chevalier, Patrick Bouchon, Julien Jaeck, Diane Lauwick, Nathalie Bardou, Alain Kattnig, Fabrice Pardo, and Riad Haïdar.  Absorbing metasurface created by diffractionless disordered arrays of nanoantennas.  Applied Physics Letters, 107(25):251108, 2015.
  * (19) Jon W Stewart, Gleb M Akselrod, David R Smith, and Maiken H Mikkelsen.  Toward multispectral imaging with colloidal metasurface pixels.  Advanced Materials, 29(6):1602971, 2017.
  * (20) Peter M Piechulla, Evgeniia Slivina, Derk Bätzner, Ivan Fernandez-Corbaton, Prerak Dhawan, Ralf B Wehrspohn, Alexander N Sprafke, and Carsten Rockstuhl.  Antireflective huygens’ metasurface with correlated disorder made from high-index disks implemented into silicon heterojunction solar cells.  ACS Photonics, 8(12):3476–3485, 2021.
  * (21) Anthony Jouanin, Jean Paul Hugonin, and Philippe Lalanne.  Designer colloidal layers of disordered plasmonic nanoparticles for light extraction.  Advanced Functional Materials, 26(34):6215–6223, 2016.
  * (22) Yidenekachew J Donie, Dominik Theobald, Somayeh Moghadamzadeh, Adrian Mertens, Ihteaz M Hossain, Ulrich W Paetzold, Uli Lemmer, and Guillaume Gomard.  Planarized and compact light scattering layers based on disordered titania nanopillars for light extraction in organic light emitting diodes.  Advanced Optical Materials, 9(14):2001610, 2021.
  * (23) Herve Bertin, Yoann Brule, Giovanni Magno, Thomas Lopez, Philippe Gogol, Laetitia Pradere, Boris Gralak, David Barat, Guillaume Demésy, and Beatrice Dagens.  Correlated disordered plasmonic nanostructures arrays for augmented reality.  ACS Photonics, 5(7):2661–2668, 2018.
  * (24) Kevin Vynck, Romain Pacanowski, Adrian Agreda, Arthur Dufay, Xavier Granier, and Philippe Lalanne.  The visual appearances of disordered optical metasurfaces.  Nature Materials, 21(9):1035–1041, 2022.
  * (25) Natalia A Loiko, Alexander A Miskevich, and Valery A Loiko.  Incoherent component of light scattered by a monolayer of spherical particles: analysis of angular distribution and absorption of light.  Journal of the Optical Society of America A, 35(1):108–118, 2018.
  * (26) Karthik Sasihithlu, Nir Dahan, Jean-Paul Hugonin, and Jean-Jacques Greffet.  A surface-scattering model satisfying energy conservation and reciprocity.  Journal of Quantitative Spectroscopy and Radiative Transfer, 171:4–14, 2016.
  * (27) T Yamaguchi, S Yoshida, and A Kinbara.  Optical effect of the substrate on the anomalous absorption of aggregated silver films.  Thin Solid Films, 21(1):173–187, 1974.
  * (28) VA Fedotov, VI Emel’yanov, KF MacDonald, and NI Zheludev.  Optical properties of closely packed nanoparticle films: spheroids and nanoshells.  Journal of Optics A: Pure and Applied Optics, 6(2):155, 2003.
  * (29) Johann Toudert, Lionel Simonot, Sophie Camelio, and David Babonneau.  Advanced optical effective medium modeling for a single layer of polydisperse ellipsoidal nanoparticles embedded in a homogeneous dielectric medium: surface plasmon resonances.  Physical Review B, 86(4):045415, 2012.
  * (30) D Bedeaux and J Vlieger.  A phenomenological theory of the dielectric properties of thin films.  Physica, 67(1):55–73, 1973.
  * (31) J Vlieger and D Bedeaux.  A statistical theory for the dielectric properties of thin island films.  Thin Solid Films, 69(1):107–130, 1980.
  * (32) Remi Lazzari and Ingve Simonsen.  Granfilm: a software for calculating thin-layer dielectric properties and fresnel coefficients.  Thin Solid Films, 419(1-2):124–136, 2002.
  * (33) Leung Tsang and Jin Au Kong.  Scattering of Electromagnetic Waves: Advanced Topics, volume 26.  John Wiley & Sons, 2004.
  * (34) Mary Carmen Peña-Gomar, Francisco Castillo, Augusto García-Valenzuela, Rubén G Barrera, and Elías Pérez.  Coherent optical reflectance from a monolayer of large particles adsorbed on a glass surface.  Applied Optics, 45(4):626–632, 2006.
  * (35) Augusto García-Valenzuela, Edahí Gutiérrez-Reyes, and Rubén G Barrera.  Multiple-scattering model for the coherent reflection and transmission of light from a disordered monolayer of particles.  Journal of the Optical Society of America A, 29(6):1161–1179, 2012.
  * (36) Chen-To Tai.  Dyadic Green functions in electromagnetic theory.  Institute of Electrical & Electronics Engineers (IEEE), 1994.
  * (37) Michael I Mishchenko, Larry D Travis, and Daniel W Mackowski.  T-matrix computations of light scattering by nonspherical particles: a review.  Journal of Quantitative Spectroscopy and Radiative Transfer, 55(5):535–575, 1996.
  * (38) Craig F Bohren and Donald R Huffman.  Absorption and scattering of light by small particles.  John Wiley & Sons, 2008.
  * (39) John David Jackson.  Classical electrodynamics, 1999.
  * (40) Lukas Novotny and Bert Hecht.  Principles of Nano-Optics.  Cambridge university press, 2012.
  * (41) Rémi Carminati and John C Schotland.  Principles of Scattering and Transport of Light.  Cambridge University Press, 2021.
  * (42) Hermann Weyl.  Ausbreitung elektromagnetischer wellen über einem ebenen leiter.  Annalen der Physik, 365(21):481–500, 1919.
  * (43) Kevin Vynck, Romain Pierrat, Rémi Carminati, Luis S Froufe-Pérez, Frank Scheffold, Riccardo Sapienza, Silvia Vignolini, and Juan José Sáenz.  Light in correlated disordered media.  arXiv:2106.13892, 2021.
  * (44) Michael Paulus, Phillipe Gay-Balmaz, and Olivier JF Martin.  Accurate and efficient computation of the green’s tensor for stratified media.  Physical Review E, 62(4):5797, 2000.
  * (45) Pochi Yeh.  Optical waves in layered media, 2nd edition.  Wiley-Interscience, 2005.
  * (46) Jean Paul Hugonin and Philippe Lalanne.  Reticolo software for grating analysis.  arXiv:2101.00901, 2021.
  * (47) Mélanie M Bay, Silvia Vignolini, and Kevin Vynck.  Pyllama: A stable and versatile python toolkit for the electromagnetic modelling of multilayered anisotropic media.  Computer Physics Communications, 273:108256, 2022.
  * (48) Tom G Mackay and Akhlesh Lakhtakia.  The transfer-matrix method in electromagnetics and optics.  Morgan and Claypool, San Rafael, CA, USA, 2020.
  * (49) D Yuk Kei Ko and JR Sambles.  Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals.  Journal of the Optical Society of America A, 5(11):1863–1866, 1988\.
  * (50) Maxime Bertrand, Alexis Devilez, Jean-Paul Hugonin, Philippe Lalanne, and Kevin Vynck.  Global polarizability matrix method for efficient modeling of light scattering by dense ensembles of non-spherical particles in stratified media.  Journal of the Optical Society of America A, 37(1):70–83, 2020.
  * (51) Amos Egel, Krzysztof M Czajkowski, Dominik Theobald, Konstantin Ladutenko, Alexey S Kuznetsov, and Lorenzo Pattelli.  Smuthi: A python package for the simulation of light scattering by multiple particles near or between planar interfaces.  Journal of Quantitative Spectroscopy and Radiative Transfer, 273:107846, 2021.
  * (52) Mathieu Langlais, Jean-Paul Hugonin, Mondher Besbes, and Philippe Ben-Abdallah.  Cooperative electromagnetic interactions between nanoparticles for solar energy harvesting.  Optics Express, 22(103):A577–A588, 2014.
  * (53) Cédric Blanchard, Jean-Paul Hugonin, Alima Nzie, and Domingos De Sousa Meneses.  Multipolar scattering of subwavelength interacting particles: Extraction of effective properties between transverse and longitudinal optical modes.  Physical Review B, 102(6):064209, 2020.