The average transmitted wave in random particulate materials

Microwave remote sensing is significantly altered when passing through clouds or dense ice. This phenomenon is not unique to microwaves; for instance, ultrasound is also disrupted when traversing through heterogeneous tissues. Understanding the average transmission in particle-filled environments is central to improve data extraction or even to create materials that can selectively block or absorb certain wave frequencies. Most methods that calculate the average transmitted field assume that it satisfies a wave equation with a complex effective wavenumber. However, recent theoretical work has predicted more than one effective wave propagating even in a material which is statistically isotropic and for scalar waves. In this work we provide the first clear evidence of these predicted multiple effective waves by using high-fidelity Monte-Carlo simulations that do not make any statistical assumptions. To achieve this, it was necessary to fill in a missing link in the theory for particulate materials: we prove that the incident wave does not propagate within the material, which is usually taken as an assumption called the Ewald–Oseen extinction theorem. By proving this we conclude that the extinction length—the distance it takes for the incident wave to be extinct—is equal to the correlation length between the particles.


Introduction
Most materials, at some length scale, are formed of a random configuration of smaller particles. Consider particles in powder for pharmaceuticals, grains in the sand, oil droplets in emulsions, and aggregates in solid composites. The wide number of these particulate materials, and engineering applications, make it worth while to develop methods to measure these materials and design them intelligently.
Background. When it comes to measurement and characterisation, the main tools are classical waves such as electromagnetic and ultrasonic or acoustic. The governing equations for these classical waves would be well understood if the material itself was known in all its details. Unfortunately, in most cases it is impossible to know in detail the microstructure of the material because it is disordered. In these scenarios, ensemble averaging and statistical assumptions need to be used to obtain solvable systems [16,30,31,44].
How classical waves interact with particulate materials (on average) is well understood within certain limits. In the long wavelength limit, where the particles appear small compared to the wavelength of the incident wave, it is well understood how to calculate effective properties [37]. In the dilute limit, where there is no multiple scattering, Mie theory has led to characterisation methods such as Dynamic Light Scattering and laser diffraction have led to a range of widely used tools 12 .
Pushing the limits. In the cases where multiple scattering is significant, and the incident wavelength is not long (compared to the microstructure), the average wave is not as simple to describe [4,15,26]. This is especially true when using exotic pair-correlations [38,39] and resonant particles [25,34]. To push the theory to these new limits, we need to clearly understand the validity of all the assumptions made. Within this context, we aim to address two significant assumptions that currently remain unanswered.
Multiple effective wavenumbers. Most of the literature assumes there is only one effective wavenumber [8,23,27,36,41,45]. As the medium is isotropic and homogeneous (after ensemble averaging) it seems reasonable to assume that there is only one effective wavenumber k ⋆ for waves travelling in a bulk material (i.e. no waveguide). However, two different theoretical methods [15,47] have predicted that there exist at least two (complex) effective wavenumbers for one fixed frequency. Here we give the first clear numerical evidence of these multiple effective wavenumbers as well as demonstrate that multiple wavenumbers are trigger by particles, and frequencies, that lead to strong scattering.
Incident wave extinction. The second assumption we address is about the presence of the incident wave inside the particulate material. It is often assumed that the incident wave does not propagate, or contribute to the total field, inside the particulate material. This assumption is called the Ewald-Oseen extinction theorem [4,9,22,30,36,41,42]. We address this by providing a proof that for any material geometry, frequency, and incident wave, the Ewald-Oseen extinction theorem holds when deep enough within the material. Specifically, we show that the incident wave does not propagate further than the correlation length between the particles. That is, the extinction length is equal to the inter particle correlation length.
Monte-Carlo. There have been several studies that use Monte-Carlo methods to validate effective wave theory. Examples include comparing Monte-Carlo with: the average scattering from a sphere filled with particles [33,50], and one effective wavenumber from the theory [5,6]. To our knowledge, there has been no Monte-Carlo validation or evidence that more than one multiple effective wavenumber exist. Here, by using precise Monte-Carlo simulations, we provide the first clear evidence that at least two effective wavenumbers, and therefore two effective waves, are present in the transmitted field. The theoretical methods that predict these multiple effective wavenumbers make several statistical assumptions, whereas our numerical simulations make no such assumptions. They therefore provide a clear validation of the theoretical predictions.
Summary of the paper. Below in Section 1.1 we provide an overview of the theory and the results of this paper.
In Section 2 we discuss what the theory predicts for the average wave in a plate filled with random particles, followed by how to easily identify when multiple wavenumbers should appear.
In Section 3 we discuss our Monte-Carlo simulations, which involve simulating waves scattered from tens of thousands of particle configurations, how we verify convergence, and how we clearly demonstrate that there are scenarios where at least two effective wavenumbers appear in the transmitted field.
In Section 4 we provide rigorous derivations that: the incident wave does not propagate in the particulate material, and the average transmitted wave is a sum of waves which satisfy effective wave equations. Our derivations are more general than just for a plate, they hold for a finite region, and in fact for any spatial dimension. So the proof we provide is also valid for three dimensional materials.

Overview of the theory
Consider a harmonic incident plane-wave u inc (x) = e ikx , satisfying Helmholtz equation with wavenumber k, so that u inc (x)e −iωt satisfies a scalar wave equation, where x is the distance of propagation. When this incident wave propagates through a random particulate medium, the general assumption has been that it will be replaced by one effective wave of the form: where k ⋆ is a complex effective wavenumber, for the fixed frequency ω, A ⋆ is the amplitude, and ⟨u(x)⟩ is the ensemble average of u(x) over all possible particle configurations [8,23,27,36,41,45].
In the low frequency limit, when the particles are small relative to the incident wavelength, there is substantial evidence to justify (1). But beyond the low frequency limit, and when using more exotic pair-correlations and distributions for the particles, there is no clear consensus. Two different methods [12,47] suggest a different form for (1) given by where the k p are complex effective wavenumbers. In the next two paragraphs we further explain the form (37). Incident wave extinction. When A inc is not zero, a part of the incident wave remains even when propagating through a random particulate. The work by Martin [28] suggests that this is the case. Assuming A inc = 0 is called the Ewald-Oseen extinction theorem [9,22,36], and it is due to the scattering and absorption of energy by the random distribution of particles. Much of the literature [2] assumes this extinction happens after the incident wave has propagated a certain distance in the material. In this paper we prove that for any: frequency, material geometry, particle distribution, and for two and three spatial dimensions, the incident wave is extinct (A inc = 0 for plane-waves) at a distance equal to particle correlation length.
Multiple effective wavenumbers. Each term in the sum of (2) represents an effective wave with a different wavenumber k p . This is highly unusual for scalar waves at a fixed angular frequency ω. Yet, two different theoretical methods have predicted the existence of at least two (p > 1) complex effective wavenumbers [12,35,47]. The evidence for this unusual prediction, beyond just theoretical, is lacking. In this work we use highly accurate simulations and Monte-Carlo method for circular cylindrical particles to demonstrate that these extra wavenumbers are present for specific frequencies. To our knowledge, this is the first clear evidence of the existence of these multiple effective wavenumbers.
Average reflected field. Beyond just a curiosity, these extra effective wavenumbers can have a significant effect on the average reflected, transmitted, or scattered wave from a particulate material [10,12]. It is also far simpler to calculate the cases where there is only one dominate wavenumber [8,15]. showing when more than one effective wavenumber is needed. The x-axis shows ka, with a being the particle radius, and k being the incident wavenumber. The y-axis is the particle volume fraction. The regions with a light colour in the background, anything less than approximately 0.5 shown in the colour bar, require more than one effective wavenumber. The height of the green curve is the scattering strength of just one particle (given by (18)). Clearly strong scatterers lead to multiple effective wavenumbers.
To understand when these multiple wavenumbers are needed, we have produced a series of phase diagrams. An example is shown in Figure 1 for sound soft particles. The regions with a light colour, anything below 0.5 shown in the colour bar, require more than one effective wavenumber to accurately describe the transmitted field. For more details see Figure 4 and the discussion around that figure. For example, we can see that only one wavenumber is needed for low particle volume fraction ϕ < 0.1, with only one exception around ka = 0.6. The green curve in the figure shows the scattering strength of just one particle by itself. Surprisingly, we see that at the frequencies at which the single particle scatters the strongest (the peaks in the green curve) is also the frequencies at which two (or more) effective wavenumbers are required, and therefore needed to accurately describe wave transmission. These results, and other phase diagrams, are further discussed in Section 3.1.

