Effective Waves for Random Three-dimensional Particulate Materials

How do you take a reliable measurement of a material whose microstructure is random? When using wave scattering, the answer is often to take an ensemble average (average over time or space). By ensemble averaging we can calculate the average scattered wave and the effective wavenumber. To date, the literature has focused on calculating the effective wavenumber for a plate filled with particles. One clear unanswered question was how to extend this approach to a material of any geometry and for any source. For example, does the effective wavenumber depend on only the microstructure, or also on the material geometry? In this work, we demonstrate that the effective wavenumbers depend on only microstructure and not the geometry, though beyond the long wavelength limit there are multiple effective wavenumbers. We show how to calculate the average wave scattered from a random particulate material of any shape, and for broad frequency ranges. As an example, we show how to calculate the average wave scattered from a sphere filled with particles.

To develop new sensors to monitor in real time, large quantities of materials, we first need to understand how waves scatter from these dense particulates, and develop efficient models to describe this scattering.
Every particle counts. Any method that uses waves to probe a particulate material needs to consider how each particle scatters waves. This is because both particle properties and positions influence the total scattered waves, as shown in Figure 1.
Although it is possible to numerically simulate scattered waves from a specific arrangement of particles, these numerical methods are computationally too intensive for most practical applications. For example, one droplet of most emulsions will contain hundreds of millions of oil particles, whose positions are unknown. The most successful methods avoid these heavy computations by replacing the material with an equivalent homogeneous material [7,44]. This equivalent homogeneous material is calculated by taking an ensemble average.
The ensemble average. Ensemble averaging not alone simplifies the calculations, it is also the route to devising measurements which do not depend on the positions of the particles, which are unknown. One way to do this is to take the average of the scattered field. This average can be taken over space or over time (for ergodic systems). Both of these types of average measurements eliminate the need to know the particle position, and so lead to reliable measurements [18,44].
If u(x) represents the transmitted wave field, measured at some distance x, then, for a plane wave source propagating along the x-axis, it is common to approximate the ensemble average as a plane wave of the form u(x) ≈ Ae iω(x/c * −t)−αx , (1.1) where c * is the (effective) wave speed, α the rate of attenuation, and A the average transmission coefficient. The process of calculating the ensemble average links the measurables c * , α, and A to the particles; it is this link which drives many sensing methods. It is common to combine c * and α into one quantity, the complex effective wavenumber: k * = ω/c * + iα.
What is known. One scenarios has been mostly clearly understood: a plane wave incident on a halfspace or plate region filled with particles. This setup has, what we call in this paper, planar symmetry. For planar symmetry, in the limits of low frequency or low volume fraction, there are explicit formulas [47,35,41,34,5,4], and an understanding on how to calculate wave reflection and transmission [42,26]. Further, the effective wavenumbers for planar symmetry have also been rigorously deduced [24] (given typical statistical assumptions), though there is often more than one effective wavenumber for the same fixed frequency [26,56,57]. One clear question that remained was how to extend this approach to a material with any geometry and for any source? For example, is the effective wavenumber k * the same for other geometries? There has even been evidence [28] that the effective properties (and wavenumber) depend on the geometry of the material. If this were true, these effective wavenumbers would not be very useful, as they would change for every sample of the same material.
In the electromagnetic community, the analysis of effective wave properties in particulate media has a long tradition. Some of the most significant contributions are collected in textbooks, e.g., [51,53,52] and journal literature [44,50]. With a few exceptions, the analysis deals again with planar symmetry.
This paper. Here we develop the theory for effective waves and wavenumbers for materials in any geometry. The key to achieve this is to use the representation: where φ p (x) is a function that satisfies ∇ 2 φ p (x) + k 2 p φ p (x) = 0. This representation allows us to deduce a dispersion equation for the k p that does not depend on the material geometry. This question of whether the geometry changes the effective wavenumbers has been raised in previous studies [28].
In this paper we present a framework for effective scalar waves in any material geometry, and then specialise to a material shaped as a sphere and a plate. This allows us to design highly efficient numerical methods for these cases.

A collection of particles
We begin with the deterministic many-particle scattering problem and use the Null-field approach [32]. Consider N different particles, where the i-th particle is centred at the location r i as shown in Figure 2. 2 The radius of the minimum circumscribed sphere, centred at r i , is a i , i = 1, 2, . . . , N . We assume that no minimum circumscribed spheres intersect. Each particle can have a different shape and material properties.
The particles are located in a homogeneous, isotropic media with wavenumber k, which is either a real number or a complex number with a positive imaginary part.
The prescribed sources are located in the region V in , which is a region disjoint to all particles, 3 and these sources generate the field u in (r) everywhere outside V in .
For a point r, outside of the circumscribed spheres of all particles, we can write the total field u(r) as a sum of the incident wave u in (r) and all scattered waves in the form [30,32,36] u(r) = u in (r) + u sc (r), u sc (r) = N i=1 n f i n u n (kr − kr i ), (2.1) where we assumed |r − r i | > a i for i = 1, 2, . . . N , the f i n are coefficients we need to determine, and for convenience we use scalar spherical waves: u n (kr) = h (1) (kr)Y n (r), (outgoing spherical waves) v n (kr) = j (kr)Y n (r), (regular spherical waves) (2.2) where r = |r|, and n denotes a multi index n = { , m}, with summation being over = 0, 1, 2, 3 . . . and m = − , − + 1, . . . , −1, 0, 1, . . . , . For more details, see Appendix A. The spherical Hankel and Bessel functions are denoted h (1) (z) and j (z), respectively. The field n f i n u n (kr − kr i ) is the wave scattered from particle-i.

Incident field
We assume the incident field is generated outside of all particles, see Figure 2, so it has an expansion in regular spherical waves u in (r) = n g n v n (kr) = nn g n V nn (kr i )v n (kr − kr i ), (2.3) where for the last equality we used a translation matrix of the regular spherical waves, V nn (kr i ), to write the incident wave in terms of spherical waves centred at r i . See Appendix B for details.

Scattered field
The coefficients f i n from (2.1) are determined by using the T-matrix to relate the field incident on the i-th particle, u(r)− n f i n u n (kr−kr i ), to the wave scattered from the i-th particle, n f i n u n (kr−kr i ), which leads to [30,32,36] f i n = n T i n V n n (kr i )g n + N j=1 j =i n T i n U n n (kr i − kr j )f j n , i = 1, 2, . . . , N, where U nn is the translation matrix of the outgoing spherical waves u n , see Appendix B. Above we have used a diagonal T-matrix which assumes a spherical particle; later we will explain how this leads to the solution for non-spherical particles whose orientation is independent of position and properties.
Equation (2.7) is very difficult to calculate when the number of particles N is large. Nevertheless, there are several software packages making substantial progress, e.g., MSTM (Multiple Sphere T Matrix) [39,21,22].
The T i n depend only on the properties of the i-th particle, while the scattering coefficients f i n depends on the positions and properties of all the particles. For example, for acoustics, and a homogeneous spherical particles, we would have [35]: where γ i = ρ i k/(ρk i ), a i is the particle radius, ρ is the background density, while ρ i and k i are the density and wavenumber of the particle.

