# Effective T-matrix of a cylinder filled with a random 2D particulate

K. Napal, P. Piva, A. L. Gower

Effective T-matrix of a cylinder filled
with a random 2D particulate

K. K. Napal1, P. S. Piva1 and A. L. Gower1

1 Department of Mechanical Engineering, University of Sheffield, UK
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

(August 28, 2024)

Abstract. When a wave, such as sound or light, scatters within a densely
packed particulate, it can be rescattered many times between the particles,
which is called multiple scattering. Multiple scattering can be unavoidable
when: trying to use sound waves to measure a dense particulate, such as a
composite with reinforcing fibers. Here we solve from first principles
multiple scattering of scalar waves, including acoustic, for any frequency
from a set of 2D particles confined in a circular area. This case has not been
solved yet, and its solution is important to perform numerical validation, as
particles within a cylinder require only a finite number of particles to
perform direct numerical simulations. The method we use involves ensemble
averaging over particle configurations, which leads us to deduce an effective
T-matrix for the whole cylinder, which can be used to easily describe the
scattering from any incident wave. In the specific case when the particles are
monopole scatters, the expression of this effective T-matrix simplifies and
reduces to the T-matrix of a homogeneous cylinder with an effective wavenumber
$k_{\star}$. To validate our theoretical predictions we develop an efficient
Monte Carlo method and conclude that our theoretical predictions are highly
accurate for a broad range of frequencies.

Subject Areas. wave motion, acoustics, statistical physics

Keywords. multiple scattering, random media, effective media, ensemble average

##  1 Introduction

Ensemble averaging. Multiple scattering is unavoidable when using waves to
characterise a particulate composite, or designing metamaterials to control
wave propagation. Further, the number of particles in most applications makes
direct numerical simulations impossible for current computing power, though
there are some notable attempts [13, 15]. Even if such simulations were
possible, the positions of the particles are impossible to know for one
particular sample. A mathematical method to solve this problem, called
ensemble average, is to take an average over all possible particle positions.
This type of averaging can occur naturally for fluids and gases, as the
particles are moving. That is, taking an average of the scattered signal over
time can be equivalent to an ensemble average111In ergodic systems, if enough
time has passed, all physically possible states of the system will have
occurred, and so that taking an average over time is equivalent to averaging
over all possible configurations.. Assuming that ensemble averaging is
equivalent to averaging over time and space, which is often the case, is
called the ergodicity assumption [22].

What is known. The scenarios that are best understood are: 1) waves in an
infinite media with no boundaries [27, 36], and 2) plane waves incidence on a
half-space or plate filled with particles [33, 34]. Both scenarios have been
considered to obtain effective wavenumbers [6, 19, 37]. Though we mention the
methods that use Lipmann-Schwinger for acoustics, they involve an extra
integral which is often omitted [20] and complicates the calculations. Our
approach in this work is valid for any type of scatterer and scalar waves.
Both scenarios of using plane waves and an infinite media have several
applications (typically when considering layered media such as planetary
atmospheres, layers in the ocean, or soils) but one significant drawback is
that it has been very challenging to numerically validate the assumptions used
for these methods. Both use statistical assumptions, such as the Quasi-
Crystalline Assumption (QCA) [17] that is not based on an asymptotic
approximation. Validation is needed to establish the range of validity of
these assumptions. However, direct numerical simulations of scattering from a
configuration of particles for both planes and infinite media require a huge
number of particles [14], [9] or the introduction of periodic boundaries which
can introduce artifacts [3].

The cylindrical setting. The methods developed to describe the average plane
wave propagating in a disordered particulate plate or half-space can now be
extended to other geometries [10]. The ideal scenario to compare theoretical
predictions with direct numerical simulations is to have cylindrical particles
inside a cylinder, as this reduces the problem to two dimensions and we need
only a finite number of particles for the direct numerical simulations. See
Figure 1 for an illustration. This is the simplest case to perform numerical
validation of a very general theory [10]. Further, we show in this work that
the effective dispersion equation for the cylindrical geometry is the same as
the plane-wave case. So numerically validating the cylindrical geometry will
also serve as numerical validation of the dispersion equation for plane waves
and all geometries. We also note that it appears that the cylindrical setting
has never been solved from a first principles approach. The formulas we
provide are also valid for any inter-particle pair correlation.

Industrial applications. Beyond numerical validation, there are industrial
applications that need a method to calculate waves scattered from a cylinder
with particles. Examples of cylinders filled with cylindrical particles
include concrete beams reinforced with iron, cables filled with wires, or
fibre-reinforced composite [1]. Applications include designing cylinders with
exotic effective properties or developing methods to measure the cylindrical
particles [24, 25]. In terms of measurement, it is likely that more
information can be extracted from waves scattered from a cylinder filled with
a particulate than just plane wave reflection from a plate filled with the
same particulate.

Effective properties. The most common approach to model the average scattering
from, say, a spherical or cylindrical region with particles is to assume the
region is homogeneous with some effective properties [31, 26, 5, 18, 30, 32,
35], and then use the standard boundary conditions such as continuity of
displacement. This approach is valid for low-frequency [31] but for finite
frequencies is incorrect in three dimensions [10], and we demonstrate the same
here for two dimensions in this work. To obtain an accurate model for a finite
frequency the boundary condition needs to be deduced from first principles,
together with an eigensystem for the effective wavenumbers [10]. Although the
process is more involved, the final expression for the average scattered wave
from a cylindrical region is simple: the average scattered wave can be
calculated from an effective T-matrix for any incident frequency, source, and
particle properties. We stress that without deducing the results from first
principles, as we do here, it would not be possible to just guess the form of
this effective T-matrix.

Monopole scatterers. One somewhat surprising result we deduce is that if the
particles scatter only monopole waves, that is waves that have radial
symmetry, then the material as a whole behaves as homogeneous, where the mass
density is the same as the background, and the bulk density is given by a
simple formula. We deduce this for particles in a cylinder and hypothesize
that it is true for any material filled with monopole scatterers, even when
including all orders of multiple scattering. Beyond just curiosity, there are
many particles that behave approximately like a monopole scatterer, and
therefore the simple formulas we deduced are appropriate. For example, in
acoustics void-like particles are approximately monopoles for a broad
frequency range, see the Dirichlet case in Figure 2. In elasticity, particles
become approximately monopole when the bulk modulus is much greater than the
shear modulus [4]. Other cases include resonators such as a split ring
resonator [28].

Overview of the method. After ensemble averaging over particle configurations
within a cylindrical region, the system inherits cylindrical symmetry. For
example, if the source has radial symmetry then the average scattered field
will also have radial symmetry. This is also true for sources with more
general rotational symmetry resulting in scattered fields of the same
rotational symmetry (cf. Figure 1). In this paper, we take advantage of this
mode-to-mode symmetry to analyze the general behaviour of the random
particulate material independently from the incident field.

After denoting with $\mathrm{V}_{n}$ and $\mathrm{U}_{n}$ the regular and
outgoing cylindrical waves of order $n$ (cf. 2.2), the cylindrical symmetry
translates as follows: when the exciting source is $\mathrm{V}_{n}$, then the
average scattered field is $\mathrm{T}_{n}\mathrm{U}_{n}$ where the complex
number $\mathrm{T}_{n}$ only depends on the properties of the random
particulate cylinder (radius and properties of the particles). Since the
scattering problem is linear, the knowledge of the $\mathrm{T}_{n}$ allows us
to describe the scattering from any incident field, after decomposing the
latter into the modes $\mathrm{V}_{n}$. Having simple expressions for
$\mathrm{T}_{n}$ is crucial to help guide methods to characterize or design
particulate materials. We do so by using the effective wave method approach
[10]. Finally, we validate our results with an adapted Monte Carlo method in
which the rate of convergence is accelerated thanks to the cylindrical
symmetry.

The mode to mode scattering
Figure 1:  Consider a circular region filled with sound-hard particles. Both
panels compare the average scattered wave $\langle u_{\mathrm{sc}}\rangle$,
for a source with only one mode: $u_{\mathrm{inc}}=\mathrm{V}_{N}(\bm{r})$,
when using our effective wave method and a brute force Monte Carlo approach.
The Monte Carlo approach simulates the scattered field from either one, $M=1$,
or two hundred, $M=200$, configurations of particles and then takes the
average over these configurations. In general, the average scattered wave is
given by (2.13). The left panel shows that when using an incident wave
$u_{\mathrm{inc}}=\mathrm{V}_{0}(\bm{r})$, with radial symmetry, then the only
scattered mode also has radial symmetry after average. The right panel shows
how a source with a $60^{\circ}$ rotation symmetry also leads to a scattered
wave with the same symmetry.

Overview of this paper

In Section 2, we first introduce the statistics of the random particulate
material and the required notations for the ensemble averaging. We then define
the T-matrix of the effective cylinder whose exact formula depends on the
solutions of the averaged Foldy-Lax equations.

The latter is solved in Section 3 by using the effective waves method, which
consists in finding solutions that are isotropic waves with a complex
wavenumber $k_{\star}$. The method leads to an eigenvalue problem called the
dispersion equation, whose eigenvalue provides $k_{\star}$.

In Section 4, we use the expression of the solutions of the averaged Foldy-Lax
equations to deduce a formula of the effective T-matrix. The latter is very
simple when the particles are monopole scatterers:

$\mathrm{T}_{n}=-\frac{{\cal C}_{n}}{{\cal
D}_{n}}\qquad\text{with}\qquad\left\\{\begin{array}[]{rcl}{\cal
C}_{n}&=&\displaystyle
k\mathrm{J}^{\prime}_{n}(k\tilde{R})\mathrm{J}_{n}(k_{\star}\tilde{R})-k_{\star}\mathrm{J}_{n}(k\tilde{R})\mathrm{J}^{\prime}_{n}(k_{\star}\tilde{R})\\\
{\cal D}_{n}&=&\displaystyle
k\mathrm{H}^{\prime}_{n}(k\tilde{R})\mathrm{J}_{n}(k_{\star}\tilde{R})-k_{\star}\mathrm{H}_{n}(k\tilde{R})\mathrm{J}^{\prime}_{n}(k_{\star}\tilde{R})\end{array}\right.$
(1.1)

where $\tilde{R}:=R-a$. This result is remarkable because the above expression
corresponds to the T-matrix of a homogeneous acoustic cylinder of radius
$\tilde{R}$, sound speed $c_{\star}=\omega/k_{\star}$ and density $\rho$ which
is equal to the background medium. Particles that are approximately monopole
scatterers appear in mainly two scenarios: either small sound soft particles
or resonators [29]. For all these cases, the T-matrix (1.1) shows us that the
effective wavenumber $k_{\star}$ suffices to describe the random material.

When the particles are not monopole scatterers, (1.1) is not exact. The exact
formula is given by

$\mathrm{T}_{n}=-\frac{\displaystyle\sum_{n^{\prime}}{\cal
C}_{n-n^{\prime}}F_{n^{\prime}}}{\displaystyle\sum_{n^{\prime}}{\cal
D}_{n-n^{\prime}}F_{n^{\prime}}}$ (1.2)

where the ${\cal C}_{n}$ and ${\cal D}_{n}$ are the same as before, and the
weights $F_{n}$ are the eigenfunctions of the dispersion equation, associated
with the effective wavenumber $k_{\star}$.

For monopole scatterers, we have that ($n^{\prime}=0$), and the above reduces
to (1.1). We note that often $F_{0}$ is the largest term, which explains why
(1.1) can give accurate results for non-monopole scatterers, see for example
the numerical results for sound soft (Dirichlet) particles shown in Figure 2
in the low-frequency regime. In this same figure, we also see how for sound
hard particles (Neumann), which are dipole scatters, the Monte Carlo results
and the formula (1.2) closely match, whereas the T-matrix for monopole
scatters doesn’t match.

Validation of 1.1 and 1.2
Figure 2: Compares various methods to calculate the component $\mathrm{T}_{0}$
of the T-matrix of a cylinder filled with particles. The solid red line is our
effective waves method (EWM) (4.21), the black points are from a brute force
Monte Carlo method (MC), and the dashed blue line is our method when only
monopole scattering is accounted for (EWM-MA) (1.1). The left (right) graph
shows the results for sound soft (hard) particles. In both cases, the general
expression of the effective T-matrix matches the Monte Carlo results. The EWM-
MA method only matches well with the Monte Carlo for sound soft scatterers for
low frequencies.

Set up of cylindrical particles inside a cylinder
$\mathbf{r}_{2}$$\bullet$$\mathbf{r}_{1}$$\bullet$$\mathbf{r}_{3}$$\bullet$$\mathbf{r}_{4}$$\bullet$$\tilde{R}=R-a_{1}$$R$$\mathcal{R}_{1}$$O$$\bullet$$a_{1}$$x$$z$$y$
Figure 3: Illustration of a possible configuration of four particles. The
particles are cylinders whose axes are aligned with the $z$-axis. The position
$\bm{r}_{i}\in\mathbb{R}$ of the $i^{\text{th}}$ particle is determined by the
intersection of its axis with the plane $xOy$. The particles are of radius $a$
and physically contained in a cylinder of radius $R$. In this specific
example, all the particles are the same and the centres $\bm{r}_{i}$ are
therefore confined in the same disc ${\cal R}_{1}$ shown in green.

##  2 Random particulate material

###  2.1 Deterministic scattering from J particles

Here we summarise some results for multiple scattering of acoustic waves, in
the time-harmonic regime, by a collection of $J$ cylinders, referred as to
particles, which axis are aligned with the $z$-axis. The centre of the
particle $i$ is identified by $\bm{r}_{i}\in\mathbb{R}^{2}$ as shown in Figure
3 and is assumed to be confined in an area denoted by ${\cal R}_{i}$.

The propagation of waves in free space is governed by the 2D Helmholtz
equation

$\Delta u+k^{2}u=0$ (2.1)

where $\Delta:=\partial_{x}^{2}+\partial_{y}^{2}$ is the 2D Laplace operator
and $k\in\mathbb{R}$ is the wavenumber of the homogeneous background. Consider
the following two independent bases of the Helmholtz equation $\mathrm{V}_{n}$
and $\mathrm{U}_{n}$ defined by

$\begin{array}[]{|rcll}\mathrm{U}_{n}(k\bm{r})&:=&\mathrm{H}_{n}(kr)\mathrm{e}^{\mathrm{i}n\theta}&\forall\bm{r}\in\mathbb{R}^{2}\setminus\\{\bm{0}\\}\\\
\mathrm{V}_{n}(k\bm{r})&:=&\mathrm{J}_{n}(kr)\mathrm{e}^{\mathrm{i}n\theta}&\forall\bm{r}\in\mathbb{R}^{2}\end{array}$
(2.2)

where $k$ is the background wavenumber and $(r,\theta)$ are the polar
coordinates of $\bm{r}$, ie $\bm{r}=(r\cos\theta,r\sin\theta)$. The specific
solutions $\mathrm{V}_{n}(k\bm{r})$ have the particularity of being smooth
while $\mathrm{U}_{n}(k\bm{r})$ have a singularity at the origin and are
outgoing solutions.

Both the incident field $u_{\mathrm{inc}}$ and the scattered field
$u_{\mathrm{sc}}$ are solutions of (2.1). While we assume that
$u_{\mathrm{inc}}$ is smooth and regular in the region that covers the
particles222This is true for any source which originated from outside this
region. the scattered field has to be a sum of outgoing fields from each
particles centered at $\bm{r}_{i}$, as a result we can write

$u_{\mathrm{inc}}(\bm{r})=\sum_{n=-\infty}^{+\infty}g_{n}\mathrm{V}_{n}(k\bm{r})$
(2.3a)
$u_{\mathrm{sc}}(\bm{r})=\sum_{i=1}^{J}\sum_{n=-\infty}^{+\infty}f_{n}^{i}\mathrm{U}_{n}(k\bm{r}-k\bm{r}_{i}).$
(2.3b)

The coefficients $g_{n}\in\mathbb{C}$ are given parameters of the system and
$f_{n}^{i}\in\mathbb{C}$ are unknowns that can be determined by following the
Foldy-Lax self-consistent method [21]. The latter leads to the following
system of equations

$f_{n}^{i}=\mathrm{T}_{n}^{i}\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{i})g_{n^{\prime}}+\mathrm{T}_{n}^{i}\sum_{j\neq
i}\sum_{n^{\prime}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{i}-\bm{r}_{j})f_{n^{\prime}}^{j}\quad\forall
n\in\mathbb{Z},\,\forall i=1\dots J,$ (2.4)

where $\mathrm{T}_{n}^{i}$ is the T-matrix of particle $i$, which can
represent a wide range of particles [7].

To give an example, the expression of a T-matrix of a homogeneous particle
with wavenumber $k_{i}$, density $\rho_{i}$ and radius $a_{i}$ is given by

$\mathrm{T}^{i}_{n}=-\frac{\rho_{i}k\mathrm{J}_{n}^{\prime}(ka_{i})\mathrm{J}_{n}(k_{i}a_{i})-\rho
k_{i}\mathrm{J}_{n}(ka_{i})\mathrm{J}_{n}^{\prime}(k_{i}a_{i})}{\rho_{i}k\mathrm{H}^{\prime}_{n}(ka_{i})\mathrm{J}_{n}(k_{i}a_{i})-\rho
k_{i}\mathrm{H}_{n}(ka_{i})\mathrm{J}_{n}^{\prime}(k_{i}a_{i})}\cdot$ (2.5)

The system (2.4) totally determines the scattering coefficients $f_{n}^{i}$,
which allows to solve the scattering problem from the $J$ given particles.
However, this result is not very useful in practice for two main reasons:
first, the position of the particles is often unknown, and second, there can
be a very large number of particles in most industrial applications [2]. The
ensemble average over all particle positions, that we summarize below, solves
both these problems from the computational standpoint.

###  2.2 Particulate distribution

We describe each particle, say particle $i$, with two random variables:
$\bm{r}_{i}$ is the probability distribution of the centre of the particle in
space, and $\lambda_{i}\in{\cal S}$ describes all other properties of the
particle (radius $a_{i}$, density $\rho_{i}$, etc.).

The centre of a particle $\bm{r}_{i}$ is restricted in a region ${\cal R}_{i}$
which depends on the particle’s properties $\lambda_{i}$. For example,
$\bm{r}_{i}$ has to be one radius away from the boundary $\partial{\cal R}$
which completely contains all particles. In this paper we assume that
$\bm{r}_{i}$ is distributed uniformly over the set ${\cal R}_{i}(\lambda_{i})$
defined by

${\cal R}_{i}:=\\{\bm{r}\in{\cal R}\,:\,\mathrm{dist}(\bm{r},\partial{\cal
R})>a_{i}\\}.$ (2.6)

As a consequence,

$\mathrm{p}(\bm{r}_{i},\lambda_{i})=\mathrm{p}(\bm{r}_{i}|\lambda_{i})\mathrm{p}(\lambda_{i})=\frac{\mathrm{p}(\lambda_{i})}{|{\cal
R}_{i}|},$ (2.7)

where $|{\cal R}_{i}|$ is the volume of ${\cal R}_{i}$.

#### Joint particle distribution

Each particle is described by two random variables $(\lambda_{i},\bm{r}_{i})$,
$i=1,\dots,J$. In general, the particle positions are correlated, for example,
no particles can overlap. The most commonly used term to describe inter-
particle correlation is the pair correlation $\mathrm{g}$, which satisfies
(cf. [16, eq. 8.1.2]):

$\mathrm{p}(\bm{r}_{1};\bm{r}_{2}|\lambda_{1};\lambda_{2})=\frac{\mathrm{g}(\bm{r}_{1},\lambda_{1};\bm{r}_{2},\lambda_{2})}{|{\cal
R}_{1}||{\cal R}_{2}|}\frac{J}{J-1},$ (2.8)

where on the left is the joint law of two particles’ position when their
properties are known.

The pair correlation $\mathrm{g}$ describes how correlated any two particles
are, when the positions and properties of all other particles are unknown. For
example, if $\mathrm{g}=1$ for all values of its arguments then both
$\bm{r}_{i}$ and $\bm{r}_{j}$ in the above are independent and uniformly
distributed over ${\cal R}_{i}$ and ${\cal R}_{j}$, respectively (in the limit
$J\to\infty$).

Finally, we introduce the density

$\mathfrak{n}(\lambda_{i}):=\frac{J}{|{\cal
R}_{i}|}\mathrm{p}(\lambda_{i})\quad\text{(number of $\lambda_{i}$ types
particles per unit volume).}$ (2.9)

Then we derive the following useful relation

$\mathrm{p}(\bm{r}_{j},\lambda_{j}|\bm{r}_{i},\lambda_{i})=\frac{\mathrm{p}(\bm{r}_{j},\lambda_{j};\bm{r}_{i},\lambda_{i})}{\mathrm{p}(\bm{r}_{i},\lambda_{i})}=|{\cal
R}_{i}|\mathrm{p}(\lambda_{j})\mathrm{p}(\bm{r}_{i};\bm{r}_{j}|\lambda_{i};\lambda_{j})=\frac{\mathfrak{n}(\lambda_{j})}{J-1}\mathrm{g}(\bm{r}_{i},\bm{r}_{j}),$
(2.10)

where we used (2.7), (2.8) and (2.9).

#### The pair correlation

In this paper, we consider that the particles have a distribution which is
isotropic and homogeneous in space. As a consequence, the pair correlation is
of the form
$\mathrm{g}(\bm{r}_{i},\lambda_{i};\bm{r}_{j},\lambda_{j})=\mathrm{g}(|\bm{r}_{j}-\bm{r}_{i}|,\lambda_{i},\lambda_{j})$.
We assume the pair correlation is of the following form:

$\mathrm{g}(r,\lambda_{1},\lambda_{2})=\begin{cases}0,\quad&r<a_{12},\\\
1+\delta\mathrm{g}(r,\lambda_{1},\lambda_{2}),\quad&a_{12}<r<b_{12},\\\
1,\quad&r>b_{12},\end{cases}$ (2.11)

where $g(r,\lambda_{1},\lambda_{2})=0$ when particles overlap, with
$a_{12}\geq a_{1}+a_{2}$ and in general $a_{12}$ can depend on $\lambda_{1}$
and $\lambda_{2}$. This form of the pair correlation moreover assumes that at
a certain distance $b_{12}$ from each other, the particles become
uncorrelated, that is $g(r,\lambda_{1},\lambda_{2})=1$ for $r>b_{12}$. This
assumption will lead to analytic simplifications, as well as being a good
approximation for most disordered materials. A typical plot of the pair
correlation is illustrated in Figure 4.