A plate filled with particles
For a medium that is isotropic and homogeneous, it seems reasonable to assume that there is only one effective wavenumber for one fixed angular frequency ω. However, recently two different theoretical models have predicted that there exist at least two (complex) effective wavenumbers [12,47], as shown in (37).
Here we design a computational experiment to give clear evidence of at least two of these effective wavenumbers. We do this by using a robust numerical method based on high fidelity Monte-Carlo simulations for two dimensional disks.
To describe the material, let P j be the disk occupied by the j-th particle, as represented by the circles in Figure 2 and shown in more detail in Figure 10. Let P = ∪ j P j be the union of all particles. For simplicity we consider circular particles of equal size. In other words, using standard set-builder notation where a is the radius of the particle, and |x| is the length of the vector x.
The particles are restricted inside the plate geometry R, given by where W is the width and H is the height of the plate as shown in Figure 2. Figure 2: Scattering of an incident plane-wave approaching from the bottom, onto one specific configuration of randomly distributed circular cylinders (or particles) Λ. The simulation directly solves the governing equations [29], and the green line shows where the field is measured.
The total field u(r) satisfies a Helmholtz equation which depends on whether r is inside a particle or not: where k = ω/c and k 0 = ω/c 0 are the real wavenumbers of the background and particles respectively. The scalars c and c 0 are, respectively, the wavespeeds in the background and particles.
The simplest scenario to numerically check for effective waves is for planar symmetry. As in this case each frequency has only one mode: the plane-wave. To achieve this, we fill a plate region with a configuration of randomly distributed cylindrical particles, as shown in Figure 2. The particles are identical, except for their positions.