Ensemble averaging
Even if the position and properties of all particles were known, it is still very challenging to solve (2.7) for a large number of particles, say, over 10 6 . Also, many sensors can not even measure f i n , but instead measure the scattered field averaged either in time or space [18,44]. For these reasons it makes sense to calculate the ensemble average scattered waves. The first step towards achieving this is to introduce a probability for the particles having certain properties and positions [36,30,51,50].

Statistical assumptions
To describe the properties and shape of the i-th particle, we will use the variable λ i , which allows us to define T n (λ i ) := T i n for every i. This means that the f i n , governed by (2.7), depend on the positions r 1 , r 2 , . . . , r N and the properties λ 1 , λ 2 , . . . , λ N of all the particles.
To ensemble average we need to assign a probability density for any configuration r 1 , r 2 , . . . , r N , and any properties λ 1 , λ 2 , . . . , λ N . The first step is consider the r i and λ i as random variables. Next we assume that the particle properties λ i are sampled from the same domain S. For example, if λ i = a i , the radius of particle-i, for every i, then we could choose S = [A 1 , A 2 ] so that all λ i ∈ S, i.e., we restrict all particle radii in some interval. For particle origins r i we can not restrict them all to the same domain because the particles may have a different sizes. So instead we choose a different domain for each, that is, for a given λ i we have that r i ∈ R i . For example, if all the particles were contained in a sphere with of radius R, then a particle with radius a i would have its origin r i restricted in a sphere of radius R − a i . That is, R i would be a sphere of radius R − a i . For more details on ensemble averaging for multi-species particles see [27].
The main parameters we use to describe the average particulate material are a ij (the minimal allowed distance |r i − r j | between particle i and particle j), (3.2) where |R i | is the volume of R i and p(λ i ) is the probability density of the particle having the property λ i . In this paper we allow the minimal distance between two particles a ij to be larger or equal to the sum of the particle radii a i + a j . Note we committed an abuse of notation for the function p, and will continue to do so.
Let p(r i , λ i ) be the probability density of having a particle centred at r i ∈ R i with λ i ∈ S, after ensemble averaging over all other particle positions and properties. If we assume that r i is equally likely to be anywhere in R i we obtain We also need to define conditional probabilities: where M is any integer smaller than the number of particles N .
To solve the ensemble average equations, the probability function for two particles p(r i , λ i ; r j , λ j ) needs to be given. To achieve this, we use an assumption called hole correction, which assumes that any two particles are equally likely to be anywhere within regions 4 , except that their minimum circumscribed spheres do not overlap [16,17]: To deduce the above for |r i − r j | ≥ a ij we used where p(r i ; r j |λ i ; λ j ) is the probability density of having one particle centred at r i , knowing that it has the property λ i , and another particle at r j , knowing that it has the property λ j . The approximation above assumes that the volume of one particle is negligible in comparison to the volume of its confining region.
To help interpret hole-correction (3.5) we will do some extra calculations. For simplicity, we assume that the particle properties λ i and λ j are independent of each other to reach where for the last approximation we used (3.5) and (3.1). An alternative way to calculate the above is to approximate p(r j , λ j |r i , λ i ) for its expected value in r i and λ i when |r i − r j | ≥ a ij , that is which when using (3.3) leads to the same conclusion as hole correction (3.7). Here, dr i is the volume measure of the region R i . Later, we show that the quasi-crystalline approximation (3.12) makes an approximation which is analogous to (3.8).
We can now define the ensemble average of f 1 n as where the above integrals are over all feasible values for the particles positions r i and properties λ i .
We also need the conditional ensemble averages, which we define as Note that in (3.10) we are holding the first particle's position r 1 and properties λ 1 fixed while averaging over the other particles. In (3.11) we are averaging f 2 n while holding the first and second particles positions r 1 , r 2 and properties λ 1 , λ 2 fixed.
For consistency and simplicity, we will use an approximation for f 2 n (r 1 , λ 1 ; r 2 , λ 2 ) which is analogous to both (3.8) and (3.5), and is called the quasi-crystalline approximation: That is, we replace f 2 n (r 1 , λ 1 ; r 2 , λ 2 ) for its expected value in r 1 and λ 1 , see [27] for a brief discussion on the topic. This is a standard approach used across statistical physics. It is called a closure approximation [33,1].
Because the particles only differ due to their position r i and properties λ i , we have that f i n (r i , λ i ) = f j n (r j , λ j ) for any i and j (all particles with the same properties are indistinguishable). This is why we now define: f n (r j , λ j ) := f j n (r j , λ j ) for j = 1, 2, . . . , N. (3.13)

Average scattered field
To calculate the ensemble average scattered field we first choose a point r outside of the material, where we want to measure the scattered field. For example, turning to Figure 3, the point r needs to be outside of R 2 and at least one particle radius a 2 away from the boundary of R 2 . Then we multiple both sides of (2.1) by p(r 1 , . . . , r N , λ 1 , . . . , λ N ) and integrate over all possible particle positions and properties to reach where u in (r) = u in (r), because the incident wave does not depend on the particle configuration, and where we used (3.3), (3.4), and (3.13). Note that to take the limit N → ∞, it normally makes sense to fix the number density n(λ i ) and the probability p(λ i ), and then allow the volume of the region |R i | to grow with N .
We can rewrite the above when |r| > |r 1 | for every r 1 ∈ R 1 . In this case, we can use the translation matrix (B.1) for u n to obtain u sc (r) = n F n u n (kr), with F n = n S n(λ 1 )

R1
V n n (−kr 1 ) f n (r 1 , λ 1 ) dr 1 dλ 1 . (3.16) The F n are then the average scattering coefficients of the whole material.

Average governing equations
To calculate f n (r 1 , λ 1 ), we need to ensemble average the governing equation (2.7). To achieve this, we set i = 1, multiple both sides of (2.7) by p(r 2 , λ 2 ; r 3 , λ 3 ; . . . ; r N , λ N |r 1 , λ 1 ), and then integrate over all feasible positions and properties while holding r 1 and λ 1 fixed. Then to transform the result into an equation where f n (r 1 , λ 1 ) is the only unknown we use (3.7) and (3.12), to obtain for all r 1 ∈ R 1 and λ 1 ∈ S, where we define used R 2 \ B(r 1 ; a 12 ) = {r ∈ R 2 : r ∈ B(r 1 ; a 12 )} andn(λ 2 ) = N −1 N n(λ 2 ). The system (3.17) can be used to solve for f n (r 1 , λ 1 ) for any given material geometry R 1 and any T-matrix T n . If all particles were the same, i.e., same shape and properties, then (3.17) would be equivalent to [36,Equation (4.13)] and [30,Equation (12)]. If we considered a two dimensional material, with different types of particles, then (3.17) would be equivalent to [27,Equation (3.6)].