Radial pair-correlation function
Figure 4: Typical plot of the pair correlation function.  The pair-correlation
$\mathrm{g}(r)$ is zero for $r<a_{12}$ because particles cannot overlap. The
local minima (maxima) indicate that there is a distance $r$ from one fixed
particle, say at $\bm{r}_{1}$, where it is less (more) likely to find other
particles. Finally, assuming $\mathrm{g}(r)=1$ for $r\geq b_{12}$ means that a
particle at $\bm{r}_{1}$ becomes uncorrelated to particles which are further
than $b_{12}$.

###  2.3 Definition of the effective T-matrix

The scattered field (2.3b) can be simplified after using Graf’s addition
theorem (A.1,ii) with $\bm{x}=k\bm{r}$ and $\bm{d}=-k\bm{r}_{i}$

$u_{\mathrm{sc}}(\bm{r})=\sum_{n=-\infty}^{+\infty}\mathfrak{F}_{n}\mathrm{U}_{n}(k\bm{r})$
(2.12a)
$\mathfrak{F}_{n}:=\sum_{i=1}^{J}\sum_{n^{\prime}=-\infty}^{+\infty}\mathrm{V}_{n^{\prime}-n}(-k\bm{r}_{i})f_{n^{\prime}}^{i}.$
(2.12b)

Note that (2.12) is only valid for $r>R$, nevertheless this is enough for the
definition of the effective T-matrix below. Taking the ensemble average of the
equation above, as defined in Apendix B, leads to

$\langle
u_{\mathrm{sc}}(\bm{r})\rangle=\sum_{n=-\infty}^{+\infty}\langle\mathfrak{F}_{n}\rangle\mathrm{U}_{n}(k\bm{r}).$
(2.13)

Since the scattering problem is linear with respect to the incident field,
$\langle\mathfrak{F}_{n}\rangle$ depends linearly on the modes $g_{n}$ of the
incident field (2.3a):

$\displaystyle\langle\mathfrak{F}_{n}\rangle=\sum_{N=-\infty}^{+\infty}{\cal
T}_{n,N}g_{N}\quad\text{(effective T-matrix definition)}.$ (2.14)

This relation defines the T-matrix ${\cal T}$ of the averaged material which
connects the modes of the incident field with the ones of the averaged
scattered field. It allows us to describe the scattering from any incident
field, provided the coefficients $g_{n}$ in (2.3a).

In the specific case when the area where the particles are confined is
circular, we have $\mathcal{T}_{n,N}=0$ if $n\neq N$ (cf. Appendix E), so that
(2.13) becomes

$\langle
u_{\mathrm{sc}}(\bm{r})\rangle=\sum_{n=-\infty}^{+\infty}\mathrm{T}_{n}g_{n}\mathrm{U}_{n}(k\bm{r}).$
(2.15)

where $\mathrm{T}_{n}:=\mathcal{T}_{n,n}$.

Computing ${\cal T}_{n,N}$ requires to compute
$\langle\mathfrak{F}_{n}\rangle$. From (2.12a), (B.1), and (B.3) we obtain

$\langle\mathfrak{F}_{n}\rangle=\int_{\cal
S}\mathfrak{n}(\lambda_{1})\int_{{\cal
R}_{1}}\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(-k\bm{r}_{1})\langle
f_{n^{\prime}}\rangle(\bm{r}_{1},\lambda_{1})\,\mathrm{d}\bm{r}_{1}\mathrm{d}\lambda_{1}.$
(2.16)

Here we used (B.7) to substitute $\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\langle
f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})$.

The function $\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ needs to be
determined before we can compute $\langle\mathfrak{F}_{n}\rangle$. In Appendix
B.1 we show how to obtain the governing equation:

$\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{1})g_{n^{\prime}}\\\
+\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime}}\rangle(\bm{r}_{1},\lambda_{1})\mathrm{g}(|\bm{r}_{1}-\bm{r}_{2}|,\lambda_{1}\lambda_{2})\,\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2},$
(2.17)

where we also used the simpler pair-correlation (2.11).

In the next section, we use the _effective wave method_ to solve (2.17). This
method introduced in [10] proved to be successful in 3D and provides a closed
formula for $\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})$. This in turn
allows us to compute $\langle\mathfrak{F}_{n}\rangle$ (2.16) and reach an
explicit formula for $\mathcal{T}_{n,N}$ by specifying $g_{n}=\delta_{n,N}$ in
(2.14).

##  3 The effective wave method

To solve the general integral governing equation (2.17), we use the effective
waves method [10] as summarised below.

HTML]E9F0E9  Overview of the effective waves method The starting point is to
assume that there exists $k_{\star}\in\mathbb{C}$ such that
$(\Delta+k_{\star}^{2})\,\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})=0.$
(3.1) This assumption will greatly simplify the governing integral equation
(2.17), from which we will be able to determine both $k_{\star}$ and $\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})$. To summarise, the method has three
steps: 1. Separate microstructure and boundary. The assumption (3.1) is used
in the governing equation (2.17) to derive two separate equations called the
ensemble wave equation and the ensemble boundary conditions. The first one
only depends on the microstructure of the random material while, in contrast,
the second one takes into account the incident field and the shape of the
random material, acting much like a boundary condition. 2. The effective
eigensystem. We decompose $\langle f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})$ in
the basis of functions $\mathrm{V}_{n}(k_{\star}\bm{r}_{1})$ and substitute
the decomposition into the ensemble wave equation. The coefficients
$\mathbf{F}(\lambda_{1})$ of the decomposition are then shown to be the
eigenfunctions of an eigensystem which eigenvalue is $k_{\star}$. 3. The
ensemble boundary condition. To determine the amplitudes of the eigenfunctions
$\mathbf{F}(\lambda_{1})$ we then use the ensemble boundary condition, which
takes into account the incident field and the shape of the random material.

###  3.1 Separate microstructure and boundary

We follow the steps described above to determine the solutions $\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ of (2.17). The first step uses (3.1),
and some algebraic manipulations shown in Appendix C.1 to rewrite the
governing equation (2.17) into two separate equations:

(3.2) (3.3)

The terms ${\cal J}_{n^{\prime}n}(\bm{r}_{1})$, ${\cal
I}_{n^{\prime}n}(\bm{r}_{1})$, and ${\cal K}_{n^{\prime}n}(\bm{r}_{1})$,
respectively defined in (C.6) and (C.3), involve the function $\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})$.

One of the key advantages of splitting the integral equation (2.17) into the
two separate equations above is that the ensemble wave equation (3.2) does not
depend on the shape of the region ${\cal R}$ nor the incident wave. As we will
see below, the ensemble wave equation (3.2) can be used to determine the
effective wavenumber $k_{\star}$, which implies that $k_{\star}$ only depends
on the microstructure: the density of the particles
$\mathfrak{n}(\lambda_{1})$, their properties provided by the T-matrix
$\mathrm{T}_{n}(\lambda_{1})$, and the pair correlation $\mathrm{g}$ which
explicitly appears in the quantity ${\cal K}_{n,n^{\prime}}(\bm{r}_{1})$
(C.3). We further discuss how to interpret $k_{\star}$ in Section 4.5.

On the other hand, the ensemble boundary condition (3.3) acts like a boundary
condition, and shows how the incident wave and material boundary affect the
overall solution.

###  3.2 The effective eigensystem

Since $\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ satisfies (3.1), it can
be decomposed into the modes

$\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\sum_{n_{1}}F_{nn_{1}}(\lambda_{1})\mathrm{V}_{n_{1}}(k_{\star}\bm{r}_{1}),$
(3.4)

where $\mathrm{V}_{n_{1}}$ is defined in (2.2).

The unknowns $k_{\star}$ and $F_{nn_{1}}(\lambda_{1})$ can be determined by
substituting (3.4) into (3.2). The details are shown in Appendix C.2 with the
resulting equation being

$F_{nn_{1}}(\lambda_{1})+\sum_{n^{\prime}n_{2}}\delta_{n_{2}-n_{1}+n^{\prime}-n}\mathrm{T}_{n}(\lambda_{1})\int_{\cal
S}{\cal
N}_{n^{\prime}-n}^{12}(k,k_{\star})F_{n^{\prime}n_{2}}(\lambda_{2})\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0,$
(3.5)

where

$\displaystyle{\cal
N}_{l}^{12}(k,k_{\star}):=2\pi\frac{\mathrm{N}_{l}(ka_{12},k_{\star}a_{12})}{k_{\star}^{2}-k^{2}}-2\pi\int_{a_{12}}^{b_{12}}\mathrm{J}_{l}(k_{\star}r)\mathrm{H}_{l}(kr)\delta\mathrm{g}(r)r\mathrm{d}r,$
(3.6)
$\displaystyle\mathrm{N}_{l}(x,y):=x\mathrm{H}^{\prime}_{l}(x)\mathrm{J}_{l}(y)-y\mathrm{H}_{l}(x)\mathrm{J}^{\prime}_{l}(y).$
(3.7)

The above is a nonlinear eigenvalue problem, which is why we refer to
$k_{\star}$ as an eigenvalue. After calculating $k_{\star}$ we can calculate
the eigenfunctions $F_{n^{\prime}n_{2}}(\lambda_{1})$ by solving the linear
system (3.5), though in practice it is far better to calculate both
$k_{\star}$ and $F_{n^{\prime}n_{2}}(\lambda_{1})$ using the modal
decomposition as we do in Section 4.

Multiple effective wavenumbers
Figure 5: Shows 11 of the eigenvalues $k_{p}$ of (4.13) for sound soft
particles of radius $a=1$ confined in a cylinder of radius $R=10$ when
$\omega=1$. $k_{\star}$ corresponds to the one with the smallest imaginary
part. The $x$-axis is the real part and $y$-axis is the imaginary part. The
other wavenumbers have a much larger imaginary part leading to evanescent
waves inside the random material.

###  3.3 The ensemble boundary condition

The eigensystem (3.5) is not enough to fully determine the $F_{nn_{1}}$, for
instance if $F_{nn_{1}}$ is a solution to the eigensystem then so is $\alpha
F_{nn_{1}}$ for any scalar $\alpha$.

To fully determine $F_{nn_{1}}$ we need to substitute (3.4) into the ensemble
boundary condition (3.3). The details are shown in Appendix D, where we deduce
the ensemble boundary conditions (D.5) without specifying the region ${\cal
R}_{2}$, followed by choosing ${\cal R}_{2}$ to be a cylinder with radius
$R-a_{2}$ which results in the boundary condition:

$\displaystyle
g_{N}+\frac{2\pi}{k_{\star}^{2}-k^{2}}\sum_{n^{\prime}}\int_{\cal
S}F_{n^{\prime}(N-n^{\prime})}(\lambda_{2})N_{N-n^{\prime}}(kR_{2},k_{\star}R_{2})\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0,$
(3.8)

where $N_{l}$ is defined by (3.6).

##  4 An effective cylinder

###  4.1 Modal decomposition of the problem

In this section, we assume that the particles are confined in a circular disc
${\cal R}$ with radius $R$. Since (2.17) is linear with respect to $g_{n}$, it
can be decomposed in simpler and independent equations: let $\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})$ be the solution when substituting
$g_{n}=\delta_{n,N}$ in (2.17), then

$\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\displaystyle\mathrm{V}_{N-n}(k\bm{r}_{1})\\\
+\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime},N}\rangle(\bm{r}_{1},\lambda_{1})\mathrm{g}(|\bm{r}_{1}-\bm{r}_{2}|,\lambda_{1},\lambda_{2})\,\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2},$
(4.1)

where we can recover the solution for any incident wave by using

$\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\sum_{N}g_{N}\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1}).$ (4.2)