Effective waves for planar symmetry
Here we summarise the results of the theory for plane-wave symmetry. The results here will be compared with a Monte-Carlo method detailed in the next section.
We consider an incident plane-wave of the form: and consider particles in a plate region R with an infinite height, whereas Figure 2 shows a truncated plate with a finite height. The theoretical methods consider an ensemble average of the total field u. To achieve this, we describe one configuration of particles with where r j is the centre position of the particle P j . Naturally, the field u depends on the particle positions. To make this explicit we use u(x; Λ). Next, to calculate the ensemble average, we need to define the probability of all possible particle configurations. To do this, we introduce the joint probability density given by p(Λ). For a brief overview on the probability density function p, see [13,15,23]. The theoretical methods then calculate and predict the ensemble average defined by where the integral is over all possible particle positions, and the fields depend only on the spatial position x as we are considering planar symmetry. It is widely assumed that the average ⟨u(x)⟩ satisfies a Helmholtz equation with a unique effective complex wavenumber k * [8,23,27,36,41,45]. In Section 4.1 we prove that ⟨u(x)⟩ is a sum of several effective waves, but only when x is deep enough within the material. To define what "deep enough" means we need to introduce the particle pair-correlation.
Pair-correlations. We assume that particles are distributed both homogeneously and isotropically, which leads to for an infinite number of particles. For details on pair-correlations see [18]. For a disordered or random configuration of particles we have that where b 12 ≥ a 12 > 2a. For the region a 12 ≤ |r 1 − r 2 | ≤ b 12 the pair correlation can take any values for our calculation below, though we expect g(r) to be continuous in r. In this work we use two different pair-correlations. The first, and the simplest, is called Hole-Correction, which assumes that b 12 = a 12 . The second is called the Percus-Yevick approximation, which more accurately approximates the pair correlation for particles that are distributed according to a uniform random probability, except no two particles can overlap [3,18,21,41]. We use the results from [1] to obtain Percus-Yevick for disks. Figure 3 shows the Percus-Yevick distributions for several different particle volume fractions ϕ. Figure 3: The Percus-Yevick approximation is a pair correlation that represents particles that are uniformly randomly placed, except particles do not overlap. That is, particles do no attract or repel each other. The particle radius a = 1.2 and ϕ is the particle volume fraction.
For disordered media we have that g(r) → 1 as r → ∞. That is, particles become uncorrelated as they are further apart. To prove that ⟨u(x)⟩ is a sum of effective waves we need a slightly stronger assumption: that there is a distance b 12 at which particles are completely uncorrelated as used in (11).
Effective plane-waves. By having a length b 12 at which particles become uncorrelated, we demonstrate in Section 4.1 that for an infinite plate geometry (H → ∞ in (4)) filled with particles where we used planar symmetry as shown in [15]. The A ± p are complex amplitudes, the k p are the complex effective wavenumber, and P is the number of effective wavenumbers. There is an infinite number of wavenumbers P , but according to theoretical calculations only a few are needed for accurate results. For a detailed discussion on these multiple effective wavenumbers see [10,12].
The dispersion equation. To calculate the wavenumbers k p we use the dispersion equation appearing in [10,12,23]. The assumptions needed to arrive at this dispersion equation are shown in Section 4. To summarise, the k p are determined by solving where G ℓ =ˆb 12 a 12 J ℓ (k ⋆ r)H ℓ (kr)(g(r) − 1)rdr, (15) and the term n is the average number of particles per area, r = |r|, δ nn ′ is the Kronecker delta, and H ℓ is the Hankel function of the first kind, while J ℓ is the Bessel function. The term T n is the T-matrix which determines how one particle scatters waves by itself. For circular homogeneous particles in acoustics we have that where γ = ρ o c o /ρc and k o = ω/c o , with ρ o being the mass density of the particles and c o being the wavespeed within the particles. There are infinitely many k ⋆ which solve det M (k ⋆ ) = 0. We denote these solutions as k 1 , k 2 , . . .. The main objective of our Monte-Carlo simulations is to check if the theoretical predictions of the wavenumbers k p are accurate, and to clearly demonstrate that there is more than one effective wavenumber appearing in the Monte-Carlo results. Before doing this, let us first explore the effective wavenumbers predicted by solving the dispersion equation (13).
Effective wavenumbers. We want to identify when the dispersion equation (13) predicts that there is more than one effective wavenumber that has a significant contribution to the average transmitted wave. It is important to understand when this occurs, as it is far simpler to calculate the average field when there is only one effective wavenumber [12,15].
Most of the scientific community at present is also not aware that more than one effective wavenumber can be excited [4,29,30], for just one scalar wave. So we will identify for which parameters we can run a heavy Monte-Carlo simulation, as detailed in the next section, to find clear numerical evidence of multiple effective wavenumbers.  Table 1. The colour is given by (17) where the lighter colours (those above 0.5 shown in the color bar) indicate that more than one effective wavenumber can be excited. The green curve shows the scattering strength of just one particle and is given by (18).  Table 1: Shows the two main particle properties used for the numerical results. Note that sound-soft (sound-hard) particles are strong (weak) scatterers.
According to theoretical results, only one effective wavenumber k 1 is needed when Im k 1 ≪ Im k p for p = 2, 3, . . . . In any other case, more than one effective wavenumber can be excited and contribute to the average transmission [10,12]. Though we note that the form of the incident wave and geometry of the material also affect how the wavenumbers are excited [15].
The first step is to sweep the parameter space by varying the frequency ω and particle volume fraction ϕ, and for each value calculate the effective wavenumbers k p by solving (13). We do this for both sound-hard and sound-soft particles by changing ρ o and c o in (16). The two main particle properties used are shown in Table 1. Next, based on the results shown in [12,10], we can estimate where more than one effective wavenumber is excited by plotting a heatmap where the colour is given by After many failed attempts, we find that the measure (17) is closely related to the scattering strength of a single particle, given by where T n is the T-matrix given by (16) for acoustics. The results of sweeping over frequency and particle volume fraction are shown in Figure 4. We call these figures phase diagrams, as we see sudden shifts from only one effective wavenumber to two or more wavenumbers. The regions with lighter shading correspond to cases where the value of (17) is low, so more than one effective wavenumber is excited. Conversely, the regions with dark shading is where only one effective wavenumber is excited. The green curves shown on-top of the phase diagrams are the scattering strength for just one particle. Clearly we see that a large scattering strength leads to more than one effective wavenumber, which is an important observation, as calculating (18) is far simpler than calculating the wavenumbers k p .
The phase diagrams in Figure 4 show a large region of the parameter space, but to see in more detail when two or more effective wavenumbers are needed it helps to plot the imaginary parts of the wavenumbers k p against frequency, which we do in Figure 5 for a particle volume fraction of 25%. For sound-soft particles (Figure 5a), there are many frequencies ka where two or even three effective wavenumbers have a similar imaginary part, meaning that these wavenumbers can be excited. In contrast, for sound-hard particles (Figure 5b), there is only one effective wavenumber that can be easily excited, as there is one curve, representing k 1 , that has a significantly lower imaginary part than all the others. Im k p (b) sound-hard scatterers Figure 5: Depicts the imaginary part of the wavenumbers k p with respect to the non-dimensional frequency ka. We only show the three wavenumbers with the smallest imaginary parts as these are the only ones which make a significant contribution to the average wave. The lowest curve represents k 1 , which is the easiest to excite. The particles occupy 25% of the material in each case and their properties for sound soft and hard particles are shown in Table (1). In the next section, we use the results presented here to identify specific scenarios where performing a Monte-Carlo simulation would be appropriate. This will allow us to verify whether the predictions of more than one effective wavenumber are indeed accurate.

The Monte-Carlo simulation
A simple way to approximate the average field (9) is to perform a simulation for each particle configuration and then take an average over all configurations. As we are focusing on uniformly distributed particles, the probability density p(Λ) is a constant if particles do not overlap, and p(Λ) = 0 if any two particles do overlap. This allows us to numerically approximate the ensemble average where each Λ s is one randomly sampled configuration of particles within the plate that depends on the parameter s. That is, we first create a configuration of particles Λ s , then simulate the scattered waves u(r; Λ s ), and then we repeat this process for S configurations of particles, until the average shown above converges. As mentioned in previous sections, ⟨u(r; Λ)⟩ should converge to a sum of plane-waves as shown by (12).
The setup for our Monte-Carlo simulations is shown in Figure 2 (note the x-axis is the vertical axis in the figure), and we used an incident wave of the form u inc (x) = e ikx . To calculate the field u(r; Λ s ) for each configuration we use a multipole expansion for each particle, together with addition translation matrices, to solve the boundary conditions [29]. Specifically, we solve [12, Equation (2.9)] for each s, calculate (19) and evaluate the field at ⟨u(x, 0; Λ s )⟩ with the x values satisfying (12). This method accurately approximates the exact solution when increasing the truncation order of the multipole expansion until reaching convergence.
An initial numerical investigation revealed that it was computationally feasible to simulate a finite number of particles within the region 0 ≤ x ≤ 20, so the width of plate W = 20 and height H = 400. For details on the methodology of our Monte-Carlo simulations, see Appendix A.