Averaging non-spherical particles
As a side, we explain how (3.17) can accommodate particles which are not exactly spherical, as shown in Figure 2.
Assume we have non-spherical particles with a T-matrix T nn (λ j , τ j ), which depends on the particle properties λ j and orientation τ j . If every particle's orientation is statistically independent from everything else a , then we can set the T n in (3.17) to equal where p(τ j ) is the probability density of particle j being rotated by a τ j angle. Note that in general τ j could represent three Euler angles. In particular, if the particle is equally likely to be oriented in any direction then its T-matrix T nn (λ j , τ j ) averaged over every angle τ j becomes diagonal [45,54].
a Including the minimal allowed distance between any two particles (3.2).

Symmetry reductions
Before solving (3.17) to determine the field f n (r 1 , λ 1 ), we first look at how to use symmetries to represent f n (r 1 , λ 1 ) in a reduced form.
We could apply symmetry reduction directly to the governing integral equation (3.17). It is, however, simpler to just impose symmetries on the average scattered wave u sc (r) (3.15) and then deduce the resulting symmetry for f n (r 1 , λ 1 ) as we demonstrate below. To omit a heavy notation, we will in this section omit the dependence of f n on λ 1 and the integrals over the species S.
Azimuthal symmetry: we expect this symmetry when the total scattered wave u sc does not change when rotating the measurement point r around the z−axis. This occurs, for example, for the incident plane wave u in (r) = e ikz and a spherical material region R 1 = {|r| ≤ R : r ∈ R 3 } centred at the origin. When azimuthal symmetry is present, we expect u sc (Pr) = u sc (r) for every |r| ∈ R 1 , where we define the operator P such that Pr is a φ 0 rotation of the vector r around the z-axis. Note that when R 1 is a sphere, then the above should hold true for r ≥ R.
To determine the consequences of this symmetry, we turn to the average scattered wave (3.15) and rewrite in the form where we used a change of variables r 1 → Pr 1 , and used (2.2) and (A.1) to substitute u n (kPr − kPr 1 ) = e imφ0 u n (kr − kr 1 ). Notice that the volume measure dPr 1 = dr 1 . Equating the above to u sc (r) and using (3.15) then suggests that f n (r 1 ) = f n (Pr 1 )e imφ0 . Then by using a spherical coordinate system (r 1 , θ 1 , φ 1 ) for r 1 , and by choosing φ 0 = −φ 1 (without loss of generality) we find that (arguments in spherical coordinates) for every φ 1 . This symmetry can now be verified by checking that the right hand-side is a solution to (3.17), though this is a longer calculation.
Planar symmetry: For an incident plane wave u in (r) = e ik·r and the material region R 1 = {r ∈ R 3 : z > 0}, we expect the average scattered wave to satisfy the planar symmetry: where r 0 = x 0x + y 0ŷ . If we then use (3.15) in the above we find that where for the second equation we changed to the variable of integration x 1 + x 0 → x 1 and y 1 + y 0 → y 1 . For the above to be equal to for every x, x 0 , y, y 0 , and z < 0 suggests that f n (r 1 ) = f n (r 1 − r 0 )e ikk·r0 , then by choosing r 0 = x 1x + y 1ŷ , we find that f n (r 1 ) = f n (z 1ẑ )e ikk·(x1x+y1ŷ) .
(3.21) This symmetry can be verified by checking that the right hand-side is a solution to (3.17).
One case that combines both azimuthal (3.19) and planar symmetry (3.21) is the incident plane-wave e ikz (k =ẑ) and material region z > 0. In this case (arguments in spherical coordinates), where the first equation is due to planar symmetry (3.21) and the second is due to azimuthal symmetry (3.19). Equation (3.22) can only be true for every φ 0 when This result will be used later to reach a simplified dispersion equation.

Effective wavenumbers 4.1 Wave decomposition
Much of the literature has focused on solving (3.17) by assuming that the unknown field f n (r 1 , λ 1 ) satisfies a wave equation for the spatial variable r 1 and for some effective wavenumber k * . This assumption implies that the average transmitted fields u(r) also satisfy a wave equation [42]. Recent results [26,24] have demonstrated that for a half-space, the exact solution for f n (r 1 , λ 1 ) is a sum of plane waves, each with a different wavenumber. Here we generalise this result by considering f n (r 1 , λ 1 ) to be a sum of isotropic waves of any type, i.e., not necessarily a plane wave, and the particulate material to occupy any region.
In general, we propose the representation where the Laplacian ∇ 2 rj is taken in terms of r j . Our first major result is to calculate the effective wavenumbers k 1 , k 2 , · · · , and to demonstrate that they depend only on the particle properties, and not on the geometry of the region enclosing the particles. Although the geometry of material and the incident wave will determine which of these wavenumbers are excited. A another major result, is that most of the effective wavenumbers are highly attenuating, which implies that the series (4.1) rapidly converges to the exact solution. In the remained of this section, we show how to deduce a system that determines the f p,n (r 1 , λ 1 ) that is decoupled from the material geometry.
To simplify the governing system (3.17) we note that by definition (B.1) the translation matrix U n n (kr 1 − kr 2 ) satisfies a wave equation in either r 1 or r 2 with wavenumber k. This and the representation (4.1) leads to , for r 1 ∈ R 1 , (4.2) Then for r 1 ∈ R 1 (a 12 ) = {r ∈ R 1 : d(r, ∂R 1 ) ≥ a 12 }, we can integrate both sides over r 2 ∈ R 2 \ B(r 1 ; a 12 ) and apply Green's second identity to reach, where dA 2 and dA are the surface elements for r 2 and r, respectively, and ν 2 and ν are outward pointing normal vectors to the surfaces ∂R 2 and ∂B(0; a 12 ), respectively (see Figure 3). To reach (4.5) we changed the integration variable to r = r 2 − r 1 .
Both integrals I p (r 1 ) and J p (r 1 ) depend on the indices n, n and on the state variable λ 1 and λ 2 , which we omit to avoid a heavy notation.