From the 3D version of these effective equations [10] we know that (3.5) can
be reduced by using symmetry. We show how to do this for the modal
decomposition below, with the result being:

$\displaystyle\langle
f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})=\alpha_{N}F_{n,N}(\lambda_{1})\mathrm{V}_{N-n}(k_{\star}\bm{r}_{1})$
(4.3)

where $\alpha_{N}\in\mathbb{C}$ is some amplitude which is introduced for
later convenience. Below is the proof of (4.3).

###### Proof.

Using the rotational symmetry of the modal source, we simplify the form of the
modal solutions $\langle f_{n,N}\rangle(\bm{r}_{1},\lambda_{1})$. To this end,
we denote by $\mathbf{R}_{\phi}$ the rotation matrix of angle $\phi$ and
replace $\bm{r}_{1}$ with $\mathbf{R}_{\phi}\bm{r}_{1}$ in (4.1):

$\langle
f_{n,N}\rangle(\mathbf{R}_{\phi}\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\displaystyle\mathrm{V}_{N-n}(k\mathbf{R}_{\phi}\bm{r}_{1})\\\
+\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\mathbf{R}_{\phi}\bm{r}_{1}-k\mathbf{R}_{\phi}\bm{r}_{2})\langle
f_{n^{\prime},N}\rangle(\mathbf{R}_{\phi}\bm{r}_{2},\lambda_{2})\mathrm{g}(|\bm{r}_{1}-\bm{r}_{2}|)\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2},$
(4.4)

where we changed the integration variable from $\bm{r}_{2}$ to
$\mathbf{R}_{\phi}\bm{r}_{2}$, which is possible for any rotation as ${\cal
R}_{2}$ is a disc. Then from (2.2) we deduce the property

$\mathrm{U}_{n}(\mathbf{R}_{\phi}\bm{r}_{1})=\mathrm{U}_{n}(\bm{r}_{1})\mathrm{e}^{\mathrm{i}n\phi}\quad\text{and}\quad\mathrm{V}_{n}(\mathbf{R}_{\phi}\bm{r}_{1})=\mathrm{V}_{n}(\bm{r}_{1})\mathrm{e}^{\mathrm{i}n\phi},$

which used in (4.4), and then multiplying both sides of the equation with
$\mathrm{e}^{-\mathrm{i}(N-n)\phi}$ leads to

$\langle
f_{n,N}\rangle(\mathbf{R}_{\phi}\bm{r}_{1},\lambda_{1})\mathrm{e}^{-\mathrm{i}(N-n)\phi}=\mathrm{T}_{n}(\lambda_{1})\displaystyle\mathrm{V}_{N-n}(k\bm{r}_{1})+\\\\[6.0pt]
\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\mathrm{e}^{-\mathrm{i}(N-n^{\prime})\phi}\langle
f_{n^{\prime},N}\rangle(\mathbf{R}_{\phi}\bm{r}_{2},\lambda_{2})\mathrm{g}(|\bm{r}_{1}-\bm{r}_{2}|)\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2}.$
(4.5)

Now note that both $\langle
f_{n,N}\rangle(\mathbf{R}_{\phi}\bm{r}_{1},\lambda_{1})\mathrm{e}^{-\mathrm{i}(N-n)\phi}$
and $\langle f_{n,N}\rangle(\bm{r}_{1})$ solve exactly the same integral
equation. So by assuming uniqueness, i.e. that there is only one solution to
the above, we conclude that

$\langle f_{n,N}\rangle(\bm{r}_{1})=\langle
f_{n,N}\rangle(\mathbf{R}_{\phi}\bm{r}_{1})\mathrm{e}^{-\mathrm{i}(N-n)\phi},$
(4.6)

for any $\bm{r}_{1}$ and $\phi$. Let $(r_{1},\theta_{1})$ be the polar
coordinates of $\bm{r}_{1}$, then, without lose of generality, we then choose
$\phi=-\theta_{1}$ which leads to

$\langle f_{n,N}\rangle(r_{1},\theta_{1},\lambda_{1})=\langle
f_{n,N}\rangle(r_{1},0,\lambda_{1})\mathrm{e}^{\mathrm{i}(N-n)\theta_{1}}.$
(4.7)

Finally, because $\langle f_{n,N}\rangle(\bm{r}_{1})$ satisfies a wave
equation (3.1) it can be written in a series of the form (3.4). Using the
symmetry above we see that only one term in this series remains shown by
(4.3). ∎

###  4.2 The modal dispersion equation

We can deduce a simpler effective eigensystem and dispersion equation by using
the symmetry (4.3).

To start we substitute the expansions (C.8) and (4.3) into the modal
decomposition (4.2) to obtain

$\sum_{n_{1}}F_{nn_{1}}(\lambda_{1})\mathrm{V}_{n_{1}}(k_{\star}\bm{r}_{1})=\sum_{N_{1}}\alpha_{N_{1}}g_{N_{1}}F_{n,{N_{1}}}(\lambda_{1})\mathrm{V}_{{N_{1}}-n}(k_{\star}\bm{r}_{1})$
(4.8)

Then since $\mathrm{V}_{n}$ form an orthogonal basis of functions,

$F_{nn_{1}}=\sum_{N_{1}}\delta_{{N_{1}}-n_{1}-n}\alpha_{N_{1}}g_{N_{1}}F_{n,{N_{1}}}.$
(4.9)

To simplify the effective eigensystem (3.5) we consider one mode at a time by
taking $g_{N_{1}}=\delta_{{N_{1}}-N}$, which used in (4.9) implies that we can
substitute $F_{nn_{1}}=\delta_{N-n_{1}-n}\alpha_{N}F_{n,N}$ into (3.5). The
result after some algebraic manipulations is

$F_{n,N}(\lambda_{1})+\sum_{n^{\prime}}\mathrm{T}_{n}(\lambda_{1})\int_{\cal
S}{\cal
N}_{n^{\prime}-n}^{12}(k,k_{\star})F_{n^{\prime},N}(\lambda_{2})\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0.$
(4.10)

The above equation is identical to the case of the eigensystem for plane waves
[12], and matches also the eigensystems for a single type of particle [11, 19]
when taking $\mathfrak{n}(\lambda)=\delta(\lambda-\lambda_{1})$. This result
is somewhat expected as the ensemble wave equation (3.2) does not depend on
the incident wave and material geometry, which also explains why the modal
index $N$ only appears in $F_{n,N}$ in the above.

Instead of solving (3.5), it is far simpler to solve the above, and then write
the general solution in the modal form (4.9). In practice, to solve (4.10), we
can discretise the integral over ${\cal S}$ as a set of reals
$\\{t_{1},\dots,t_{S}\\}$. Then define a block vector $\mathbf{F}$ containing
the entries $F_{n}(t_{s})$ for $n=-M,-M+1,\ldots,M-1,M$, for some finite $M$,
and for $s=1,\dots,S$, so that the eigensystem becomes

$(\mathbf{I}+\mathbf{M})\cdot\mathbf{F}=\bm{0}$ (4.11)

where

$M_{nn^{\prime},ss^{\prime}}(k_{\star})=\mathrm{T}_{n}(t_{s})\,{\cal
N}_{n^{\prime}-n}^{12}(k,k_{\star})\,\mathfrak{n}(t_{s^{\prime}}).$ (4.12)

The parameter $k_{\star}$ is then obtained by solving the equation

$\mathrm{det}[\mathbf{I}+\mathbf{M}(k_{\star})]=0.$ (4.13)

1\. Multiple wavenumbers The dispersion equation (4.10) has infinitely many
eigenvalues $k_{p}$ with $p=1,2,\ldots$. Consequently, $\langle
f_{n}\rangle(\bm{r}_{1})$ can more generally be written as a sum over all the
eigenvalues $k_{p}$ and their corresponding eigenfunctions [8], which leads to
more accurate solutions. However, only a small difference in comparison to
using just the eigenvalue $k_{\star}$ with the smallest imaginary part is
observed (cf. [11, 8, 10] for details). For this reason, and for simplicity,
we only account for the one wavenumber $k_{\star}$ in this paper. See Figure 5
for a typical distribution of the many eigenvalues of (4.10).

###  4.3 The modal ensemble boundary condition

To determine the $\alpha_{N}$ that appear in (4.9) we need to use the ensemble
boundary condition. The simplest way to do this is again to take
$g_{N_{1}}=\delta_{N_{1}-N}$ in (4.9) and then substitute the result into
(3.8) to obtain

$\displaystyle
1+\frac{2\pi\alpha_{N}}{k_{\star}^{2}-k^{2}}\sum_{n^{\prime}}\int_{\cal
S}F_{n^{\prime},N}(\lambda_{2})N_{N-n^{\prime}}(kR_{2},k_{\star}R_{2})\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0,$
(4.14)

which we can use to determine:

$\alpha_{N}=-\frac{k_{\star}^{2}-k^{2}}{2\pi}\left(\sum_{n^{\prime}}\int_{\cal
S}F_{n^{\prime},N}(\lambda_{2})N_{N-n^{\prime}}(kR_{2},k_{\star}R_{2})\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}\right)^{-1}.$
(4.15)

###  4.4 The effective T-matrix

As discussed Section 4, the effective T-matrix $\mathcal{T}_{n,N}$ can easily
describe the average scattered wave for any incident wave through:

$\langle\mathfrak{F}_{n}\rangle=\sum_{N=-\infty}^{+\infty}{\cal
T}_{n,N}g_{N},$ (4.16)

where the $\langle\mathfrak{F}_{n}\rangle$, given by (2.16), are the average
coefficients of the waves scattered from the whole cylinder ${\cal R}$. Note
the above holds for any choice of $g_{N}$ and ${\cal T}_{n,N}$ does not depend
on $g_{N}$.

To calculate $\mathcal{T}_{n,N}$ we substitute the modal decomposition (4.2)
into (2.16) to obtain

$\langle\mathfrak{F}_{n}\rangle=\sum_{N}g_{N}\int_{\cal
S}\mathfrak{n}(\lambda_{1})\int_{{\cal
R}_{1}}\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(-k\bm{r}_{1})\langle
f_{n^{\prime},N}\rangle(\bm{r}_{1},\lambda_{1})\,\mathrm{d}\bm{r}_{1}\mathrm{d}\lambda_{1}.$
(4.17)

Comparing the above with the definition of the effective T-matrix (4.16), it
is clear that

${\cal T}_{n,N}=\int_{\cal S}\mathfrak{n}(\lambda_{1})\int_{{\cal
R}_{1}}\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(-k\bm{r}_{1})\langle
f_{n^{\prime},N}\rangle(\bm{r}_{1},\lambda_{1})\,\mathrm{d}\bm{r}_{1}\mathrm{d}\lambda_{1}.$
(4.18)

To calculate the above, we follow the same steps shown in Appendix C.
Specifically, use use Green’s second identity (C.4), the regular expansion
(4.3), and the orthogonality of the $\mathrm{V}_{n}$ functions to conclude
that

${\cal
T}_{n,N}=\delta_{N-n}\frac{2\pi\alpha_{n}}{k_{\star}^{2}-k^{2}}\int_{\cal
S}\mathfrak{n}(\lambda_{1})\sum_{n^{\prime}}F_{n^{\prime},n}(\lambda_{1}){\cal
Q}_{n-n^{\prime}}(kR_{1},k_{\star}R_{1})\mathrm{d}\lambda_{1}$ (4.19)

where $R_{1}$ is the radius of the disc ${\cal R}_{1}$ and we defined:

${\cal
Q}_{l}(x,y):=x\mathrm{J}_{l}^{\prime}(x)\mathrm{J}_{l}(y)-y\mathrm{J}_{l}(x)\mathrm{J}_{l}^{\prime}(y).$
(4.20)

Finally, we substitute $\alpha_{n}$ given by (4.15), which results in

${\cal T}_{n,N}=-\delta_{N-n}\frac{\displaystyle\int_{\cal
S}\mathfrak{n}(\lambda_{1})\sum_{n^{\prime}}F_{n^{\prime},n}(\lambda_{1}){\cal
Q}_{n-n^{\prime}}(kR_{1},k_{\star}R_{1})\mathrm{d}\lambda_{1}}{\displaystyle\int_{\cal
S}\mathfrak{n}(\lambda_{1})\sum_{n^{\prime}}F_{n^{\prime},n}(\lambda_{1})N_{n-n^{\prime}}(kR_{1},k_{\star}R_{1})\,\mathrm{d}\lambda_{1}}\cdot$
(4.21)

Note that once the effective T-matrix is known, the scattering from any
incident field can be computed with (2.15) after decomposing the incident in
modes (2.3a). For example, in Figure 6, we have plotted the total pressure
field $u:=u_{\mathrm{inc}}+u_{\mathrm{sc}}$ resulting from an incident plane
wave and a point source.

Average field for two different sources
Figure 6: Average pressure field when the incident field is a plane wave
(left) and a point source (right). The material is made of sound-hard
particles of radius one confined in a disc of radius 10.

###  4.5 Monopole particles only

The effective T-matrix (4.21) resembles the T-matrix for a homogeneous
cylinder, see for example (2.5). In fact, it is a weighted average of the
factors of a homogeneous T-matrix, as explained in the introduction. From this
observation we see that if the particles are monopole scatterers we obtain a
significant simplification.

Let us assume here that the particles scatter only monopole waves, in which
case the scattered field (2.3b) becomes

$u_{\mathrm{sc}}(\bm{r})=\sum_{i=1}^{J}f_{0}^{i}\mathrm{H}_{0}(kr-kr_{i}),$
(4.22)

which leads to $\langle f_{n,N}\rangle(\bm{r}_{1})=0$ if $n\neq 0$ and, as a
result of (4.3), $F_{n,N}=0$ if $n\neq 0$. Substituting this result into the
effective T-matrix (4.21), and assuming that the radius $R_{1}$ of the region
${\cal R}_{1}$ is the same for every type of particle $\lambda_{1}$, we
obtain:

$\displaystyle{\cal T}_{n}^{\mathrm{M}}=-\frac{\displaystyle
k\mathrm{J}^{\prime}_{n}(k\tilde{R})\mathrm{J}_{n}(k_{\star}\tilde{R})-k_{\star}\mathrm{J}_{n}(k\tilde{R})\mathrm{J}^{\prime}_{n}(k_{\star}\tilde{R})}{\displaystyle
k\mathrm{H}^{\prime}_{n}(k\tilde{R})\mathrm{J}_{n}(k_{\star}\tilde{R})-k_{\star}\mathrm{H}_{n}(k\tilde{R})\mathrm{J}^{\prime}_{n}(k_{\star}\tilde{R})}\quad\text{(monopole
scatterers).}$ (4.23)

This corresponds to the T-matrix of a homogeneous cylinder of radius
$\tilde{R}=R-a$, sound speed $c_{\star}=\omega/k_{\star}$ and mass density
$\rho_{\star}=\rho$ where $\rho$ is the density of the host medium (cf.
(2.5)).

##  5 Numerical results

The numerical results obtained in this paper use the open-source package
EffectiveTMatrix.jl written in Julia [23]. More details on the simulations
below can be found on the GitHub page https://github.com/Kevish-
Napal/EffectiveTMatrix.jl.

###  5.1 Optimized Monte Carlo simulations

We use a Monte Carlo method to validate our theoretical results (4.21) and
(4.23). To develop an efficient Monte Carlo method we rely on the following
symmetry of the modes:

$u_{\mathrm{inc}}(\bm{r})=\mathrm{J}_{N}(kr)\mathrm{e}^{\mathrm{i}N\theta}\implies\langle
u_{\mathrm{sc}}\rangle(\bm{r})=\mathrm{T}_{N}\mathrm{H}_{N}(kr)\mathrm{e}^{\mathrm{i}N\theta}.$
(5.1)

This result is easily obtained from (2.15) with the specific choice
$g_{n}=\delta_{n,N}$, which substituted into (2.14) leads

$\mathrm{T}_{N}=\langle\mathfrak{F}_{N}\rangle.$ (5.2)

In other words, $\mathrm{T}_{N}$ can be numerically estimated by simulating
the waves scattered from one particle configuration at a time by using (2.4),
and then taking the average of $\mathfrak{F}_{N}$ (2.12b) over many different
particles configurations.

To illustrate the efficiency of this Monte Carlo method we compare it with
another method which directly computes $\langle u_{\mathrm{sc}}\rangle$, which
has been commonly used in the literature [16]. For this second method we use
(2.3b) to compute $\langle u_{\mathrm{sc}}\rangle(R,0)$. Then from (5.1) we
can also compute $\mathrm{T}_{N}$ with

$\mathrm{T}_{N}=\mathrm{H}_{N}^{-1}(kR)\langle u_{\mathrm{sc}}\rangle(R,0).$
(5.3)

The two methods (5.2) and (5.3) are compared in Figure 7. The standard
deviation of the mean of the second method is larger than the first one,
resulting in a slower convergence. The reason is that (5.3), in contrary to
(5.2), includes all the modes of each scattered field computed for a specific
particles configuration:

$\langle
u_{\mathrm{sc}}\rangle(R,0)=\bigl{\langle}\sum_{n}\mathfrak{F}_{n}\mathrm{H}_{n}(kR)\bigr{\rangle}.$
(5.4)

While the terms $n\neq N$ of the sum vanish on average, they significantly
contribute to the standard deviation of the mean in (5.3).

Optimized Monte Carlo simulations
Figure 7: Computation of $\mathrm{T}_{0}$ using the two different methods
(5.2) (modal) and (5.3) (naive). We computed 5000 realisations of these
quantities for different configurations of particles, the plots on the top and
bottom respectively correspond to the real and imaginary parts of the results.
The distribution of the results is reported on the histograms on the left. The
plots on the right correspond to the cumulative average of the realisations.
Both methods converge to the same limit, however, the modal method converges
faster and presents a lower standard deviation of the mean.

###  5.2 Validation of the Effective Waves Method

We validate the Effective Waves method (4.18) against the Monte Carlo method
(5.2) for several frequencies $\omega$ ranging in $\Omega:=[0.05:0.015:1.5]$.
To this end we define the relative error, averaged over frequencies:

$\epsilon_{n}:=\frac{1}{|\Omega|}\sum_{\omega\in\Omega}\frac{|\mathrm{T}_{n}^{\text{MC}}(\omega)-\mathrm{T}_{n}^{\text{EWM}}(\omega)|}{|\mathrm{T}_{n}^{\text{MC}}(\omega)|}$
(5.5)

where $\mathrm{T}_{n}^{\text{MC}}(\omega)$ is obtained following the Monte
Carlo method (5.2) and $\mathrm{T}_{n}^{\text{EWM}}(\omega)$ following the
Effective Waves method (4.18). A few plots of
$\mathrm{T}_{n}^{\text{MC}}(\omega)$ and $\mathrm{T}_{n}^{\text{EWM}}(\omega)$
are provided in Figures 2 and 8. The values of $\epsilon_{0}$, $\epsilon_{1}$,
$\epsilon_{2}$, $\epsilon_{3}$, $\epsilon_{4}$ for the cases of sound soft and
sound hard particles are reported in Table 1.

Monte Carlo results: sound-soft particles
Figure 8: Compares various methods to calculate the components
$\mathrm{T}_{1}$, $\mathrm{T}_{2}$, $\mathrm{T}_{3}$ and $\mathrm{T}_{4}$ of
the T-matrix of a cylinder filled with sound soft particles. The solid red
line is our effective waves method (EWM) (4.21), the black points are from the
Monte Carlo method (MC) (5.2), and the dashed blue line is our method when
only monopole scattering is accounted for (EWM-MA) (4.23). The general
expression of the effective T-matrix matches the Monte Carlo results. The EWM-
MA method only matches well with the Monte Carlo for low frequencies.

| $\epsilon_{0}$ | $\epsilon_{1}$ | $\epsilon_{2}$ | $\epsilon_{3}$ | $\epsilon_{4}$
---|---|---|---|---|---
sound soft particles | $6.75\mathrm{e}^{-2}$ | $4.59\mathrm{e}^{-2}$ | $5.40\mathrm{e}^{-2}$ | $6.45\mathrm{e}^{-2}$ | $7.45\mathrm{e}^{-2}$
sound hard particles | $2.66\mathrm{e}^{-2}$ | $2.43\mathrm{e}^{-2}$ | $2.33\mathrm{e}^{-2}$ | $2.35\mathrm{e}^{-2}$ | $2.37\mathrm{e}^{-2}$
Table 1: Relative errors $\epsilon_{n}$ defined by (5.5) in the cases of sound
soft and sound hard. The particles are of radius 1 and confined in a circular
area of radius 20. The volume fraction is 0.05 which corresponds to 19.7
particles on average.

##  6 Conclusion

Main goal. Our main goal was to describe how an incident wave is scattered
from a cylinder filled with smaller cylinders, which we have called particles,
that are placed in a disordered but correlated way. We describe this
correlation through the inter-particle pair correlation, see (2.8). The
literature so far has focused on plane waves scattered from a halfspace or
plate filled with a particulate [14]. There has been at least one paper on
solving this scenario, but used an ad hoc method, whereas here everything is
deduced from first principles making only two assumptions: such as the Quasi-
Crystalline Approximation (QCA) [17], and expressing the average field as a
sum of effective waves, which has been shown to be the analytic solution [8].
One of the key advantages of describing the scattering from a cylinder with 2D
particles is that it is far easier to validate this scenario with direct
numerical simulations. Validation is still necessary as the theory requires
the use of QCA whose validity is not clearly established. Further, validating
this scenario also serves to validate the predicted effective wavenumbers from
any material geometry [10].

Modal scattering. A method that allowed us to greatly simplify both the
theoretical and Monte Carlo calculations which we use to validate this work,
was to make use of all symmetries present. We achieved this by solving for
each polar mode of the incident wave separately. This simple, but effective
technique, leads us to an effective T-matrix given by equation (4.21) that can
be used to calculate the average scattered wave from any incident wave, see
section 3 for a brief overview. We note that we were able to describe the
scattered field without calculating the average transmitted field. In future
work, this may be interesting to do, for example, to clearly identify when
different effective wavenumbers are excited [14].

Effective T-matrix. The result of our theoretical work is summarised by the
T-matrix (4.21). Beyond using just QCA, to calculate this T-matrix, we also
assumed that only one effective wavenumber $k_{\star}$ is excited. This is
true for a wide range of parameters but it is not always the case. In
particular it appears that very strong scattering at moderate frequencies can
trigger more than one effective wavenumber to be excited [14, 11, 10]. One
possible extension to our work is to include the effect of more than one
effective wavenumber.