The Monte-Carlo results
Based on the results of the previous section, we choose several cases to simulate the average field with a Monte-Carlo method. We want to identify cases where two or even three wavenumbers can be excited for frequencies as low as possible, because increasing the frequencies leads to a significantly larger computational cost for the Monte-Carlo simulations. Which is why we only performed Monte-Carlo simulations for three different frequencies, beyond the low frequency limit. See Appendix A for details, and descriptions on computational cost.  (13) predicts there is only one effective wavenumber k p with a lower imaginary part, and therefore this is the only wavenumber that can be excited. The green and red triangles represent effective wavenumbers predicted by the (13) when using either the Hole-Correction or Percus-Yevick pair-correlation. The purple dotted points represent the effective wavenumber which best fits the Monte-Carlo simulations when using the formula (12).
One effective wavenumber. From Figure 5b we see that for sound-hard particles there is a broad range of frequencies where there is only one effective wavenumber with an imaginary part which is far smaller than all others. This means that is easy to excite this wavenumber, but very difficult to excite the others. This message is also confirmed by the phase diagram shown in Figure 4. To verify there is only one effective wavenumber, we perform Monte-Carlo simulations for ka = 0.36 for sound-hard particles, and fit the formula (12). The results shown in Figure 6b confirm that there is only one wavenumber present in the Monte-Carlo simulations.
There is a small discrepancy between the effective wavenumbers from the dispersion equation (13) and the fitted effective wavenumbers, which could be due to errors introduced in the Monte-Carlo simulations due to truncating an infinite region.
On the other hand, for sound-soft particles we see from Figure 5a that only for lower frequencies 0 ≤ ka ≤ 0.2, only one effective wavenumber with a smaller imaginary part exists. Monte-Carlo simulations, shown in Figure 6a, confirm that there is only one effective wavenumber as predicted from the theory.
Fitting for several effective wavenumbers. Increasing the frequency for sound-soft particles leads to many frequencies where two, or more, effective wavenumbers have a lower imaginary part. For example, for 0.25 < ka < 0.63 there are two wavenumbers with a smaller imaginary part, which means it is possible to excite two effective wavenumbers in this frequency range. To exemplify, we choose two different frequencies to perform Monte-Carlo simulations: 1) ka = 0.36, where we aim to excite two effective wavenumbers, and 2) ka = 0.62 where we are aim to excite three effective wavenumbers. Figure 7: The graphs above compare the average field from a Monte-Carlo simulation for particles in a plate, as shown in Figure 2, with two types of fitted waves. The frequency ka, particle properties and volume fraction ϕ used are shown below each figure. The Dominant Wave is the result of fitting for just one effective wavenumber when using the formula (12), and is currently believed to be accurate by most working in the field. We see here that it is not possible to fit for just one wavenumber. The Fitted wave is a result of fitting the formula (12) with P = 2 for the top two graphs and P = 3 for the bottom two, where the Extended Wave shows what the formula (12) predicts outside of the fitted region. Figure 7 shows how the formula (12) fitted to the Monte-Carlo data clearly matches the Monte-Carlo data, whereas the graphs on the right show how fitting for only one effective wavenumber, termed the "Dominate wave", does not match the Monte-Carlo data. The two panels on the left also show extending the fitted wave, by using the formula (12), beyond the fitted region can still closely follow the Monte-Carlo results, but is expected that as the extended field gets closer to x = 0, the match with the Monte-Carlo field will get worse due to a boundary layer where many effective wavenumbers become significant [10,12]. The presence of this boundary layer, and our theoretical results in Section 4, have all guided how best to perform, and fit to, the Monte-Carlo results as discussed in Appendix A.
Avoid overfitting. Having an accurate fit does not necessarily give strong evidence that the formula (12) is correct. This is specially true for higher frequencies where the Monte-Carlo simulation has a higher standard error of the mean, shown by the shaded yellow region. A larger standard error means that there is a range of effective wavenumbers which can fit this data, and still be closer to the Monte-Carlo results than the standard error. Further, fitting a sum of plane-waves, such as shown by (12), can lead to over-fitting when using too many wavenumbers, and even become ill-posed. We explain more details on this in the Monte-Carlo methodology in Appendix A.
To verify that we avoid over-fitting, and check if small changes in the field lead to large changes in the fitted wavenumbers, we develop a method, called the projection method, which fits the formula (12) for all possible combinations of the effective wavenumbers in a region. Details are given in Appendix A.
At least two effective wavenumbers. The result of fitting for two effective wavenumbers for ka = 0.36 is shown in Figure 8. The main conclusion is that this is the first clear evidence that there are two complex effective wavenumbers by using Monte-Carlo simu-lations which are highly accurate. In more detail: the figure shows that there are two separate effective wavenumbers: one within the blue dash curve on the left and the other, necessarily, within the blue dashed curve on the right. Wavenumbers within these dashed regions lead to fitting errors which are less than the error committed by the Monte-Carlo simulation. We also see the dispersion equation (13), together with the Percus-Yevick pair correlation, predicts at least one wavenumber within the blue dash curve. The predicted wavenumbers when using the Hole-Correction pair-correlation is also shown. The dispersion equation (13) predicts an infinite number of effective wavenumbers, so we present only the two wavenumbers with smallest imaginary part. Figure 8: Shows how two wavenumbers are needed to fit the formula (12) to the Monte-Carlo results, for the particle properties ρ 0 = c 0 = 0.30, ϕ = 25% and frequency ka = 0.36. When using the two best fits, shown by Projection method, we obtain the fitting shown in Figure 7a. The density plot shows what regions of complex wavenumbers that best fit the Monte-Carlo results. When using one wavenumber k 1 in the dashed blue curve on the left, there exists another wavenumber k 2 within the dashed blue region on the right that together to a fitting error which is smaller than the error of the Monte-Carlo simulation (the standard error of the mean). For more details see Appendix A.
Three effective wavenumbers. The Figure 9 shows in steps our attempt to identify three effective wavenumbers for ka = 0.62. The summary is that the amplitude A 3 for the third wavenumber in (12) is too small, |A 3 | ≈ 0.004, and as a result, we can not reliable say whether the Monte-Carlo simulations show that there are indeed three effective wavenumbers. This is because the expected errors from the Monte-Carlo simulation are on the order of around 0.003. However, we do conclude again that one effective wavenumber is not enough, at least two are needed.
In more detail, the top two graphs of Figure 9 show the result of fitting the formula (12) with P = 2 to the Monte-Carlo results. Although there are clearly choices of two wavenumbers which fit the data well, there are no possible choices which lead to fitting errors that are less than the standard error of the mean, as shown by Figure 9a which uses the two wavenumbers with the best fit. The bottom two graphs Figure 9 show how the fitting errors decrease when adding a third wavenumber, i.e. using P = 3 in the formula (12). By using P = 3, instead of P = 2, the fitting error decreased from 6% to 3.8%. Although the fitting error does decrease to below the standard error of the mean, it is only a small decrease, and close to the errors inherent in the Monte-Carlo data. We note that the two Projection Method wavenumbers shown in Figure 9b are equal to two of the Projection Method wavenumbers shown in Figure 9b. Figure 9: The top two graphs show the result of using two effective wavenumbers in the formula (12) to fit to the Monte-Carlo results for the particle properties ρ 0 = c 0 = 0.30, ϕ = 25% and frequency ka = 0.62. The graph 9b shows a density plot over the effective wavenumbers, with light regions indicating that those wavenumbers better fit the data. However, all possible choices of two wavenumbers lead to fitting errors which are less than the standard error of the mean of the Monte-Carlo simulations. The two wavenumbers with the best fit are denoted by the Projection method, and lead to the field shown in the graph 9a. The bottom two graphs use three effective wavenumbers in the formula (12) to fit to the Monte-Carlo results. In this case, we find 4 sets of wavenumbers, all close to each other, that have a fit error less than the standard error of the mean. The result of using the three wavenumbers with the best fit is shown in 9c. However, in 9d we see that there are many choices for the wavenumbers which lead to small fitting errors. In particular, the Projection method wavenumber with the smaller imaginary part is sensitive to small changes in the the Monte-Carlo results.