The boundary layer
The region R 1 \ R 1 (a 12 ) is often called a boundary layer. The simplification (4.3) only occurs when r 1 is not in this boundary layer, that is when r 1 ∈ R 2 (a 12 ). This is because then ∂R 2 ∩ ∂B(r; a 12 ) = ∅, and therefore the boundary of R 2 \ B(r 1 ; a 12 ) becomes ∂R 2 ∪ ∂B(r 1 ; a 12 ).
In this paper we will not discuss how to evaluate the system (3.17) for r 1 ∈ R 1 \ R 1 (a 12 ). Evaluating in this boundary layer is used to combine the different fields f p,n , but this is only needed for very strong multiplescattering [24,26]. We also note that for a slab geometry, it is possible to evaluate the complex integrals when r 1 is located in the boundary layer, see [30,31].
By substituting equations (4.1)-(4.5) into the governing equation (3.17), and assuming r 1 ∈ R 2 (a 12 ), we obtain The key to simplifying (4.6) is to note that both f p,n (r 1 ) and J p (r 1 ) satisfy a wave equation with wavenumber k p and spatial position r 1 , whereas V n n (kr 1 ) and I p (r 1 ) satisfy a wave equation with wavenumber k. This enables us to use Theorem C.1 to conclude that both valid for r 1 ∈ R 1 (a 12 ). Now it is clear that (4.7) is independent of both the region of particles R j and the incident field. We will show how the effective wavenumbers k p can be completely determined from (4.7). On the other hand, equation (4.8) depends on both the region and incident wave, and will lead to a restriction on how to combine the f p,n (r, λ). Equation (4.8) is sometimes called the extinction equation.
Both (4.7) and (4.8) can be further simplified by expanding the fields in terms of orthonormal functions, which we do in the section below.

The dispersion equation
The effective wavenumbers k p and much about the fields f p,n (r 1 , λ 1 ), can be calculated just from the ensemble wave equation (4.7). Depending on the symmetries of f p,n (r 1 , λ 1 ) we can reach different dispersion equations, the most general of which just assumes that f p,n (r 1 , λ 1 ) is a smooth field.

Effective regular waves
Here, we determine the effective wavenumbers k p from the ensemble wave equation (4.7). To do so, we use an origin O, for our coordinate system, located in R 2 (a 12 ), and we expand f p,n (r 1 , λ 1 ) in a series of regular spherical functions of the wave equation in the domain R 2 (a 12 ).
This regular series takes the form: where the coefficients F p,nn1 (λ 1 ) are to be determined, and we used the translation matrix V n1n2 (k p r) in the second expansion. The convergence of this series depends on the behaviour of the expansion coefficients F p,nn1 (λ 1 ), which in turn depend on both the confining geometry and the incident field. At this stage, we assume the series is convergent. In Appendix D we show how substituting the above into the dispersion equation (4.7) leads to where the c n nn1 are numbers which are defined in Appendix B, and (5.5) By using the orthonormal property of spherical harmonics (5.3) reduces to (the regular eigen-system) (5.6) As the above is a linear system of equations for the unknowns F p,n n1 (λ 1 ), we can rewrite it in the form 5 where the second equations holds for a non-zero F. This determinant equation can be used to find all effective wavenumbers k p for any geometry. Although (5.7) contains all possible effective wavenumbers it is computational simpler to solve the planar dispersion equation which also contains all viable effective wavenumbers, as we will show below. Solving (5.7) can be numerically difficult for two reasons: 1) the roots of (5.7) are multiple roots with different multiplicities, and 2) there are many spurious roots, as discussed in the optional box below.
The plane wave dispersion, and other reduced dispersion equations, are calculated by restricting the form of f p,n (r 1 , λ 1 ) through the use of symmetry reductions as shown in Section 3.4.

Spurious wavenumbers k p
One problem with using (5.7) to find the wavenumbers k p is that (5.7) has spurious roots. That is, it has solutions k p which are not solutions to the ensemble wave equation (4.7). Figure ? shows some of these spurious solutions.
These spurious solutions appear when truncating the index n 2 in F p,nn2 and then solving (5.7). Calculating the eigenvectors F p,nn2 of these spurious solutions k p would then lead the series in (5.1) to rapidly diverge, which is physically not viable. That is, when truncating the system (5.6) for 2 ≤ L 2 , it would also be natural to truncate 1 ≤ L 2 . However, neglecting the terms F p,n n1 with 1 > L 2 on the right hand-side of (5.6) is only approximately correct if |F p,n n1 | is small or at least getting smaller when increasing 1 . These spurious roots k p on the other hand lead to |F p,n n1 | which increase with 1 .

Effective azimuthal waves
When both the incident wave and material region share a rotational symmetry around the z−axis we can reach a reduced representation of f p,n (r 1 , λ 1 ) by using (3.19). For example, this occurs when the incident wave is u in (r) = e ikz and the material region is a sphere centred at the origin.
Combining the symmetry (3.19) with the representation (5.1) leads to the form . In other words, in terms of the representation (5.1) we have that F p,nn1 = 0 unless m 1 = −m and 1 ≥ |m|.

Effective plane-waves
Here, we restrict f p,n (r 1 , λ 1 ) by imposing planar symmetry (3.21). This will allow us to deduce simpler dispersion equations, as well as deduce reflection and transmission from a plate.
By combining the symmetry (3.21) with the wave equation (4.1) we first conclude that where k p · r 1 is simply a sum of element wise multiplication without conjugation, even though k p is a complex vector.
In the appendix E we show that substituting the above into (4.7) leads to and that, when considering only one type of particle, the above reduces to an equation which is found in much of the literature [35,16,51,38,13,27]. 12) and the angles θ p and φ p can be complex numbers, meaning that we may have |k p | = 1, but we do have thatk p ·k p = 1 for the real inner product.
Equation (5.11) can be turned into a determinant equation, much like (5.7), from which we can calculate effective wavenumbers k p : The form of the eigensystem (5.11) seems to suggest that the k p depend onk p . However, a direct (though cumbersome) evaluation of the resulting (truncated) determinant system would confirm thatk p has no contribution. Further, the more general eigensystem (5.6) does not depend onk p . As a sanity check, we can explicitly show that every solution to (5.11) is also a solution to (5.6). To achieve this we rewrite the solution (5.10) in the form (5.1) by using a plane-wave expansion (A.2) in (5.10) to obtain: Because the above is in the form (5.1) and the field f p,n (r 1 , λ 1 ) satisfies the general dispersion equation (4.7) (when (5.11) is satisfied) then F p,nn1 and k p must also satisfy the regular eigensystem (5.6).
As the k p are independent ofk p , we can choose anyk p to calculate the k p . We exemplify for a single species: takê k p =ẑ so that Y n1 (k p ) = √ 2 1 + 1/ √ 4πδ m10 , then the k p must satisfy for a single species 6 , wheren = n(N − 1)/N and n is the number density of particles. This equation can be even further simplified when there is azimuthal symmetry, as is the case for the incident plane-wave e ikz and material region z > 0. In this case we can apply the symmetry (3.22) to (5.15) and reach All the k p that satisfy (5.16) also satisfy (5.15), however, there are solutions to (5.15) which do not satisfy (5.16), see Figure 4. That is, it is not possible to excite all effective wavenumbers when considering only direct incidence u in (r) = e ikz . In other words, one type of experiment (one type of incident wave and material geometry) can only excite a portion of all the effective wavenumbers 7 .