Monopole Scatterers. One surprising result, is that if the particles only
scatter monopole waves, then the effective T-matrix greatly simplifies and
becomes (4.23). This form is exactly the same as the T-matrix for a homogenous
cylinder, one with constant material parameters. We hypothesise that any
material filled with monopole scatterers would, on average, respond like a
homogeneous material. Monopole scatterers are a good approximation for many
types of resonant particles [4]. Figure 8 compares the results for the
monopole scattering approximation with Monte Carlo results for sound-soft
particles, which does not assume the particles scatter like monopoles.

Monte Carlo. Beyond deducing an effective T-matrix for a particulate cylinder,
we also developed an efficient Monte-Carlo method, which matched our
theoretical predictions very accurately see Figures 2 and 8. Our numerical
validation was for a broad frequency range, but we did not cover a broad range
of parameters. Doing this would be valuable future work, and could help
clearly identify the limits of QCA and different approximations for the inter-
particle pair-correlation.

A prototype for new materials. The setting we deduce is the ideal case to test
new disordered particulates. That is, to use exotic inter-particle pair-
correlation to achieve effects such as band-gaps, or impedance matching. The
scenario of a cylinder filled with 2D particles is ideal for testing new types
of particles, and inter-particle pair-correlations, because they can be easily
validated with Monte-Carlo methods. When designing new materials with exotic
responses, and stretching the limits of the theory, we need to have a way to
validate those predictions. And this paper provides that.

##  Appendix A Bessel functions and translation matrices

Given two points $\bm{x},\,\bm{y}\in\mathbb{R}^{2}$, we have the following
identities where $\bm{d}=\bm{x}-\bm{y}$

$\begin{array}[]{|lrcll}\text{i)}&\mathrm{V}_{n}(\bm{y})&=&\displaystyle\sum_{n^{\prime}=-\infty}^{+\infty}\mathrm{V}_{n-n^{\prime}}(\bm{d})\mathrm{V}_{n^{\prime}}(\bm{x}),&\text{for
all }\bm{x},\bm{d}\in\mathbb{R}^{2}\\\
\text{ii)}&\mathrm{U}_{n}(\bm{y})&=&\displaystyle\sum_{n^{\prime}=-\infty}^{+\infty}\mathrm{V}_{n-n^{\prime}}(\bm{d})\mathrm{U}_{n^{\prime}}(\bm{x}),&\text{for
all }|\bm{x}|>|\bm{d}|\\\
\text{iii)}&\mathrm{U}_{n}(\bm{y})&=&\displaystyle\sum_{n^{\prime}=-\infty}^{+\infty}\mathrm{U}_{n-n^{\prime}}(\bm{d})\mathrm{V}_{n^{\prime}}(\bm{x}),&\text{for
all }|\bm{x}|<|\bm{d}|.\end{array}$ (A.1)

The above formulas are direct consequences of Graf’s theorem (see [21, Th.
2.11-2.12] for instance).

##  Appendix B Ensemble averaging

Here we give a brief overview of ensemble averaging so that this paper is more
self-contained. For more details, see [10, 21] and the references within. To
simplify computations, we represent one particle configurations with:

$\displaystyle\Lambda=\bm{r}_{1},\lambda_{1}\dots\bm{r}_{J},\lambda_{J},\quad\quad\Lambda^{(1)}=\bm{r}_{2},\lambda_{2}\dots\bm{r}_{J},\lambda_{J},$
$\displaystyle\Lambda^{(1,j)}=\bm{r}_{2},\lambda_{2},\dots,\bm{r}_{j-1},\lambda_{j-1},\bm{r}_{j+1},\lambda_{j+1},\dots,\bm{r}_{J},\lambda_{J}.$

Using these definitions, we define the ensemble average, and conditional
ensemble averages, of a quantity $\mathrm{A}$, which can depend on the
positions and properties of all the particles, as

$\displaystyle\langle\mathrm{A}\rangle:=\int\mathrm{A}(\Lambda)\,\mathrm{p}(\Lambda)\mathrm{d}\Lambda\qquad\langle\mathrm{A}\rangle(\bm{r}_{i},\lambda_{i}):=\int\mathrm{A}(\Lambda)\,\mathrm{p}(\Lambda^{(i)}|\bm{r}_{i}\lambda_{i})\mathrm{d}\Lambda^{(i)},$
(B.1)
$\displaystyle\langle\mathrm{A}\rangle(\bm{r}_{i},\lambda_{i};\bm{r}_{j},\lambda_{j}):=\int\mathrm{A}(\Lambda)\,\mathrm{p}(\Lambda^{(i,j)}|\bm{r}_{i},\lambda_{i};\bm{r}_{j},\lambda_{j})\mathrm{d}\Lambda^{(i,j)},$
(B.2)

where the domain of integration for $\Lambda$ is over all possible particles
positions and properties. The term
$\langle\mathrm{A}\rangle(\bm{r}_{i},\lambda_{i};\bm{r}_{j},\lambda_{j})$ is
the ensemble average of $\mathrm{A}$ conditional to
$(\bm{r}_{1},\lambda_{1},\bm{r}_{2},\lambda_{2})$.

Taking an ensemble average (B.1) on both sides of (2.12a) leads to

$\displaystyle\langle\mathfrak{F}_{n}\rangle$
$\displaystyle=\sum_{i=1}^{J}\sum_{n^{\prime}=-\infty}^{+\infty}\int(-1)^{n-n^{\prime}}\mathrm{V}_{n-n^{\prime}}(k\bm{r}_{i})\langle
f_{n^{\prime}}^{i}\rangle(\bm{r}_{i},\lambda_{i})\mathrm{p}(\bm{r}_{i},\lambda_{i})\,\mathrm{d}\bm{r}_{i}\mathrm{d}\lambda_{i}$
(B.3)
$\displaystyle=J\sum_{n^{\prime}=-\infty}^{+\infty}\int(-1)^{n-n^{\prime}}\mathrm{V}_{n-n^{\prime}}(k\bm{r}_{i})\langle
f_{n^{\prime}}^{1}\rangle(\bm{r}_{1},\lambda_{1})\mathrm{p}(\bm{r}_{1},\lambda_{1})\,\mathrm{d}\bm{r}_{1}\mathrm{d}\lambda_{1}$

where we used the definition of conditional probability:
$\mathrm{p}(\Lambda)=\mathrm{p}(\bm{r}_{i},\lambda_{i})\mathrm{p}(\Lambda^{(i)})$,
the definition of conditional average (B.2) to introduce the term $\langle
f_{n^{\prime}}^{i}\rangle(\bm{r}_{i},\lambda_{i})$, and then used that
particles are indistinguishable333Said in another way, the variables of
integration $\bm{r}_{i}$ and $\lambda_{i}$ are just dummy variables which can
be all changed to $\bm{r}_{1}$ and $\lambda_{1}$.. Finally, using (2.7) and
(2.9) leads to formula (2.16).

###  B.1 Average governing equation

Here we briefly show how to reach the averaging governing equation by using
just one assumption, the Quasi-Crystalline Approximation (QCA). For more
details see [10].

Taking the conditional average (B.1)2 of (2.4) with $i=1$ gives

$\langle
f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{1})g_{n^{\prime}}\\\
+\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{j\neq
1}\sum_{n^{\prime}}\int\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{j})f_{n^{\prime}}^{j}(\bm{r}_{1},\lambda_{1})\mathrm{p}(\Lambda^{(1)}|\bm{r}_{1},\lambda_{1})\mathrm{d}\Lambda^{(1)}.$
(B.4)

We can simplify the above by using the definition of the conditional average
(B.2) of $f_{n^{\prime}}^{j}(\bm{r}_{1},\lambda_{1};\bm{r}_{j},\lambda_{j})$,
using (2.10), and that particles are indistinguishable to obtain

$\langle
f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{1})g_{n^{\prime}}\\\
+\mathrm{T}_{n}(\lambda_{1})\displaystyle\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime}}^{2}\rangle(\bm{r}_{2},\lambda_{2},\bm{r}_{1},\lambda_{1})\mathrm{g}(\bm{r}_{1},\lambda_{1};\bm{r}_{2},\lambda_{2})\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2}.$
(B.5)

However, this equation is not a closed form equation for $\langle
f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})$ and an extra assumption is required
to proceed further.

A standard solution found in the literature to tackle the problem mentioned
above is to use the quasi-crystalline approximation (QCA). It is stated as
follows:

$\langle
f_{n}^{2}\rangle(\bm{r}_{1},\lambda_{1};\bm{r}_{2},\lambda_{2})\approx\langle
f_{n}^{2}\rangle(\bm{r}_{2},\lambda_{2}),\quad|\bm{r}_{1}-\bm{r}_{2}|\geq
a_{12}\quad\text{(QCA).}$ (B.6)

See [12] for a brief discussion on this approximation.

Finally we use that particles are indistinguishable which implies that
$\langle f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})=\langle
f_{n}^{2}\rangle(\bm{r}_{2},\lambda_{2})$ when $\bm{r}_{1}=\bm{r}_{2}$ and
$\lambda_{1}=\lambda_{2}$ to substitute

$\langle f_{n}^{2}\rangle(\bm{r}_{2},\lambda_{2})=\langle
f_{n}\rangle(\bm{r}_{2},\lambda_{2})\quad\text{and}\quad\langle
f_{n}^{1}\rangle(\bm{r}_{1},\lambda_{1})=\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ (B.7)

into (B.5), which together with QCA (B.6) leads to the average governing
equation (2.17).

##  Appendix C The effective waves method

###  C.1 Derivation of two ensemble equations

Here we show how to use the effective wave assumption to rewrite the governing
equation (4.1) into two separate equations: the effective wave equation and
the effective boundary conditions. To achieve this we define using set builder
notation:

$D(\bm{y},a):=\\{\bm{x}\in\mathcal{R}:\;|\bm{x}-\bm{y}|\leq a\\}.$

Using the decomposition of the pair correlation function (2.11), we split the
domain of integration in the governing equation (2.17) into two integrals: one
over $D(\bm{r}_{1},a_{12})$, another one over ${\cal R}_{2}\setminus
D(\bm{r}_{1},a_{12})$ in the form

$\displaystyle\langle$ $\displaystyle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\sum_{n^{\prime}}\mathrm{V}_{N-n}(k\bm{r}_{1})g_{n^{\prime}}+$
(C.1) $\displaystyle\mathrm{T}_{n}(\lambda_{1})\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{{\cal R}_{2}\setminus
D(\bm{r}_{1},a_{12})}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime}}\rangle(\bm{r}_{2},\lambda_{2})\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2}+$
$\displaystyle\mathrm{T}_{n}(\lambda_{1})\sum_{n^{\prime}}\int_{\cal
S}\mathfrak{n}(\lambda_{2})\int_{D(\bm{r}_{1},a_{12},b_{12})}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime}}\rangle(\bm{r}_{2},\lambda_{2})\delta\mathrm{g}(|\bm{r}_{1}-\bm{r}_{2}|,\lambda_{1},\lambda_{2})\mathrm{d}\bm{r}_{2}\mathrm{d}\lambda_{2}.$

where $D(\bm{r}_{1},a_{12},b_{12}):=D(\bm{r}_{1},b_{12})\setminus
D(\bm{r}_{1},a_{12})$ and the annulus $D(\bm{r}_{1},a_{12},b_{12})$ is
completely contain within ${\cal R}_{2}$ when

$\mathrm{dist}(\bm{r}_{1},\partial{\cal R}_{2})\geq b_{12}.$ (C.2)

In this section, and in the whole paper, we only solve (4.1) for $\bm{r}_{1}$
that satisfies the above. This avoids the boundary layer [11, 8] which greatly
complicates the solution and is only needed when there is a large particle
volume fraction, moderate frequencies, and strongly scattering particles.