Deducing the average transmitted wave
In [15,12,10] the authors demonstrated that there exists several effective wavenumbers, however it was not clear how these appear in the average transmitted field. Taking inspiration from [28], the average transmitted field should be a sum of waves, each with a different effective wavenumber. In this section we demonstrate this for any incident field and material geometry. We use the same notation given in [15] and combine the methods shown in [15] and [28] to show that the average transmitted field is a sum of the incident field plus several effective fields.
We note that although the paper [15] is written for three dimensional particles, the results that lead up to [15,Section 5] are valid for any dimension as long as we appropriately define the spherical waves u n and v n and the translation matrices V nn ′ and U nn ′ . In the case of two dimensions and scalar waves these terms are v n (kr) = J n (kr)e imθ , u n (kr) = H n (kr)e imθ , V nn ′ (r) = v n−n ′ (r), and U nn ′ (r) = u n−n ′ (r), where J n and H n are the Bessel function and Hankel function of the first kind, and (r, θ) are the polar coordinates of r. In our calculations here, and below, we will use the general notation rather than substitute the specific form for the two dimensions as shown above. This way, the proofs we present are valid for any dimension. For just one configuration of particles, the way we represent the total field u(r) depends on whether r is inside a particle or not. We choose to write the field in the form [19,20,28] u(r) = u inc (r) + u sc (r), for r ∈ R \ P, where u inc (r) is the incident wave, u sc (r) is a sum of all the scattered waves, and u j in (r) is the field inside particle j. We use P j to denote the region occupied by particle j, whereas P = ∪ j P j is the union of all particles. See Figure 10 for an illustration. Figure 10: A two-dimensional region R filled with equal-sized disks P j , which represent the particles. The region P = ∪ j P j depicts the region inside the particles while the shaded region (R \ P) depicts the region outside the particles. The region R completely contains all the particles, while the region R 1 contains only the particle centres. As the particles are at least one radius a away from the boundary of R, we have that R 1 is smaller than R. Note that as the particles do not overlap, we have that R \ P = R \ P ℓ − j̸ =ℓ P j for every ℓ.
The sum of the scattered waves u sc (r), shown in (22), is given by The field inside the j-th particle can be written in terms of a regular radial waves expansion: and, as a reminder, T n is given by (16). For a light introduction on the T-matrix and multiple scattering see this note and this note. The ensemble average of any field f, and the conditional ensemble average, is defined as ⟨f⟩(r 1 ) =ˆR J−1 1 f p(r 2 , . . . , r J |r 1 )dr 2 · · · dr J , where R J 1 represents that the integration domain of all J integrals is R 1 , and p(r 1 , r 2 , . . . , r J ) is the probability density of having particles centered at r 1 , r 2 , · · · , r J , while p(r 2 , . . . , r J |r 1 ) is the conditional probability density and can be defined as p(r 2 , . . . , r J |r 1 ) := p(r 1 , . . . , r J ) p(r 1 ) .
To calculate the ensemble average of the transmitted field, it is helpful to write the field in the following form where By taking the ensemble average of (29) we obtain ⟨u(r)⟩ = ⟨u inc (r)χ R\P (r)⟩ + ⟨u sc (r)χ R\P (r)⟩ + J⟨u 1 in (r)χ P 1 (r)⟩, where each term of the sum in (29) is the same after ensemble averaging because they are all integrated over the same domain, which means that all particles are indistinguishable from each other. To calculate (31) we need to make a few assumptions. Isotropy and homogeneity assumption. We assume isotropy and homogeneity, which means that p(r 1 ), the probability density of one particle being at r 1 , is a constant. Then, because the integral of p(r 1 ) over R 1 has to equal 1, we conclude that where |R 1 | is the volume of R 1 , and n is the number density of particles defined by with ϕ being the particle volume fraction, and |P 1 | being the volume of a particle. Isotropy and homogeneity also imply that [40] the pair-correlation g had the form (11). For convenience, we now use the form p(r 2 |r 1 ) = n J − 1 g(|r 1 − r 2 |).
The term (J −1) appears now, rather than J, due to there being a finite number of particles, and it ensures that g(|r 1 − r 2 |) → 1 when particles become uncorrelated, as confirmed by [18,Equation (8.1.2)]. For more details on the pair-correlation see [18]. Correlation distance assumption. We can only resolve the integrals appearing in (31) when r satisfies min that is, when the distance of r to the boundary ∂R 1 is greater than b 12 + a. The distance b 12 appears in (11) and is sometimes called the correlation length. Our analysis shows that when the above is true, the transmitted field (31) is a sum of effective waves, and the incident wave is no more.

Transmitted effective waves
In this section, we use the effective wave assumption to demonstrate that (36) is composed of functions which satisfy effective wave equations. Below we show that where ∇ 2 w inc (r) + k 2 w inc (r) = 0 and ∇ 2 w p (r) + k 2 p w p (r) = 0.
We also show that w inc (r) = 0 when (35) holds, and when using the effective boundary condition which is deduced from first principals in [15]. That is, the incident wave is not present inside the material 3 .
For what follows, to keep the notation concise, we define the ball region using standard set-builder notation The first term in (36) is given by (1 − ϕ)u inc (r), which clearly contributes to the term w inc (r) in (37), because the incident wave satisfies the Helmholtz equation (38) 1 .
The second term J⟨u 1 sc (r)χ R\P 1 (r)⟩ is more involved. In Appendix B.3 we prove this term has the reduced form (70). To calculate this term, first we use (23) together with (27) to obtain ⟨u 1 sc (r)⟩(r 1 ) = n ⟨f n ⟩(r 1 )u n (kr − kr 1 ).
Following the method shown in [15], we use a series expansion of the form where each f p,n satisfies a Helmholtz equation: ∇ 2 f p,n (r 1 ) + k 2 p f p,n (r 1 ) = 0, and the k p and f p,n are determined from the governing equation of ⟨f n ⟩(r 1 ). In [15] it is shown how such a series solution can be used for any material region, where [10] proofs the above form for plane-waves.
As shown in [15,Section 4], we can use Green's second identity to writê where from which we can see that ∇ 2 I p,n (r) + k 2 I p,n (r) = 0 and ∇ 2 J p,n (r) + k 2 p J p,n (r) = 0.
So, clearly, I p,n (r) contributes to w inc (r), while J p,n (r) contributes to w p (r) in (37).
The third term in (36) is given by (76) in Appendix B.3.1, which after using (41) and (40) becomes where G(x 1 ) is defined just below (76), although it is not required for our goals here. The first of the above integrals is analogous to (43), so leads to terms of the form (37). The second of these integrals only has r dependence in f p,n (r−x 1 ), and therefore contributes to w p (r) in (37).
The fourth term in (36) is given by (62), which after using (24), (41) and the change of variables from r 1 to x 1 = r − r 1 becomes which can only contribute to terms of the form w p (r) in (37).