Plane-wave dispersion has all viable effective wavenumbers
Because the representation (5.1) is more general than a plane-wave representation (5.10), we know that all solutions k p to the plane-wave dispersion (5.13) must also satisfy the more general regular dispersion (5.6). There is even an explicit conversion from plane-wave solutions to the regular solutions (5.14). However, it is not at all obvious that all viable solutions k p of the regular dispersion (5.6) must satisfy plane-wave dispersion (5.13).
To show that all viable wavenumbers k p that satisfy the regular dispersion (5.7) must also satisfy the plane-wave dispersion (5.13), we need to rewrite the expansion (5.1) in terms of plane-waves. We achieve this by using where Ω q is the solid angle of the radial unit vectorq. The above can be verified by using a plane-wave expansion (A.2) for e ikpq·r1 , to write Written in this form, f p,n (r 1 , λ 1 ) is now a superposition of plane waves all with the same wavenumber k p but with different directionsq. Note that the F p,nn1 depend on k p but are independent of theq. To find a dispersion equation we can repeat the same steps that led to (5.11) to reach: As the mapq → F p,n (q, λ 1 ) is smooth, we can conclude that the integrand in the above is also a smooth function ofq, in which case, the above can only be zero for every r 1 when the integrand is zero for everyq. That is, the k p and F p,n (q, λ 1 ) have to satisfy the plane-wave dispersion equation (5.11), with the same k p for everyq, which in turn implies that k p has to satisfy the determinant equation (5.13). 8

Effective properties in the long wavelength limit
By taking the limit where the incident wavelength is long compared to the particle diameter, we can calculate the effective properties directly from any of the dispersion equations. This procedure is explained in detail in [47,43,27].
It is particularly interesting to calculate the effective properties from the regular dispersion equation (5.7), because this equation holds for any material geometry, which then gives us confidence that the effective properties are truly properties of the material's microstructure and medium, and not its geometry.
Here we calculate the effective properties for spherical particles (2.8) and acoustics. To achieve this, we need to consider the limit k → 0, starting with the T-matrix coefficients which scale with k in the form and with β j = ρ j ω 2 /k 2 j and β = ρω 2 /k 2 being the bulk modulus of the j-th particle and of the background medium, respectively. For particles with any shape, there is a similar result for their T-matrix when assuming that they scatter only monopole and dipole waves [55]. 8 In more detail, let Ω k f (k)e ikk·r dΩ k = 0, ∀r ∈ R, and expand f (k) in spherical harmonics, i.e., f (k) = n fnYn(k). We assume this relation holds in a ball of radius R, centred at the origin. Then, using the transformation (5.17), we obtain n i fnvn(kr) = 0, ∀r ∈ R. Orthogonality of the spherical harmonics implies fnj (kr) = 0, r ∈ [0, R], ∀n, and fn = 0, since the zeros of the spherical Bessel functions are isolated points on the real axis.
To facilitate the next steps, we rewrite the effective wavenumber k * and expand where c * is the constant effective phase speed.
By substituting the above into (5.7), and then expanding up to leading order in small k, we find three possible solutions for the effective phase speed c * . Two are these solutions are non-physical, because they do not satisfy the plane-wave dispersion, as discussed in Section 5.4 and at the end of Section 5.1. The only remaining physically viable solution is where, just for this section, we define and likewise for ∆ρ j , where ϕ(λ j ) dλ j is the volume fraction of particles with the properties λ j , so that S ϕ(λ j ) dλ j = ϕ, the total particle volume fraction. Note that for spheres ϕ(λ j ) = 4πa 3 j 3 n(λ j ). By writing c 2 * = β * /ρ * we can now identify the effective bulk modulus β * and density ρ * as which is in fact the multi-species version of a classical formula [2, equation (9)] (in the absence of viscosity) and many others [43]. When performing this same procedure for the plane-wave dispersion equation (5.11) or for azimuthal symmetry (5.9) we recover the same effective properties. Also note, that for a single species, with the properties β 1 and ρ 1 , we recover the correct limits: when ϕ → 0 we get β * → β and ρ * → ρ, and when ϕ → 1 we get β * → β 1 and ρ * → ρ 1 .

Numerical effective wavenumbers
Here we numerically explore the effective wavenumbers which solve the dispersion equations: (5.7) with azimuthal symmetry (5.8), planar symmetry (5.15), and the combined planar with azimuthal symmetry (5.16). The material properties we use are shown in the Table 1. ϕ = 30% (particle volume fraction) ka o = π/8 (non-dimensional particle radius) ρ o /ρ = c o /c = 0.1 (void particle properties) ρ o /ρ = c o /c = 10.0 (solid particle properties) Table 1: For numerical experiments we use identical spherical particles with the T-matrix (2.8) and the material properties given above, where c = ω/k and c o = ω/k o are the background and particle wave speed, respectively. Note that we give the properties of the particles relative to the background properties. The properties of the void particle are an example of strong scatterers, while that of the solid particle are an example of weak scatterers.

Algorithm -effective wavenumbers
A minimal algorithm to calculate the effective wavenumbers k p for any material region. 1. Choose the particle statistics by choosing: 1.1. a wavenumber k, where ω = ck and c is the background wave speed. 1.2. the coefficients of the T-matrix T n . One example is given by (2.8).
1.3. the function p(λ 1 ) detailed in Section 3.1. For just one species p(λ 1 ) = δ(λ 1 − λ * ). 1.4. the exclusion distance a 12 , with a 12 = 1.001(a 1 + a 2 ) being a common choice. 2. Choose a truncation ≤ L, for the in n = ( , m) and in F p,nn1 or F p,n based on how the T n decay. 3. Calculate the effective wavenumbers k p by solving (5.16) or (5.15).
In Figures 4-5 we show the result of using Algorithm 4 above to calculate the different effective wavenumbers when using the properties in Table 1. In both cases, the wavelength λ is sixteen times larger than the particle radius. The important messages to take-away from these figure are: 1. The more general regular dispersion equation has spurious roots, which are the ones that do not satisfy the planar dispersion, as discussed in Section 5.4. 2. There can be two effective wavenumbers with lower imaginary part (and are not spurious roots), as shown in Figure 4. These two will dominate calculate the ensemble average transmission and scattering, as the other wavenumbers will be very difficult to excite. For weaker scatterers there tends to be only one wavenumber with lower imaginary part as shown in Fgiure 5. 3. The simpler combined planar-azimuthal dispersion (5.16) equation contains the two most important wavenumbers. This seems to hold in general.
We remark that in most cases we find that there is only one wavenumber with a low imaginary part, and knowing this wavenumber is often enough to accurately calculate the ensemble average transmission and scattering. This is exactly what we do in the next sections.   Table 1 for void particles and incident wavenumber times particle radius ka o = π/8, which means the wavelength is 16 times longer than the particle radius.