The last integral in (C.1) can be simplified by changing the variable of
integration to $\bm{r}=\bm{r}_{2}-\bm{r}_{1}$ which leads to the integral

${\cal
K}_{n^{\prime}n}(\bm{r}_{1},\lambda_{2}):=\int_{D(\bm{0},a_{12},b_{12})}\mathrm{U}_{n^{\prime}-n}(-k\bm{r})\langle
f_{n^{\prime}}\rangle(\bm{r}+\bm{r}_{1},\lambda_{2})\delta\mathrm{g}(r,\lambda_{1},\lambda_{2})\,\mathrm{d}\bm{r}.$
(C.3)

The first integral over $\bm{r}_{2}$ in (C.1) can be further simplified by
using Green’s theorem to replace the volume integral over ${\cal
R}_{2}\setminus D(\bm{r}_{1},a_{12})$ by surface integrals: given any two
function smooth functions $u,v$ which satisfy

$\Delta u(\bm{r})+k_{\star}u(\bm{r})=0\quad\text{and}\quad\Delta
v(\bm{r})+kv(\bm{r})=0,$

over a set $\Omega$ we have that

$(k^{2}-k_{\star}^{2})\int_{\Omega}uv\,\mathrm{d}\bm{r}=\int_{\Omega}\left(\Delta
uv-u\Delta
v\right)\,\mathrm{d}\bm{r}=\int_{\partial\Omega}\left(\partial_{\bm{\nu}}uv-u\partial_{\bm{\nu}}v\right)\,\mathrm{d}s(\bm{r}).$
(C.4)

With $u(\bm{r}_{2})$ substituted for $\langle
f_{n^{\prime},N}\rangle(\bm{r}_{2},\lambda_{2})$ and $v(\bm{r}_{2})$
substituted for $\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})$ we can
use the above to deduce:

$\int_{{\cal R}_{2}\setminus
D(\bm{r}_{1},a_{12})}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\langle
f_{n^{\prime}}\rangle(\bm{r}_{2},\lambda_{2})\mathrm{d}\bm{r}_{2}=\frac{{\cal
I}_{n^{\prime}n}(\bm{r}_{1})-{\cal
J}_{n^{\prime}n}(\bm{r}_{1})}{k^{2}-k_{\star}^{2}},$ (C.5)

where we defined

$\displaystyle{\cal I}_{n^{\prime}n}(\bm{r}_{1})$
$\displaystyle:=\int_{\partial{\cal
R}_{2}}\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})\frac{\partial\langle
f_{n^{\prime}}\rangle(\bm{r}_{2},\lambda_{2})}{\partial{\bm{\nu}}_{2}}-\frac{\partial\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})}{\partial{\bm{\nu}}_{2}}\langle
f_{n^{\prime}}\rangle(\bm{r}_{2},\lambda_{2})\mathrm{d}A_{2},$ (C.6)
$\displaystyle{\cal J}_{n^{\prime}n}(\bm{r}_{1})$
$\displaystyle:=\int_{\partial
D(\bm{0},a_{12})}\mathrm{U}_{n^{\prime}-n}(-k\bm{r})\frac{\partial\langle
f_{n^{\prime}}\rangle(\bm{r}+\bm{r}_{1},\lambda_{2})}{\partial{\bm{\nu}}}-\frac{\partial\mathrm{U}_{n^{\prime}-n}(-k\bm{r})}{\partial{\bm{\nu}}}\langle
f_{n^{\prime}}\rangle(\bm{r}+\bm{r}_{1},\lambda_{2})\mathrm{d}A.$

Finally, substituting (C.3) and (C.5) into the governing equation (2.17) gives

$\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\mathrm{T}_{n}(\lambda_{1})\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{1})g_{n^{\prime}}\\\
+\sum_{n^{\prime}}\mathrm{T}_{n}(\lambda_{1})\int_{\cal S}\left[\frac{{\cal
I}_{n^{\prime}n}(\bm{r}_{1})-{\cal
J}_{n^{\prime}n}(\bm{r}_{1})}{k^{2}-k_{\star}^{2}}+{\cal
K}_{n^{\prime}n}(\bm{r}_{1})\right]\mathfrak{n}(\lambda_{2})\mathrm{d}\lambda_{2}.$
(C.7)

The above now can be split into two separate equations by noting that the
functions $\langle f_{n}\rangle(\bm{r}_{1})$, ${\cal
J}_{n^{\prime}n}(\bm{r}_{1})$ and ${\cal K}_{n^{\prime}n}(\bm{r}_{1})$ satisfy
the wave equation with wavenumber $k_{\star}$, while
$\mathrm{V}_{n}(k\bm{r}_{1})$ and ${\cal I}_{n^{\prime}n}(\bm{r}_{1})$ satisfy
the wave equation with wavenumber $k$. Since solutions of the Helmholtz
equation with different wavenumbers are independent, see [10] for details,
(C.7) can be split into the ensemble wave equation (3.2) containing the terms
with wavenumber $k_{\star}$ and the ensemble boundary conditions (3.3)
containing the terms with wavenumber $k$.

###  C.2 The effective eigensystem

Here we deduce a general eigensystem which can be used to determine the
effective wavenumber $k_{\star}$ and write $\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ in terms of eigenfunctions.

Since $\langle f_{n}\rangle(\bm{r}_{1},\lambda_{1})$ satisfies the wave
equation (3.1), it can be decomposed into the modes

$\langle
f_{n}\rangle(\bm{r}_{1},\lambda_{1})=\sum_{n_{1}}F_{nn_{1}}(\lambda_{1})\mathrm{V}_{n_{1}}(k_{\star}\bm{r}_{1}),$
(C.8)

where $\mathrm{V}_{n_{1}}$ is defined in (2.2).

The unknowns $k_{\star}$ and $F_{nn_{1}}(\lambda_{1})$ can be determined by
substituting (C.8) into (3.2), which requires the term

$\displaystyle\langle f_{n}\rangle(\bm{r}+\bm{r}_{1},\lambda_{2})$
$\displaystyle=\sum_{n_{1}n_{2}}F_{nn_{1}}(\lambda_{2})\mathrm{V}_{n_{1}-n_{2}}(k_{\star}\bm{r})\mathrm{V}_{n_{2}}(k_{\star}\bm{r}_{1}),$
(C.9)

where the right side is a result of using Graf’s addition theorem (A.1, i) in
(3.4).

By substituting (C.9) in ${\cal K}_{n^{\prime}n}(\bm{r}_{1})$ (C.3) we can use
the orthogonality of the cylindrical Bessel functions to remove the sum over
$n_{2}$, because only the cases $(n_{1}-n_{2})=(n-n^{\prime})$ are non-zero.

Likewise, we can perform the same simplification by substituting (C.9) in
${\cal J}_{n^{\prime}n}(\bm{r}_{1})$ (C.6). The simplifications result in

$\displaystyle{\cal K}_{n^{\prime}n}(\bm{r}_{1})$
$\displaystyle=2\pi\mathrm{W}_{n^{\prime}-n}(k,k_{\star})\sum_{n_{1}}F_{n^{\prime}n_{1}}(\lambda_{2})\mathrm{V}_{n_{1}+n^{\prime}-n}(k_{\star}\bm{r}_{1}),$
(C.10) $\displaystyle{\cal J}_{n^{\prime}n}(\bm{r}_{1})$
$\displaystyle=-2\pi\mathrm{N}_{n^{\prime}-n}(ka_{12},k_{\star}a_{12})\sum_{n_{1}}F_{n^{\prime}n_{1}}(\lambda_{2})\mathrm{V}_{n_{1}+n^{\prime}-n}(k_{\star}\bm{r}_{1}),$

where we introduced the notations

$\displaystyle\mathrm{W}_{l}(k,k_{\star})$
$\displaystyle:=\int_{a_{12}}^{b_{12}}\mathrm{H}_{l}(kr)\mathrm{J}_{l}(k_{\star}r)\delta\mathrm{g}(r,\lambda_{1},\lambda_{2})r\,\mathrm{d}r,$
(C.11) $\displaystyle\mathrm{N}_{l}(x,y)$
$\displaystyle:=x\mathrm{H}^{\prime}_{l}(x)\mathrm{J}_{l}(y)-y\mathrm{H}_{l}(x)\mathrm{J}^{\prime}_{l}(y).$

Finally, substituting (C.10) in the ensemble wave equation (3.2), and again
using the orthogonality of the cylindrical Bessel functions, we reach equation
(3.5) where the following term appears:

${\cal
N}_{l}^{12}(k,k_{\star})=2\pi\frac{\mathrm{N}_{l}(ka_{12},k_{\star}a_{12})}{k_{\star}^{2}-k^{2}}-2\pi\mathrm{W}_{l}(k,k_{\star}).$
(C.12)

##  Appendix D Boundary condition for effective waves

The eigensystem (3.5) is not enough to fully determine the $F_{nn_{1}}$, to do
so we need to substitute (3.4) into the ensemble boundary condition (3.3). To
achieve this the first step is to use Graf’s addition theorem (A.1, iii) with
$\bm{x}=\bm{r}_{1}$ and $\bm{y}=-\bm{r}_{2}$ to obtain

$\mathrm{U}_{n^{\prime}-n}(k\bm{r}_{1}-k\bm{r}_{2})=\sum_{n_{3}}\mathrm{V}_{n_{3}}(k\bm{r}_{1})\mathrm{U}_{n^{\prime}-n-n_{3}}(-k\bm{r}_{2})=\sum_{n_{2}}\mathrm{V}_{n_{2}-n+n^{\prime}}(k\bm{r}_{1})\mathrm{U}_{-n_{2}}(-k\bm{r}_{2}),$
(D.1)

where we used the change of variable $n_{2}=n+n_{3}-n^{\prime}$, and that
$|\bm{r}_{1}|<|\bm{r}_{2}|$ because for ${\cal I}_{n^{\prime}n}$ the variable
$\bm{r}_{2}$ is on the boundary $\partial{\cal R}_{2}$, whereas $\bm{r}_{1}$
satisfies (C.2). Substituting (3.4) in ${\cal I}_{n^{\prime}n}$ defined by
(C.6) and using the above (D.1) gives

${\cal
I}_{n^{\prime}n}(\bm{r}_{1})=\sum_{n_{2}n_{1}}F_{n^{\prime}n_{1}}(\lambda_{2})\mathcal{B}_{n_{1}n_{2}}\mathrm{V}_{n_{2}-n+n^{\prime}}(k\bm{r}_{1})$
(D.2)

where

$\mathcal{B}_{n_{1}n_{2}}=(-1)^{n_{2}}\int_{\partial{\cal
R}_{2}}\left[\mathrm{U}_{-n_{2}}(k\bm{r}_{2})\frac{\partial\mathrm{V}_{n_{1}}(k_{\star}\bm{r}_{2})}{\partial{\bm{\nu}}_{2}}-\frac{\partial\mathrm{U}_{-n_{2}}(k\bm{r}_{2})}{\partial{\bm{\nu}}_{2}}\mathrm{V}_{n_{1}}(k_{\star}\bm{r}_{2})\right]\mathrm{d}A_{2}$
(D.3)

then substituting (D.2) in (3.3) leads to