The average of the incident field
In the previous section we demonstrated that (31) is a sum of terms which satisfy the background and effective wave equations as shown in (37). The term w inc , that satisfies the background wave equation, can be seen as what remains of the incident field. In much of the literature [4,9,18,22,36] it is simply assumed that w inc := 0. This is often called the "extinction theorem", despite it being an assumption. In papers such as [28] that calculate the average field from first principals, it is not clear that w inc = 0. Here we remove any doubt by proving that when sufficiently inside the material, given by condition (35), we have that w inc := 0 for any incident field, any material region R, any frequency, and all types of particles. Using the results from Section 4.1, we collect the terms in (36) that satisfy the background wave equation to obtain Below we show that the right side of the above is zero by using the ensemble boundary condition given by [15,Equation (4.8)] where with r 2 being the variable of integration and ν 2 the normal to the boundary ∂R 1 , and the g n are the coefficients of the incident wave u inc (r) = n g n v n (kr).
The above assumes the source of the incident wave is outside of the region where the particles are [15]. The other assumptions needed to deduce the ensemble boundary conditions (49) are the same assumptions we have used for the calculations in this paper, except the boundary condition is only valid when min To start the demonstration, we multiply both sides of (49) by v n (kr − kr 1 ) and sum over n to obtain nn ′ V n ′ n (kr 1 )g n ′ v n (kr − kr 1 ) + n nn ′ p I p,n ′ n (r 1 ) k 2 − k 2 p v n (kr − kr 1 ) = 0. (53) The above can be simplified by using the fundamental property of translation matrices that Using the property of the translation matrices above and (51), we see that in accordance with [15,Equation (2.3)]. Next, by choosing r such that (35) is satisfied, it is then possible to choose r 1 so that the condition (52) is true and such that |r − r 1 | < |r 1 − r 2 | for every r 2 ∈ ∂R 1 This enables us to use the translation property (54) in (50) to obtain n I p,n ′ n (r 1 )v n (kr − kr 1 ) = Substituting (55) and (56) into (53) leads to Finally, substituting (57) into (48), we conclude the extinction theorem w inc (r) = 0 for r that satisfies (35). That is, there is no term in the average transmitted wave that satisfies the background wave equation.

Conclusions
The initial goal of this work was to find clear evidence that there exists at least two effective wavenumbers in an averaged particulate material. It is highly unusual to have two different wavenumbers for an isotropic homogeneous media supporting only scalar waves. However theoretical works [10,12,15,47,48,49] have predicted the existence of at least two effective wavenumbers, and their presence changes the overall transmitted and reflected waves. Monte Carlo results. To verify the existence of multiple effective wavenumbers we used very precise simulations that calculated scattered waves from different particle configurations and then took an average over the different particle configurations. This turned out to be far more computational expensive then we expected, and required extensive and careful analysis. To summarise, Figure 8 clearly shows that there are two separate wavenumbers that contribute to the field, and that the wavenumbers predicted by the Monte-Carlo method are similar to the wavenumbers predicted by the theory.
When it matters. A natural question that appeared during this work was how to find the parameters that led to multiple effective wavenumbers. That is, for what scenarios will the classical theory that uses only one effective wavenumber [4,23,24,28,30,45] be inaccurate? Previous work [12,15] demonstrated that there is a dispersion equation (13) which provides the effective wavenumbers k p , and that if there is only one wavenumber k 1 with an imaginary part much small than all others Im k 1 ≪ Im k p for p = 2, . . ., then the classical theory will be accurate.
However, solving the dispersion equation (13) can be time consuming, specially for higher frequencies and wide range of parameters. When plotting the regions where multiple wavenumbers appear we saw a clear pattern shown in Figure 1: particles which are strong scatterers lead to multiple effective wavenumbers. In Figure 1 the green curves show the scattering strength of just one particle (18). Finding the parameters that lead one particle to be a strong scatterer is far more practical than solving the dispersion equation (13), and proved to be a surprisingly good measure. To summarise: strong multiple scattering triggers multiple effective wavenumbers.
Resonators. In the field of metamaterials, strong scatterers such as resonators are often used to tailor the overall behaviour of the material [7]. Using this strategy for disordered, or random, particulates will lead to multiple effective wavenumbers, and will complicate how to predict the overall properties of the material. To truly understand the effect of these resonators it is first necessary to plot their dispersion diagrams by solving the dispersion equation (13) with the T-matrix T n depending on the type of particle used. An example of such a diagram is given by Figure 5.
The theoretical results. When deciding how best to sample the transmitted field, we realised that within the theoretical formulation for ensemble averaging particulates, it was not clear that the transmitted wave is a sum of waves with effective wavenumbers. This led us to derive the missing results and provide a general proof about the incident and transmitted waves.
Proof of extinction. It is often assumed that the average field inside a random material does not contain any remnant of the incident wave. This is called the Ewald-Oseen extinction theorem, but as far as the authors are aware, there is no proof of this conjecture for particulate materials. In this work, we were able to prove this extinction theorem for any particulate (for scalar isotropic waves), any frequency, and material geometry. The proof is given in Section 4, with the final equation that proves extinction being (57). The proof also provides the extinction length: the distance required for the incident wave to travel with the material until it is extinct. We proved that the extinction length is equal to the correlation length plus the particle radius b 12 + a, see equation (11) where these quantities are relative to the pair-correlation.
Proof of transmitted effective waves. The same proof for extinction also served to prove that the average transmitted field is a sum of effective waves, when the distance from the material boundary is greater than the extinction length. The proof is shown in Section 4. We note, that particularly in the field of continuously varying random media [4,46], it is assumed that the average transmitted field satisfies an effective wave equation. By proving this for particulates, from a microscopic approach, we provide a link between the two approaches.
Future avenues. This work shows that using Monte-Carlo simulations to approximate a semi-infinite media, such as a plate, filled with particles, is still computational challenging. We feel that future work focused on validating effective theories for particulates should focus on finite materials (in the computational sense), such as a cylinder filled with cylindrical particles and a sphere filled with spherical particles. There is a theoretical framework to validate against [15]. In terms of theoretical developments, our work has shown a connection between the particulate microscopic approach to effective waves [15,30,31,36,42] and approaches for continuously varying random media [4,46]. That is, we demonstrate the effective wave series used in the macroscopic approach given by (41) does lead to the average transmitted wave being a sum of effective waves, as illustrated by (37). We believe the calculations we provide now paves the way to answer the follow open question: are the two approaches equivalent?
Author contributions AK conceived of study, drafted the manuscript, wrote all the code for the numerical calculations, developed the theoretical calculations, and produced all the figures. ALG helped conceive the study, edited the manuscript, assisted with and verified the theoretical calculations. PSP assisted with the projection method and with the figures related with this method, verified the calculations, and edited the manuscript.