Scattering from a sphere filled with particles
If one effective wavenumber k 1 has a significantly smaller imaginary part than the other wavenumbers, e.g. Im k 1 Im k p for p = 2, 3, . . ., then f n (r 1 , λ 1 ) ≈ f 1,n (r 1 , λ 1 ), where f 1,n (r 1 , λ 1 ) is the wavemode associated with k 1 as shown  Note that the wavenumber closest to the x-axis will be the most important, as it has a smaller imaginary part. We used the material properties in Table 1 for solid particles and incident wavenumber times particle radius ka o = π/8.
in the representation (4.1). This occurs for a number of scenarios including: weak scattering, low frequency, or low volume fraction. In this case, we can explicitly calculate the average scattered and transmitted waves for many different material geometries. To achieve this, for each material geometry, we specialise the average scattered wave (3.15) to the material geometry, then use the average boundary conditions (4.8) to restrict the wavemode f 1,n (r 1 , λ 1 ). With the wavemode we can calculate both the average scattered and transmitted wave, although transmission requires some extra steps [42].
In this section, we calculate the average scattered wave from a sphere filled with particles. To our knowledge, the sphere case has never been analytically calculated in all its details, though there have been approximate methods [46] and numerical methods that simulate a large number of configurations [39,40,38].

The average boundary conditions
We assume that all particles are confined in a sphere of radius R which implies that the particle origins R 1 = {r 1 ∈ R 3 : |r 1 | ≤ R − a 1 }, and let the centre of the sphere be the origin of the coordinate system for r 1 . Choosing a simple geometry allows us to explicitly calculate the average boundary conditions (4.4). Assume that the k p have been determined from (5.15) and that the F p,nn1 , up to a multiplying constant α p , have been determined from (5.7). The results below can be used to completely determine, or just restrict, the α p .

Sparse boundary conditions
We find that numerically (6.3) has a unique solution for the α p when using only one effective wavenumber k 1 . This occurs even though (6.3) often has far more equations than unknowns α p . That is, the number of equations is larger than the number of effective eigenvectors P . Nonetheless, both the left and right-hand side are sparse, and mostly filled with zeros. This indicates that it is possible to transform (6.3) into a smaller equivalent system to determine the α p .

Average scattered field
Assume the particles are confined in a spherical region of radius R. For the particles to fit in this region the particle origins need to be contained within R 1 = {r 1 ∈ R 3 : |r 1 | ≤ R − a 1 }, where we let the centre of the sphere be the origin of the coordinate system for r 1 . In this case, by taking |r| > R we can use the wave representation (4.1), Green's second identity, and then (5.1), to reduce the average scattered wave (3.15): where we used (B.1-B.3) to write u n (kr − kr 1 ) = n n2 c n nn2 v * n2 (k * r 1 )u n (kr) followed by the orthogonality of the spherical harmonics, and we defined and M (x, y) = xj (x)j (y) − yj (x)j (y). The F n (λ 1 ) are the scattering coefficients of the whole spherical region (3.16). In conclusion, substituting the above into (3.15) leads to u sc (r) = n u n (kr) S F n (λ 1 )n(λ 1 )dλ 1 , r > R, (average scattered field) (6.7) where we averaged over particle rotations.
To help piece together the equations, we provide Algorithm below to calculate the average scattered wave above when only one effective wavenumber k 1 has a smaller imaginary part than the others as shown in Figure 5.
6. Algorithm -scattering from a sphere filled with particles.
This algorithm assumes one effective wavenumber has a imaginary part significantly smaller than the others. 1. Use Algorithm 4 to calculate the wavenumber k 1 with smallest imaginary part. 2. Choose a truncation 1 ≤ L 1 for the n 1 in F p,nn1 based on how the M 1 in (3.16) decay. The range of all other indices can now be determined from these truncations and the properties (B.6). 3. Calculate the multiple eigenvectors F p,nn1 (λ 1 ) of k 1 by solving (5.7). From (4.1) and (5.1) we now have that f n (r 1 , λ 1 ) = n1 v n1 (k 1 r 1 ) P p=1 α p F p,nn1 (λ 1 ), where the α p need to be determined. 4. Solve (6.3) by: This causes the number of unknowns α p to be equal to the number of coefficients g n , which leads to a unique solution α p .

setting
5. Finally, calculate the scattering coefficients of the whole spherical region (6.6).
a For azimuthal symmetry this becomes L = P − 1

A plate filled with particles
In this section, we calculate the average reflected wave from a plate region R = {r ∈ R 3 : Z 1 ≤ z ≤ Z 2 } filled with particles, and the average transmitted wave that passes through to the other side of the plate. In this case, the particle origins are confined to the region where Z 2 needs to be large enough so that Z 1 + a 1 + a 12 < Z 2 − a 1 − a 12 for every particle radius a 1 and minimum inter-particle distance a 12 . We now assume that the incident wave is a plane wave with wave-vector k = (k x , k y , k z ), where k z > 0 (the incident wave impinges the plate from below). Further, by planar symmetry (3.21) we can also assume that k p x = k x and k p y = k y .

Average transmission
We start with the transmitted field using the wave representation (3.14), (3.15), (4.1), and (5.10). By assuming that z > Z 2 , which is the side of the plate where the transmitted wave will appear, we can use Green's second identity, as we did in Section 4.1, to reduce the average (total) transmitted wave: u(r) = e ik·r + u sc (r) = e ik·r + n,p S n(λ 1 )F p,n (λ 1 ) R1 e ikp·r1 u n (kr − kr 1 ) dr 1 dλ 1 . (7.1) By using Green's second identity we can reduce the integral in r 1 to surface integrals: where we used a change of integration variable from r 1 → (r 1 − r), u n (−kr 1 ) = u n (kr 1 )(−1) , and the definition of L n which is given by (F.1). As both Z 2 − a 1 − z and Z 1 + a 1 − z are negative real numbers, we can use the result (F.6) to evaluate the above and reach u(r) = T plate e ik·r , (average transmitted field) (7.3) where and k p = (k x , k y , k pz ) due to planar symmetry (3.21), for some complex number k pz , and we used k 2 − k 2 p = k 2 z − k 2 pz .

(7.4)
By using Green's second identity we can reduce the integral over R 1 as done in (7.2). As the plate is thicker than any one particle, we then have that Z 2 − a 1 − z > 0 and Z 1 + a 1 − z > 0 for z < Z 1 , which allows us to pick the positive argument in (F.5), evaluate (7.2), and reduce (7.4) to u sc (r) = R plate e ik ref ·x , (average reflected field) (7.5) where R plate is the reflection coefficient: we define k ref = (k x , k y , −k z ), used k p = (k x , k y , k pz ) due to planar symmetry (3.21), and that k 2 − k 2 p = k 2 z − k 2 pz .

