Effective T-matrix of a cylinder filled with a random two-dimensional particulate

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 fibres. Here, we solve from first principles multiple scattering of scalar waves, including acoustic, for any frequency from a set of two-dimensional 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 scatterers, the expression of this effective T-matrix simplifies and reduces to the T-matrix of a homogeneous cylinder with an effective wavenumber k⋆. 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.


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 average1 .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

INTRODUCTION
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,

INTRODUCTION
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 V n and 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 V n , then the average scattered field is T n U n where the complex number 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 T n allows us to describe the scattering from any incident field, after decomposing the latter into the modes V n .Having simple expressions for 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 ⟨u sc ⟩, for a source with only one mode: u inc = V N (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 inc = V 0 (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 • 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 ⋆ .The method leads to an eigenvalue problem called the dispersion equation, whose eigenvalue provides k ⋆ .
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: where R := R − a.This result is remarkable because the above expression corresponds to the T-matrix of a homogeneous acoustic cylinder of radius R, sound speed c ⋆ = ω/k ⋆ and density ρ 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 ⋆ suffices to describe the random material.When the particles are not monopole scatterers, (1.1) is not exact.The exact formula is given by where the C n and 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 ⋆ .
For monopole scatterers, we have that (n ′ = 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 nonmonopole 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.

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 r i ∈ R 2 as shown in Figure 3 and is assumed to be confined in an area denoted by R i .
The propagation of waves in free space is governed by the 2D Helmholtz equation where ∆ := ∂ 2 x +∂ 2 y is the 2D Laplace operator and k ∈ R is the wavenumber of the homogeneous background.Consider the following two independent bases of the Helmholtz equation V n and

RANDOM PARTICULATE MATERIAL
Validation of 1.1 and 1.2 Figure 2: Compares various methods to calculate the component 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.
U n defined by U n (kr where k is the background wavenumber and (r, θ) are the polar coordinates of r, ie r = (r cos θ, r sin θ).The specific solutions V n (kr) have the particularity of being smooth while U n (kr) have a singularity at the origin and are outgoing solutions.Both the incident field u inc and the scattered field u sc are solutions of (2.1).While we assume that u inc is smooth and regular in the region that covers the particles2 the scattered field has to be a sum of outgoing fields from each particles centered at r i , as a result we can write The coefficients g n ∈ C are given parameters of the system and f i n ∈ 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 where T i n is the T-matrix of particle i, which can represent a wide range of particles [7].

RANDOM PARTICULATE MATERIAL
Set up of cylindrical particles inside a cylinder Illustration of a possible configuration of four particles.The particles are cylinders whose axes are aligned with the z-axis.The position r i ∈ R of the i 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 r i are therefore confined in the same disc R 1 shown in green.
To give an example, the expression of a T-matrix of a homogeneous particle with wavenumber k i , density ρ i and radius a i is given by The system (2.4) totally determines the scattering coefficients f i n , 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.

Particulate distribution
We describe each particle, say particle i, with two random variables: r i is the probability distribution of the centre of the particle in space, and λ i ∈ S describes all other properties of the particle (radius a i , density ρ i , etc.).
The centre of a particle r i is restricted in a region R i which depends on the particle's properties λ i .For example, r i has to be one radius away from the boundary ∂R which completely contains all particles.In this paper we assume that r i is distributed uniformly over the set R i (λ i ) defined by where

Joint particle distribution
Each particle is described by two random variables (λ i , r i ), i = 1, . . ., 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 g, which satisfies (cf.[16, eq.8.1.2]): where on the left is the joint law of two particles' position when their properties are known.
The pair correlation g describes how correlated any two particles are, when the positions and properties of all other particles are unknown.For example, if g = 1 for all values of its arguments then both r i and r j in the above are independent and uniformly distributed over R i and R j , respectively (in the limit J → ∞).

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 g(r i , λ i ; r j , λ j ) = g(|r j − r i |, λ i , λ j ).We assume the pair correlation is of the following form: where g(r, λ 1 , λ 2 ) = 0 when particles overlap, with a 12 ≥ a 1 + a 2 and in general a 12 can depend on λ 1 and λ 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, λ 1 , λ 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.

Definition of the effective T-matrix
The scattered field (2.3b) can be simplified after using Graf's addition theorem (A.1,ii) with x = kr and d = −kr i Radial pair-correlation function The pair-correlation 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 r 1 , where it is less (more) likely to find other particles.Finally, assuming g(r) = 1 for r ≥ b 12 means that a particle at r 1 becomes uncorrelated to particles which are further than b 12 .
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 Since the scattering problem is linear with respect to the incident field, ⟨F n ⟩ depends linearly on the modes g n of the incident field (2.3a): This relation defines the T-matrix 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 (2.15) where T n := T n,n .

THE EFFECTIVE WAVE METHOD
The function ⟨f n ⟩(r 1 , λ 1 ) needs to be determined before we can compute ⟨F n ⟩.In Appendix B.1 we show how to obtain the governing equation: 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 ⟨f n ⟩(r 1 , λ 1 ).This in turn allows us to compute ⟨F n ⟩ (2.16) and reach an explicit formula for T n,N by specifying g n = δ n,N in (2.14).

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

Overview of the effective waves method
The starting point is to assume that there exists k ⋆ ∈ C such that This assumption will greatly simplify the governing integral equation (2.17), from which we will be able to determine both k ⋆ and ⟨f n,N ⟩(r 1 , λ 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 ⟨f n,N ⟩(r 1 , λ 1 ) in the basis of functions V n (k ⋆ r 1 ) and substitute the decomposition into the ensemble wave equation.The coefficients F(λ 1 ) of the decomposition are then shown to be the eigenfunctions of an eigensystem which eigenvalue is k ⋆ .
3. The ensemble boundary condition.To determine the amplitudes of the eigenfunctions F(λ 1 ) we then use the ensemble boundary condition, which takes into account the incident field and the shape of the random material.

Separate microstructure and boundary
We follow the steps described above to determine the solutions ⟨f n ⟩(r 1 , λ 1 ) of (2.17).The first step uses (3.1), and some algebraic manipulations shown in Appendix C.1 to rewrite the

THE EFFECTIVE WAVE METHOD
governing equation (2.17) into two separate equations: 3) The terms J n ′ n (r 1 ), I n ′ n (r 1 ), and K n ′ n (r 1 ), respectively defined in (C.6) and (C.3), involve the function ⟨f n,N ⟩(r 1 , λ 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 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 ⋆ , which implies that k ⋆ only depends on the microstructure: the density of the particles n(λ 1 ), their properties provided by the T-matrix T n (λ 1 ), and the pair correlation g which explicitly appears in the quantity K n,n ′ (r 1 ) (C.3).We further discuss how to interpret k ⋆ 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.

The effective eigensystem
Since ⟨f n ⟩(r 1 , λ 1 ) satisfies (3.1), it can be decomposed into the modes where V n 1 is defined in (2.2).The unknowns k ⋆ and F nn 1 (λ 1 ) can be determined by substituting (3.4) into (3.2).The details are shown in Appendix C.2 with the resulting equation being where The above is a nonlinear eigenvalue problem, which is why we refer to k ⋆ as an eigenvalue.After calculating k ⋆ we can calculate the eigenfunctions F n ′ n 2 (λ 1 ) by solving the linear system (3.5),though in practice it is far better to calculate both k ⋆ and F n ′ n 2 (λ 1 ) using the modal decomposition as we do in Section 4.

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 αF nn 1 for any scalar α.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 R 2 , followed by choosing R 2 to be a cylinder with radius R − a 2 which results in the boundary condition:

Multiple effective wavenumbers
where N l is defined by (3.6).

An effective cylinder 4.1 Modal decomposition of the problem
In this section, we assume that the particles are confined in a circular disc R with radius R. Since (2.17) is linear with respect to g n , it can be decomposed in simpler and independent equations: let ⟨f n,N ⟩(r 1 , λ 1 ) be the solution when substituting g n = δ n,N in (2.17), then where we can recover the solution for any incident wave by using

AN EFFECTIVE CYLINDER
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: where α N ∈ 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 ⟨f n,N ⟩(r 1 , λ 1 ).To this end, we denote by R ϕ the rotation matrix of angle ϕ and replace r 1 with R ϕ r 1 in (4.1): where we changed the integration variable from r 2 to R ϕ r 2 , which is possible for any rotation as R 2 is a disc.Then from (2.2) we deduce the property which used in (4.4), and then multiplying both sides of the equation with e −i(N −n)ϕ leads to Now note that both ⟨f n,N ⟩(R ϕ r 1 , λ 1 )e −i(N −n)ϕ and ⟨f n,N ⟩(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 for any r 1 and ϕ.Let (r 1 , θ 1 ) be the polar coordinates of r 1 , then, without lose of generality, we then choose ϕ = −θ 1 which leads to Finally, because ⟨f n,N ⟩(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).

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 Then since V n form an orthogonal basis of functions, To simplify the effective eigensystem (3.5) we consider one mode at a time by taking g N 1 = δ N 1 −N , which used in (4.9) implies that we can substitute The result after some algebraic manipulations is 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 n(λ) = δ(λ − λ 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 S as a set of reals {t 1 , . . ., t S }.Then define a block vector F containing the entries F n (t s ) for n = −M, −M +1, . . ., M −1, M , for some finite M , and for s = 1, . . ., S, so that the eigensystem becomes where The parameter k ⋆ is then obtained by solving the equation det[I + M(k ⋆ )] = 0. (4.13)

Multiple wavenumbers
The dispersion equation (4.10) has infinitely many eigenvalues k p with p = 1, 2, . ... Consequently, ⟨f n ⟩(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 ⋆ 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 ⋆ in this paper.See Figure 5 for a typical distribution of the many eigenvalues of (4.10).

The modal ensemble boundary condition
To determine the α 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 = δ N 1 −N in (4.9) and then substitute the result into (3.8) to obtain which we can use to determine: 4 AN EFFECTIVE CYLINDER

The effective T-matrix
As discussed Section 4, the effective T-matrix T n,N can easily describe the average scattered wave for any incident wave through: where the ⟨F n ⟩, given by (2.16), are the average coefficients of the waves scattered from the whole cylinder R. Note the above holds for any choice of g N and T n,N does not depend on g N .
To calculate T n,N we substitute the modal decomposition (4.2) into (2.16) to obtain Comparing the above with the definition of the effective T-matrix (4.16), it is clear that 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 V n functions to conclude that where R 1 is the radius of the disc R 1 and we defined: Finally, we substitute α n given by (4.15), which results in 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 inc + u sc resulting from an incident plane wave and a point source.

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

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).material is made of sound-hard particles of radius one confined in a disc of radius 10.
which leads to ⟨f n,N ⟩(r 1 ) = 0 if n ̸ = 0 and, as a result of (4.3), F n,N = 0 if n ̸ = 0. Substituting this result into the effective T-matrix (4.21), and assuming that the radius R 1 of the region R 1 is the same for every type of particle λ 1 , we obtain:

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

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: This result is easily obtained from (2.15) with the specific choice g n = δ n,N , which substituted into (2.14)leads In other words, 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 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 ⟨u sc ⟩, which has been commonly used in the literature [16].For this

NUMERICAL RESULTS
second method we use (2.3b) to compute ⟨u sc ⟩(R, 0).Then from (5.1) we can also compute T N with T N = H −1 N (kR)⟨u sc ⟩(R, 0). ( 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: (5.4) While the terms n ̸ = 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 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.

Validation of the Effective Waves Method
We validate the Effective Waves method (4.18) against the Monte Carlo method (5.2) for several frequencies ω ranging in Ω := [0.05 : 0.015 : 1.5].To this end we define the relative error, averaged over frequencies: where T MC n (ω) is obtained following the Monte Carlo method (5.2) and T EWM n (ω) following the Effective Waves method (4.18).A few plots of T MC n (ω) and T EWM n (ω) are provided in Figures 2 and 8.The values of ϵ 0 , ϵ 1 , ϵ 2 , ϵ 3 , ϵ 4 for the cases of sound soft and sound hard particles are reported in Table 1.

CONCLUSION
Monte Carlo results: sound-soft particles Figure 8: Compares various methods to calculate the components T 1 , T 2 , T 3 and 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.

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 6 CONCLUSION ϵ 0 ϵ 1 ϵ 2 ϵ 3 ϵ 4 sound soft particles 6.75e −2 4.59e −2 5.40e −2 6.45e −2 7.45e −2 sound hard particles 2.66e −2 2.43e −2 2.33e −2 2.35e −2 2.37e −2 Table 1: Relative errors ϵ 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.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 ⋆ 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.

Appendices
A Bessel functions and translation matrices Given two points x, y ∈ R 2 , we have the following identities where The above formulas are direct consequences of Graf's theorem (see [21,] for instance).
Using these definitions, we define the ensemble average, and conditional ensemble averages, of a quantity A, which can depend on the positions and properties of all the particles, as where the domain of integration for Λ is over all possible particles positions and properties.The term ⟨A⟩(r i , λ i ; r j , λ j ) is the ensemble average of A conditional to (r 1 , λ 1 , r 2 , λ 2 ).Taking an ensemble average (B.1) on both sides of (2.12a) leads to where we used the definition of conditional probability: p(Λ) = p(r i , λ i )p(Λ (i) ), the definition of conditional average (B.2) to introduce the term ⟨f i n ′ ⟩(r i , λ i ), and then used that particles are indistinguishable3 .Finally, using (2.7) and (2.9) leads to formula (2.16).

C THE EFFECTIVE WAVES METHOD
over R 2 \ D(r 1 , a 12 ) in the form In this section, and in the whole paper, we only solve (4.1) for 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 r = r 2 − r 1 which leads to the integral The first integral over r 2 in (C.1) can be further simplified by using Green's theorem to replace the volume integral over R 2 \ D(r 1 , a 12 ) by surface integrals: given any two function smooth functions u, v which satisfy ∆u(r) + k ⋆ u(r) = 0 and ∆v(r) + kv(r) = 0, over a set Ω we have that With u(r 2 ) substituted for ⟨f n ′ ,N ⟩(r 2 , λ 2 ) and v(r 2 ) substituted for U n ′ −n (kr 1 − kr 2 ) we can use the above to deduce: where we defined (C.6) Finally, substituting (C.3) and (C.5) into the governing equation (2.17) gives

D BOUNDARY CONDITION FOR EFFECTIVE WAVES
The above now can be split into two separate equations by noting that the functions ⟨f n ⟩(r 1 ), J n ′ n (r 1 ) and K n ′ n (r 1 ) satisfy the wave equation with wavenumber k ⋆ , while V n (kr 1 ) and I n ′ n (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 ⋆ 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 ⋆ and write ⟨f n ⟩(r 1 , λ 1 ) in terms of eigenfunctions.Since ⟨f n ⟩(r 1 , λ 1 ) satisfies the wave equation (3.1), it can be decomposed into the modes where V n 1 is defined in (2.2).The unknowns k ⋆ and F nn 1 (λ 1 ) can be determined by substituting (C.8) into (3.2),which requires the term where the right side is a result of using Graf's addition theorem (A.1, i) in (3.4).By substituting (C.9) in K n ′ n (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 ′ ) are non-zero.
Likewise, we can perform the same simplification by substituting (C.9) in J n ′ n (r 1 ) (C.6).The simplifications result in  We can further simplify the above by using the orthogonality of the functions V n to obtain which holds for every N .When all particles are in a disk, then ∂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 R 2 , then n 2 = n 1 in (D.3) which reduces to where N n 1 is defined by (C.11) 2 .Substituting this into (D.5) leads to (3.8).

E Elementary proof that thee effective T-matrix is diagonal
We provide an elementary proof that T is diagonal when the particles are confined in a disk of radius R. To this end, we consider the scattering from the modal source u N inc obtained for g n = δ n,N : The notation F n,N is the corresponding F n to the specific incident field with g n = δ n,N (compare with (2.12a)).We then denote by f i n,N (σ) the resulting solution of (2.4) for the specific configuration σ = r 1 , . . ., r J .The rotation by angle ϕ of the particles r 1 , . . ., r J correspond to another valid configuration (because the random material is cylindrical), for which the solutions are given by f i n,N (R ϕ σ) = e i(N −n)ϕ f i n,N (σ) (E.1) This tells us how the rotation of a configuration modifies the coefficient F n,N , using (2.12a): F n,N (R ϕ σ) = This analysis proves that only the diagonal terms of the effective T-matrix are nonzero and can be estimated by T n,n = ⟨F n,n ⟩.

Figure 4 :
Figure 4: Typical plot of the pair correlation function.

Figure 5 :
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 ω = 1.k ⋆ 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.

. 23 )
This corresponds to the T-matrix of a homogeneous cylinder of radius R = R − a, sound speed c ⋆ = ω/k ⋆ and mass density ρ ⋆ = ρ where ρ is the density of the host medium (cf.(2.5)).