A The Monte-Carlo methodology
Here we present more details on how we performed the Monte-Carlo simulations, and analysed the results. For a reference on Monte-Carlo methods we refer to the book [43].
To compute the ensemble average wave (9) with Monte-Carlo simulations, the waves scattered by particles within a plate geometry (as shown in Figure 2) have to be simulated tens of thousands of times, with each simulation having hundreds of particles, before the standard error of the mean converges [29]. For each simulation we calculate exactly how the incident wave u inc (x) = e ikx scatters from all the particles. For these reasons, careful considerations are needed to determine how to perform the Monte-Carlo simulations.
In order, we explain how we created each particle configuration, how we determined the plate width and height, including consideration of convergence, and finally, how we analyzed the data.
Sequential addition. To place the particles we use the strategy of Sequential addition as described in Chapter 8, Section 2 of [43]. In essence, we place one particle at a time according to a random uniform distribution. If the particle overlaps with another particle it is rejected. The process is repeated until we obtain a desired particle volume fraction ϕ.
The plate width. Choosing an appropriate width W for the plate R was based on two factors: • The result from the theory shown by (12) predicts the plate needs to have width W > 2a+2b 12 for ⟨u(x)⟩ to be exactly equal to a sum of effective waves. The minimum value for b 12 is a 12 , which we found to be accurate enough for our tolerances. Using b 12 = a 12 > 2a implies that we need a plate with W > 6a = 7.2, as we used a = 1.2 for all numerical experiments.
• The plate width W can not be too wide, otherwise the average wave ⟨u(x)⟩ will be completely attenuated, which is a computational waste. Also, in the region where the wave is completely attenuated it is impossible to estimate the k p by fitting the formula (12). Materials, and frequencies, that lead to (12) needing more than one effective wavenumber k 1 to accurately approximate ⟨u(x)⟩ are highly attenuating materials. See Figure 7 for an illustration of the region where we fit the formula (12).
The plate height. If the plate filled with particles, as shown in Figure 2, was infinite in height, and the particles were excited by a plane-wave, then the average wave ⟨u(x)⟩ would be exactly a sum of plane-waves given by (12). See [15] for details. In practice, it is not of course possible to exactly simulate the wave scattered from an infinite plate filled with one specific arrangement of particles. The approximation often used is to have a cell filled with a random set of particles, and then to use periodic tilling of this cell [6]. To avoid the artifacts produced by periodic tilling we perform a convergence study to determine at what height a plate filled with particles behaves approximately like a plate of infinite size.
Let H be the height of the plate, which is illustrated in Figure 2. To determine the size needed for H, we first choose a large value H = 600, with the plate width W = 20, and then filled this plate with one configuration of particles Λ, according to the sequential addition method. We then calculate the total wave at to create a function U 600 (x). To determine the influence of H, we then reduce the height. For example, we use H = 590 and remove from Λ any particles that are now outside of the box with the reduced height. We then take the updated configuration of particles Λ and recalculate the total wave in the same region to create the function U 590 (x). For a range of heights we compute the relative error By calculating theses errors for a range of heights we can plot the error against the height as shown in Figure 11 below. For sound-soft particles, see Table 1 for details, and frequency ka = 0.3, we see that the errors have converged. This means that U 600 (x) is approximately the same as U ∞ (x), and we can estimate that the height of the plate has to be approximately H = 400 for the scattered waves to have an error of less than 1% in comparison to U ∞ (x).
For sound-hard particles, see Table 1 for details, and ka = 0.3, the errors have not converged as can be seen from Figure 11, even for very larger heights. This means it is unclear what is the error relative to an infinite plate U ∞ . For this reason, we do not focus on this case, and only performed one simulation with H = 400. Our hypothesis is that the plane-wave when scattered from the corner of the plate leads to a transmitted wave that travels down inside the plate relatively unobstructed, as these types of particles are weak scatterers.
When fitting for the effective wavenumbers, we need to consider how the errors due to truncating the height of the plate may affect the fitting and the acceptable fit errors.
where ϵ(x) is a small error that falls inside the standard error of the mean. For the cases where only one effective wavenumber was predicted, the average wave ⟨u(x, 0; Λ)⟩ is fitted well by (1), using non-linear optimization libraries in Julia [32]. For every wavenumber k p there are potentially two waves: one travelling towards the positive x -direction and another travelling in the negative x -direction. For the cases where the average wave is completely attenuated when it reaches the edge of the plate, x = 20, we should have A − p = 0, suppressing the wave travelling on the negative x-direction. From the Monte-Carlo data for cases with more than one effective wavenumber, we find that ⟨u(x, 0; Λ)⟩ is below the standard error of the mean for x ≥ 15. For this reason, we can take A − p = 0 and fit for only for A + p and k p . That is we fit to the Monte-Carlo data functions of the form: with Im[k p ] > 0. In the case presented in Figure 7a the transmitted wave attenuates quicker, therefore we focus the fitting on a narrower range, 4 ≤ x ≤ 10. Another reason to use (60), is that fitting for A + p , A − p , and k p in the case of multiple effective wavenumbers can become ill-posed. This is related to how inverting a Laplace transform is ill-posed, see Section 2.1 of [17]. That is, the more terms included in the sum shown in (12) the more ill-posed is the problem of recovering the k p , A + p , and A − p from data of the average wave ⟨u(x, 0; Λ)⟩.
For the cases with more than one effective wavenumber, non-linear optimization [32] is unstable. Therefore, we implement a method, the Projection method, to both fit the formula (60) and estimate the sensitivity of the parameters k p to the data.
Projection Method. To fit Monte-Carlo data for the cases with more than one effective wavenumbers, we design an algorithm that sweeps over all possible values of the wavenumbers k p and for each case performs a linear fit, by using least-squares, to predict the amplitudes A + p in (60). For that, we consider a mesh in the complex plane C ⊂ C from which we sample values of k p . A sketch of the algorithm for the case P = 3 is given below in Algorithm 1.
The algorithm calculates the minimum error ε(k 1 ) for every possible k 1 , and it is this error which is shown in Figures 8, 9b and 9d. We perform a further analysis to establish if the optimal wavenumbers are located at the local minima shown, and whether the fitted curves when using (60) are within the standard error of the mean of the Monte-Carlo simulation.