The average boundary conditions
When using only one effective wavenumber k 1 the equation (4.8) can be used to fully determine the field (5.10), like a boundary condition. Note we can use planar symmetry (3.20) and the form (5.10) because both the incident wave and material region share a planar symmetry.
The first step is to simplify (4.4): To explicitly calculate the above integrals, we use the translation matrices in Appendix B, followed by changing the integration variable to r = r 2 − r 1 and then using the definition (F.1) to obtain where factor (−1) 1 appeared when substituting u n1 (kr 1 − kr 2 ) = (−1) 1 u n1 (kr 2 − kr 1 ) in the integrals.
We can use the formula (F.6) to easily calculate L n1 by noting that Z 2 − a 2 − z 1 > 0 and Z 1 + a 2 − z 1 < 0. These inequalities are a result of r 1 ∈ R 1 (a 12 ) which implies that Z 1 + a 1 + a 12 ≤ z 1 ≤ Z 2 − a 1 − a 12 . Substituting (F.6) into (7.8) then leads to where we replaced (−1) m1 = (−1) m −m , (−1) 1 = (−1) + , and c n nn1 = (−1) +m c n1( ,−m)n by using the properties of c n nn1 shown in Appendix B. These replacements allow us to simplify (7.9) by applying the contraction rule (B.12) and some rearrangement to reach: where and we also used k p x = k x and k p y = k y . Substituting the above into (4.8) leads to where we used that which holds for incident plane waves with the coefficients (2.5) and can be shown by using (A.2) and the contraction rule (B.12). For (7.12) to hold for every r 1 ∈ R 1 (a 12 ) leads to two equations: one for the term multiplying e i(kx,ky,−kz)·r1 and another for the terms multiplying e ik·r1 . These two equations can be written in the form: where I 1 and I 2 are given in (7.11).
For a finite plate, both (7.13) and (7.14) need to be enforced to restrict the F p,n (λ 1 ). If the sum over p has only two terms, using one a forward propagating and the other a backward propagating mode, then these equations can be used to obtain a unique solution for the F p,n . This typically occurs when using only one effective wavenumber k 1 . For reflection from a halfspace, only (7.13) should be enforced, which is the multi-species three dimensional version of [42,Equation (20)].

Numerical results: plane-wave incident on a particulate sphere
This paper is the first, to our knowledge, to provide analytic solutions for the average wave scattered from particles within a spherical region R as given by (6.3) and (6.7). The methods used previously [38,46] have approximated the scattered field by assuming that the ensemble averaged sphere behaves like a homogeneous sphere occupying the region R with some effective properties. For these reasons, in this section we numerically compare these approaches.
For all the results below, we avoid combining a high particle volume fraction with a high frequency, as this regime triggers multiple wavenumbers with low imaginary parts, as shown in Figure 4. Whereas the calculations below rely on equation (6.3) giving a unique solution, which only occurs when using just one effective wavenumber. Using only one wavenumber is an excellent approximation when its imaginary part is much smaller than all the other wavenumbers, as shown in Figure 5. See [24,26] for details on how to calculate reflection and transmission when multiple effective wavenumbers have a small imaginary part.
mass density wave speed radius volume % Solid particles ρ s /ρ = 10 c s /c = 10 a s /R = 1/20 ϕ s = 15% Void particles Table 2: Material properties used for the numerical experiments. Note that R is the radius of the spherical region R.
The analytic scattered field. For each (angular) frequency ω, we calculate the effective wavenumber k 1 with the smallest imaginary part. We then follow Algorithm 4 and Algorithm 6 to calculate the scattered field. As a reminder, this method does not assume the spherical region behaves like some homogeneous sphere; instead these results are from careful homogenisation of all the scattered waves.
Two different homogeneous spheres. We can approximately calculate the scattered wave from the sphere R by assuming that R is filled with some homogeneous material. Below, we choose two different ways to approximate the density and sound speed of this homogeneous material we use to fill R. Note that there are many possible chooses for the density and sound speed and no clear "best choice". R R/λ = 0.2 R/λ = 0.7 Figure 6: An illustration of the size of the spherical region R, the size of particles inside, and the typical length of the incident wavelengths λ used for the results in Figures 7 -9. Note that R is the radius of R, and numerical parameters used in this section are in Table 2. 1. Hom. Low Freq.: we assume the sphere has the effective density ρ * and effective bulk module β * given by (5.23), which results in the sound speed c * = β * /ρ * . 2. Hom. Complex k 1 : we assume the sphere has the same complex wavenumber k 1 used for the analytic solution, which then implies it has sound speed c 1 = ω/k 1 . For the effective density, we again choose ρ * given by (5.23).
After choosing one of these approximations, we can calculate the scattering coefficients F n by using the T-matrix (2.8).
Taking the origin to be the centre of the spherical region R, we can then express the scattered field in the form for a sphere of radius R. For the analytic solution we have F n = S F n (λ 1 )n(λ 1 )dλ 1 from (6.7), and for the numerical results we approximate the integral as a sum.
Frequency sweep. We begin with a frequency sweep and use the particle properties given in Table 2. For each frequency, we calculate the scattering cross section for the three methods described above: the analytic and the two homogeneous spheres. The results are shown in Figure 7.
We define the non-dimensional scattering cross section by [55] (scat. cross section) = 1 2π(kR) 2 n |F n | 2 , where the | · | represents the absolute value. The above is dimensionless and the natural way to compare with the geometrical cross section of the sphere [32]. In the standard notation, σ s often denotes the scattering cross section, in which case our non-dimensional scattering cross section is equal to σ s /(2πR 2 ).
As expected the three methods converge for low frequencies, as shown in Figure 7. For R/λ > 0.05 the Homogeneous Low Frequency sphere quickly diverges from the other two solutions, and then has far more resonant frequencies. The two methods that use the same effective wavenumber, k 1 , stay closer together, but are significantly different even before reaching R/λ = 0.6. Around R/λ = 0.72 we see that both the analytic and the Homogeneous Complex k 1 methods hit a resonant frequency, but display very different responses. To further investigate this, we plot the full scattered field for the three methods in Figure 8. These fields show contour maps for the slice y = 0. The main difference between the methods is that the analytic solution has a weaker scattered field and also has a smaller shadow region.
Varying the particle volume fraction. The effects of multiple scattering between particles vary significantly with the volume fraction of the particles as shown in Figure 9. To produce these results we used a fixed frequency which corresponds to R/λ = 0.133. This frequency was chosen as it is relatively low and is the local minimum of the analytic scattering cross-section in Figure 7. As this is a relatively low frequency, it avoids the need to use multiple effective wavenumbers even for large volume fraction, as described in the beginning of this section.
For a moderate volume fraction, Hom. Complex k 1 is qualitatively a good approximation (except close to resonant frequencies) as shown in Figure 7. However, when increasing the volume fraction we see a clear drift between Hom. Complex k 1 and the Analytic method in Figure 9. Again we notice that Hom. Low Freq. has more resonant frequencies, and they are more extreme.   Table 2. The Hom. Complex k 1 exhibits a strong resonance, with the peak climbing above 13, which is not shown to avoid zooming out too far. At a similar frequency, the Analytic solution exhibits the opposite, where scattering is very weak. Note we did not show the range 0.25 < R/λ < 0.6 as it is less interesting.