$\sum_{n^{\prime}}\mathrm{V}_{n^{\prime}-n}(k\bm{r}_{1})g_{n^{\prime}}+\sum_{n^{\prime}n_{2}n_{1}}\int_{\cal
S}F_{n^{\prime}n_{1}}(\lambda_{2})\mathrm{V}_{n_{2}-n+n^{\prime}}(k\bm{r}_{1})\frac{\mathcal{B}_{n_{1}n_{2}}}{k^{2}-k_{\star}^{2}}\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0.$
(D.4)

We can further simplify the above by using the orthogonality of the functions
$\mathrm{V}_{n}$ to obtain

$\displaystyle g_{N}+\sum_{n^{\prime}n_{1}}\int_{\cal
S}F_{n^{\prime}n_{1}}(\lambda_{2})\frac{\mathcal{B}_{n_{1}(N-n^{\prime})}}{k^{2}-k_{\star}^{2}}\mathfrak{n}(\lambda_{2})\,\mathrm{d}\lambda_{2}=0,$
(D.5)

which holds for every $N$.

When all particles are in a disk, then $\partial{\cal R}_{2}$ is a circle and
the above simplifies. This is the only case we completely resolve in this
paper. Let $R_{2}$ be the radius of the disk ${\cal R}_{2}$, then
$n_{2}=n_{1}$ in (D.3) which reduces to

$\mathcal{B}_{n_{1}n_{2}}=-2\pi\delta_{n_{1}-n_{2}}N_{n_{1}}(kR_{2},k_{\star}R_{2}),$
(D.6)

where $N_{n_{1}}$ is defined by (C.11)2. Substituting this into (D.5) leads to
(3.8).

##  Appendix E Elementary proof that thee effective T-matrix is diagonal

We provide an elementary proof that ${\cal T}$ is diagonal when the particles
are confined in a disk of radius $\mathrm{R}$. To this end, we consider the
scattering from the modal source $u_{\mathrm{inc}}^{N}$ obtained for
$g_{n}=\delta_{n,N}$:

The notation $\mathfrak{F}_{n,N}$ is the corresponding $\mathfrak{F}_{n}$ to
the specific incident field with $g_{n}=\delta_{n,N}$ (compare with (2.12a)).
We then denote by $f_{n,N}^{i}(\sigma)$ the resulting solution of (2.4) for
the specific configuration $\sigma=\bm{r}_{1},\dots,\bm{r}_{J}$. The rotation
by angle $\phi$ of the particles $\bm{r}_{1},\dots,\bm{r}_{J}$ correspond to
another valid configuration (because the random material is cylindrical), for
which the solutions are given by

$f_{n,N}^{i}(\mathbf{R}_{\phi}\sigma)=\mathrm{e}^{\mathrm{i}(N-n)\phi}f_{n,N}^{i}(\sigma)$
(E.1)

This tells us how the rotation of a configuration modifies the coefficient
$\mathfrak{F}_{n,N}$, using (2.12a):

$\mathfrak{F}_{n,N}(\mathbf{R}_{\phi}\sigma)=\sum_{i=1}^{J}\sum_{n^{\prime}=-\infty}^{+\infty}\overline{\mathrm{V}_{n-n^{\prime}}(k\bm{r}_{i})}\mathrm{e}^{\mathrm{i}(N-n)\phi}f_{n^{\prime},N}^{i}(\sigma).$
(E.2)

Consequently,

$\langle\mathfrak{F}_{n,N}\rangle=\int_{\sigma}\mathfrak{F}_{n,N}(\sigma)\mathrm{p}(\sigma)\,\mathrm{d}\sigma=\int_{\sigma}\frac{1}{2\pi}\int_{0}^{2\pi}\mathfrak{F}_{n,N}(\mathbf{R}_{\phi}\sigma)\mathrm{p}(\sigma)\,\mathrm{d}\sigma\mathrm{d}\phi=\delta_{n,N}\int_{\sigma}\mathfrak{F}_{n,N}(\sigma)\mathrm{p}(\sigma)\,\mathrm{d}\sigma$
(E.3)

Finally, we deduce

$\displaystyle{\cal
T}_{n,N}=\delta_{n,N}\int\mathfrak{F}_{n,N}(\bm{r}_{1},\dots,\bm{r}_{J})\mathrm{p}(\bm{r}_{1},\dots,\bm{r}_{J})\,\mathrm{d}\bm{r}_{1},\dots,\mathrm{d}\bm{r}_{J}.$
(E.4)

This analysis proves that only the diagonal terms of the effective T-matrix
are nonzero and can be estimated by ${\cal
T}_{n,n}=\langle\mathfrak{F}_{n,n}\rangle$.

## References

  * [1] SK Bose and AK Mal.  Longitudinal shear waves in a fiber-reinforced composite.  International Journal of Solids and Structures, 9(9):1075–1085, 1973.
  * [2] RE Challis, MJW Povey, ML Mather, and AK Holmes.  Ultrasound techniques for characterizing colloidal dispersions.  Reports on progress in physics, 68(7):1541, 2005.
  * [3] Mathieu Chekroun, Loïc Le Marrec, Bruno Lombard, and Joël Piraux.  Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers.  Waves in Random and Complex media, 22(3):398–422, 2012.
  * [4] Philip A Cotterill, David Nigro, and William J Parnell.  Deeply subwavelength giant monopole elastodynamic metacluster resonators.  Proceedings of the Royal Society A, 478(2263):20220026, 2022.
  * [5] J Dubois, C Aristégui, O Poncelet, and A L Shuvalov.  Coherent acoustic response of a screen containing a random distribution of scatterers: Comparison between different approaches.  Journal of Physics: Conference Series, 269:012004, January 2011\.
  * [6] Leslie L Foldy.  The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers.  Physical review, 67(3-4):107, 1945.
  * [7] Mahadevan Ganesh and Stuart C Hawkins.  Algorithm 975: Tmatrom—a t-matrix reduced order model software.  ACM Transactions on Mathematical Software (TOMS), 44(1):1–18, 2017\.
  * [8] Artur L Gower, I David Abrahams, and William J Parnell.  A proof that multiple waves propagate in ensemble-averaged particulate materials.  Proceedings of the Royal Society A, 475(2229):20190344, 2019.
  * [9] Artur L Gower, Robert M Gower, Jonathan Deakin, William J Parnell, and I David Abrahams.  Characterising particulate random media from near-surface backscattering: A machine learning approach to predict particle size and concentration.  Europhysics Letters, 122(5):54001, 2018.
  * [10] Artur L Gower and Gerhard Kristensson.  Effective waves for random three-dimensional particulate materials.  New Journal of Physics, 23(6):063083, 2021.
  * [11] Artur L Gower, William J Parnell, and I David Abrahams.  Multiple waves propagate in random particulate materials.  SIAM Journal on Applied Mathematics, 79(6):2569–2592, 2019.
  * [12] Artur L Gower, Michael JA Smith, William J Parnell, and I David Abrahams.  Reflection from a multi-species material and its transmitted effective wavenumber.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2212):20170864, 2018.
  * [13] Nail A Gumerov and Ramani Duraiswami.  Computation of scattering from clusters of spheres using the fast multipole method.  The Journal of the Acoustical Society of America, 117(4):1744–1761, 2005.
  * [14] Aris Karnezis, Paulo S. Piva, and Art L. Gower.  The average transmitted wave in random particulate materials.  2023\.
  * [15] S Koc and Weng Cho Chew.  Calculation of acoustical scattering from a cluster of scatterers.  The Journal of the Acoustical Society of America, 103(2):721–734, 1998.
  * [16] Jin Au Kong, Leung Tsang, Kung-Hau Ding, and Chi On Ao.  Scattering of electromagnetic waves: numerical simulations.  John Wiley & Sons, 2004.
  * [17] Melvin Lax.  Multiple scattering of waves. ii. the effective field in dense systems.  Physical Review, 85(4):621, 1952.
  * [18] C. M. Linton and P. A. Martin.  Multiple scattering by random configurations of circular cylinders: Second-order corrections for the effective wavenumber.  J. Acoust. Soc. Am., 117(6):3413, 2005.
  * [19] CM Linton and PA Martin.  Multiple scattering by random configurations of circular cylinders: Second-order corrections for the effective wavenumber.  The Journal of the Acoustical Society of America, 117(6):3413–3423, 2005.
  * [20] Paul A Martin.  Acoustic scattering by inhomogeneous obstacles.  SIAM Journal on Applied Mathematics, 64(1):297–308, 2003.
  * [21] Paul A Martin.  Multiple scattering: interaction of time-harmonic waves with N obstacles.  Number 107. Cambridge University Press, 2006.
  * [22] Michael I Mishchenko, Larry D Travis, and Andrew A Lacis.  Multiple scattering of light by particles: radiative transfer and coherent backscattering.  Cambridge University Press, 2006.
  * [23] Kevish K. Napal.  Julia package EffectiveTMatrix.jl, August 2023.
  * [24] Yaxian Ni, Lei Gao, and Cheng-Wei Qiu.  Achieving invisibility of homogeneous cylindrically anisotropic cylinders.  Plasmonics, 5:251–258, 2010.
  * [25] John B Pendry and David R Smith.  Reversing light with negative refraction.  Physics today, 57(6):37–43, 2004.
  * [26] Adrien Rohfritsch, Jean-Marc Conoir, Régis Marchiano, and Tony Valier-Brasier.  Numerical simulation of two-dimensional multiple scattering of sound by a large number of circular cylinders.  The Journal of the Acoustical Society of America, 145(6):3320–3329, 2019.
  * [27] Ping Sheng and Bart van Tiggelen.  Introduction to wave scattering, localization and mesoscopic phenomena., 2007.
  * [28] Michael JA Smith and I David Abrahams.  Tailored acoustic metamaterials. part i. thin-and thick-walled helmholtz resonator arrays.  Proceedings of the Royal Society A, 478(2262):20220124, 2022.
  * [29] MJA Smith, PA Cotterill, D Nigro, WJ Parnell, and ID Abrahams.  Asymptotics of the meta-atom: plane wave scattering by a single helmholtz resonator.  Philosophical Transactions of the Royal Society A, 380(2237):20210383, 2022.
  * [30] Victor P. Tishkovets, Elena V. Petrova, and Michael I. Mishchenko.  Scattering of electromagnetic waves by ensembles of particles and discrete random media.  Journal of Quantitative Spectroscopy and Radiative Transfer, 112(13):2095–2127, September 2011.
  * [31] Daniel Torrent, Andreas Håkansson, Francisco Cervera, and José Sánchez-Dehesa.  Homogenization of two-dimensional clusters of rigid rods in air.  Physical review letters, 96(20):204302, 2006.
  * [32] L. Tsang, J. A. Kong, and T. Habashy.  Multiple scattering of acoustic waves by random distribution of discrete spherical scatterers with the quasicrystalline and Percus–Yevick approximation.  The Journal of the Acoustical Society of America, 71(3):552–558, March 1982.
  * [33] Victor Twersky.  On scattering of waves by random distributions. i. free-space scatterer formalism.  Journal of Mathematical Physics, 3(4):700–715, 1962.
  * [34] Victor Twersky.  On scattering of waves by random distributions. ii. two-space scatterer formalism.  Journal of Mathematical Physics, 3(4):724–734, 1962.
  * [35] V. K. Varadan, V. N. Bringi, V. V. Varadan, and A. Ishimaru.  Multiple scattering theory for waves in discrete random media and comparison with experiments.  Radio Science, 18(3):321–327, May 1983.
  * [36] 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 preprint arXiv:2106.13892, 2021.
  * [37] Peter Cary Waterman and Rohn Truell.  Multiple scattering of waves.  Journal of mathematical physics, 2(4):512–537, 1961.