Algorithm 1: Algorithm for the Projection Method
Data: ⟨u(x, 0; Λ)⟩ for all values of k 1 ∈ C do for all values of k 2 , k 3 ∈ C do if k 1 ̸ = k 2 ̸ = k 3 then Determine the values of A + p that best fit ⟨u(x, 0; Λ)⟩ using least-squares ; Using (60), compute the error ε(k 1 , k 2 , k 3 ) = ∥h(x) − ⟨u(x, 0; Λ)⟩∥ ; Store ε(k 1 , k 2 , k 3 ); end end Store ε(k 1 ) = min k 2 ,k 3 ε(k 1 , k 2 , k 3 ) end Result: ε(k 1 ) The blue region in Figure 8, contains all wavenumbers for which the fitted curves are within the standard error of the mean. There are clearly two distinct disconnected regions for possible values of k p , which shows that the transmitted wave is composed of a sum of two effective waves.
Computational time. Our study utilizes high-fidelity Monte-Carlo simulations, whose computational cost was significant. This limits the number of cases we could investigate with Monte-Carlo. For example, we needed to execute 40,000 simulations for scenarios depicted in Figures 7, 8 and 9. Each simulation involve hundreds of particles, and are computationally demanding especially at higher frequencies denoted by ka. To manage this, we employed parallel processing across multiple processors. For instance, simulating sound-soft particles at a frequency of ka = 0.36 requires about 72 hours for completion on an 11th Gen. Intel Core i7 with 8 cores. So without parallel processing, the runtime would increase by at least eight-fold.

B The ensemble average transmission
Below we reduce the terms in (31) to reach equation (36). For this section we use the notation and assumptions introduced in Section 4.
In many calculations throughout this section, for any function f which depends on particle configuration, we use that where we used, in order, the definitions (26) and (28) and (32). If the f only depends on r 1 , then we further have that ⟨f⟩(r 1 ) = f becausê R J−1 1 p(r 2 , . . . , r J |r 1 )dr 2 · · · dr J = 1, which is true for any joint probability density.

B.1 The transmitted internal field
We start by calculating the simplest term. Using (24), the definitions of the ensemble average (27) and (26), then (61) leads to where the ball B(r; a) is defined by (39), and we use the assumption (35).

B.2 The transmitted incident field
The next simplest computation is the ensemble average of the incident wave term in (29). To do this we will demonstrate the following equalities: where ϕ is the particle volume fraction defined by (33). First, we use the ensemble average (26), then take u inc (r) outside of the integrals, as it does not depend on the particle positions, to reach ⟨u inc (r)χ R\P (r)⟩ = u inc (r)⟨χ R\P (r)⟩.
To calculate the ensemble average on the right we use where we used the definition that integrating a probability density function p over all its variables gives 1, and which holds because if any two particles overlap, we have that p(r 1 , r 2 , . . . , r J ) = 0. Therefore, if χ P j (r) = 1, then χ P ℓ (r) = 0 for ℓ ̸ = j. Next, we use that particles are indistinguishable, except for their positions, so after ensemble averaging ⟨χ P j (r)⟩ = ⟨χ P 1 (r)⟩ for every j, which together with (61) leads to ⟨χ P (r)⟩ = J j=1 ⟨χ P j (r)⟩ = J⟨χ P 1 (r)⟩ = nˆR 1 χ P 1 (r)dr 1 = nˆB (r,a) where we used (39) and (35).
To simplify the above we first use the Quasi-Crystalline Approximation (QCA) (72) ⟨f n ⟩(r 1 , r 2 ) ≈ ⟨f n ⟩(r 1 ), which is needed to deduce effective wavenumbers [15,23,26]. Before we show how to simplify (71) for a general pair-correlation g, which we shown in Appendix B.3.1, we first deduce the results for the simplest pair correlation called Hole-Correction. It is far easier to understand this case first. The Hole-Correction approximation is the result of taking b 12 = a 12 in the general pair-correlation (11).
The above implies that the region of integration for r 1 is just R 1 \ B(r; a 12 − a). We now further split this region of integration into two disjoint regions: the first is R 1 \ B(r; a 12 + a) and the second is B(r; a 12 + a) \ B(r; a 12 − a).
For the first region r 1 ∈ R 1 \ B(r; a 12 + a) implies that |r 1 − r| ≥ a 12 + a, which together with r 2 ∈ B(r; a) leads to |r 1 − r 2 | ≥ |r 1 − r| − |r − r 2 | > a 12 , due to the triangle inequality. In other words, the region of integration for r 2 becomes r 2 ∈ B(r; a) \ B(r 1 ; a 12 ) = B(r; a).
For the second region r 1 ∈ B(r; a 12 + a) \ B(r; a 12 − a), which implies that a < a 12 − a ≤ |r 1 − r| ≤ a 12 + a.
The above guarantees that the two spheres B(r; a) and B(r 1 ; a 12 ) will intersect. Let V be this region of intersection, then the region of integration of r 2 becomes B(r; a)\B(r 1 ; a 12 ) = B(r; a) \ V. This is useful as V is formed of two spherical caps whose volume is easy to calculate 4 . Using the split of these two regions for r 1 we obtain J(J − 1)⟨u 1 sc (r)χ P 2 (r)⟩ = n 2ˆR where V cap (d) is the volume of V and d = |r − r 1 |. By using the formulas for spherical caps 5 we can calculate that V cap (d) = π 12d (a + a 12 − d) 2 (d 2 + 2(a + a 12 )d − 3(a − a 12 ) 2 ).

B.3.1 An isotropic pair correlation
In the previous section we choose a simple pair-correlation to simplify the integral (73). Here we show how to reduce this integral when assuming a more general form for the isotropic pair correlation give by (11).
In the pair-correlation (11) we assume there is a value b 12 for which g(r) = 1 when r ≥ b 12 . This is an approximation, but it is essential for the results in this section. this distance b 12 , which is also called the correlation length, dictates at what distance inside the material the incident wave will be extinct.
Following closely the steps that led to (74), we now split the integral over r 1 into two regions R 1 \ B(r; b 12 + a) and B(r; b 12 + a) \ B(r; a 12 − a), which leads to J(J − 1)⟨u 1 sc (r)χ P 2 (r)⟩ = n 2ˆR where we use (35) to guarantee that the ball B(r; b 12 + a) is completely contained within the region R 1 . Without the condition (35) it does not seem possible to show that the incident wave becomes extinct, so we hypothesise that this is a necessary condition, as well as sufficient.
The integral on the right of the first line of (75) was already resolved in the previous section, except now we replace a 12 with b 12 . For the integrals on the second line of (75), we use the change of variables from r 2 to r 21 = r 2 − r 1 and r 1 to x 1 = r − r 1 to obtain J(J − 1)⟨u 1 sc (r)χ P 2 (r)⟩ = nϕˆR 1 \B(r;b 12 +a) ⟨u 1 sc (r)⟩(r 1 )dr 1 + n 2ˆB and r 21 = |r 21 |. This concludes the calculations in this section.