Discussion
Much has already been understood about a plate, or half-space, filled with a random mix of particles, including how to calculate, and make sense of, the effective wavenumbers, reflection, and transmission [57,56,6,26,24,35,51,44]. These results are now used to probe emulsions, colloids, and slurry [7,20] with sound, and planetary systems with light [44], among other applications.
A question that remained was: how to make sense of other regions R not shaped like a plate? For example, like a droplet filled with a particulate.
Effective wavenumbers. One milestone of this paper was to show that any region R filled with the same particulate material will have the same effective wavenumbers, and these effective wavenumbers are given by solving (5.15) or (5.16), with the low frequency properties given by (5.23). In fact, (5.16) is simpler to solve then the dispersion equations previously presented in the literature.
A key step we used was to represent the average wave as a sum of wave potentials each with a different effective wavenumber, as shown by (4.1). This representation is useful because the sum converges. This happens because most of the wave potentials decay rapidly due to their wavenumbers having a large imaginary part as shown in Figures 4 and  5. For more details see [24].  Figure 7 shows the scattering cross-section for these three methods over a large frequency range.  Multiple effective wavenumbers. In this work, we concentrated on scenarios where all the effective wavenumbers, except one, lead to wave modes which decay rapidly. That is, we make use of only one effective wavenumber. This scenario, which occurs for most frequencies and particle properties simplifies the equations. See [26] for an example where many effective wavenumbers are used. It remains an open challenge to find a simply way to incorporate all effective wavenumbers for scenarios such as a sphere filled with particles.
The ensemble wave equation. Our results have enabled us to take effective wavenumbers from a halfspace, or a plate, and use them to calculate the average scattered wave from a sphere filled with particles. To our knowledge, we are the first to provide a clear first-principals approach to achieve this. Beyond the examples we present in this paper, like a sphere filled with particles, our ensemble wave equation (4.7) and ensemble boundary conditions (4.8) can be used to calculate the average field for regions of any shape. Though, depending on the shape, this may require considerable work.
Numerical results. To both demonstrate that our method can completely describe the average scattered field, and to compare with previous approaches, we present some numerical results for a sphere filled with particles in Section 8. We compared our method with approaches which assume the region R is made of some homogeneous material with effective properties. As expected, the different methods converge for low-frequency, as shown in Figures 7 and 9, though there are significant differences for finite frequencies. For one specific frequency, the difference between the methods is illustrated by a field plot in Figure 8.
Validation. The next natural step is to validate our models. Numerical validation would be ideal, as there are robust numerical methods for multiple scattering [23,44,38]. Numerical methods can also clarify the assumptions used in the modelling, such as the choice of pair-correlation and the quasi-crystalline approximation. However, a major issue, that has prevented substantial validation, is that these numerical methods have struggled to simulate an infinite halfspace or infinite plate required by most of the available theoretical predictions [25,8,9]. Now, with our framework, numerical validation for finite sized sphere filled with particles should be straight-forward. This will allow a clear way to verify the statistical assumptions used, and the range of their validity.
Electromagnetism and Elastodynamics. Our framework deals with the scalar wave equation. There exist in the literature clear routes on how to extend effective wave theory from the scalar version to elastodynamics [10], thermovisco-elasticity [37,49] and electromagnetism [50,14,13], though each requires extra algebraic manipulation. In light scattering, it is far easier to measure the average of the scattered intensity [44], though it requires the average of the scattered field, which is what we calculate in this work. Extending our framework to calculate the average intensity should enable accurate models for scattering from spheres and other compact objects.

Acknowledgements
Gerhard wishes to gratefully thank the UK Acoustic Network funded by EPSRC (EP/R005001/1) for a generous travel support which made it possible for Gerhard to visit Sheffield in the fall 2019. The authors would also like to acknowledge the late Michael Mishchenko for putting the authors in touch, which ultimately led to this paper, and for his amazing contribution to the field of scattering. The authors are also thankful to Thomas Wriedt for organising the Bremen Workshop on Light Scattering.
The spherical harmonics satisfy the parity relation and complex conjugate , and Y n (r) are orthonormal over the real unit sphere Ω, that is Plane waves can be expanded in terms of spherical harmonics by using: where both x and y can be complex vectors, and we use the dot product to mean (x 1 , x 2 , x 3 ) · (y 1 , y 2 , y 3 ) = x 1 y 1 + x 2 y 2 + x 3 y 3 with no conjugation.

B Translation matrices
The translation properties of the spherical waves are instrumental for the formulation and the solution of the scattering problem of many individual particles. These translation properties are well know, and we refer to, e.g., [3,19] for details. Some of their properties are reviewed in this appendix and a simple proof of these matrices are given in the supplementary material.
Let r = r + d, then the translation matrices for a translation d are Translation in the opposite direction is identical to the Hermitian conjugate of the translation matrices [48], i.e., The translation matrix V nn (kd) is identical to U nn (kd) but with h (1) λ (k|d|) replaced with j λ (k|d|). Notice that the translation matrices V nn (kd) and U nn (kd) have the form where the summation over the multi-index The explicit values of the coefficients c nn n1 are, see the supplementary material where the last is the contraction rule, or the linearisation formula [41]. For real θ and φ the linearisation formula can be deduced by multiplying both sides of (B.12) by Y * n2 (r), then integrating overr, and applying the definition (B.4).

F Integrals of spherical and plane waves
When dealing with effective plane-waves, we need to evaluate the following integral: L n (z) = These integrals converge when Im k ≥ |Im k p x | + |Im k p y |, where k p = (k p x , k p y , k p z ). This inequality holds when using planar symmetry (3.21), which together with k = (k x , k y , k z ) implies that k p x = k x and k p y = k y .
Substituting the above representation into L n (z) leads to i(k p z − q z )Y n (q)e i(kp+q)·r dq x dq y kq z dx dy, z > 0, i(k p z + q z )Y n (q)e i(kp−q)·r dq x dq y kq z dx dy, z < 0.
then integrating in q x and q y we get L n (z) = Y n (k + p ) 2πi i kq z (−1) m (k p z − q z )e i(kp z +qz)z , z > 0, (−1) (k p z + q z )e i(kp z −qz)z , z < 0, (F. 5) where k + p = (k p x , k p y , q z ) and q z = (k 2 − k p 2 x − k p 2 y ) 1/2 , and we used that Y n (k − p ) = (−1) m Y n (k + p ) where k − p = (−k p x , −k p y , q z ).
In most cases where we use plane-waves, we will assume the material occupies the region R 1 = {z > a 1 : r ∈ R 3 }.
In this case we have that k p x = k x and k p y = k y , due to planar symmetry (3.21), which implies thatk + p =k and q z = (k 2 − k 2 x − k 2 y ) 1/2 = k z . Substituting these results in (F.5) then leads to L n (z) = Y n (k) 2πi i kk z (−1) m (k p z − k z )e i(kp z +kz)z , z > 0, (−1) (k p z + k z )e i(kp z −kz)z , z < 0, (F. 6) where k pz is the z component of k p . The case z < 0 gives the same result obtained in [35,Equation B.5].