Multiple Waves Propagate in Random Particulate Materials

We show that in general there is not only one effective wave, but there is a series of effective waves in an ensemble averaged particulate material. By particulate material we mean a medium filled with randomly placed particles and consider waves governed by the scalar wave equation. Although most of these waves decay rapidly, they make a significant contribution to reflection and transmission beyond the low-frequency regime. In two spatial dimensions, we develop an efficient method to calculate all these waves, which gives highly accurate predictions when compared to a finite-difference scheme. This method is, to the authors knowledge, the first of its kind to give accurate predictions across a broad frequency range and particle volume fraction.

1. Introduction. Materials comprising small particles, inclusions or defects, randomly distributed inside an otherwise uniform host medium, are ubiquitous. Commonly occurring examples include composites, emulsions, dense suspensions, complex gases, polymers, and foods. Understanding how electromagnetic, elastic, or acoustic waves propagate through such media is crucial to characterize these materials and also to design new materials that can control wave propagation. For example, we may wish to use wave techniques to determine statistical information about the material, e.g., volume fraction of particles, particle radius distribution, etc.
The exact positions of all particles are usually unknown. The common approach to deal with this, which we adopt here, is to ensemble average over such unknowns. In certain scenarios, such as light scattering [38], it is easier to measure the average intensity of the wave, but these methods often need the ensemble-averaged field as a first step [14,54,53].
1.1. Historical perspective. The seminal work in this field is Foldy's 1945 paper [14], which introduced the Foldy closure approximation in order to deduce a single "effective wavenumber" k \ast in the form k \ast = k 0 -\phi g, where \phi is the volume fraction of particles and g is the scattering coefficient associated with a single particle. Foldy introduced the notion of ensemble averaging the field, but the expression deduced for k \ast was restricted to dilute dispersions and isotropic scattering. Lax improved on this by incorporating a higher-order closure approximation [25,26], now known as the quasi-crystalline approximation (QCA), and by including pair-correlation functions, which represent particle distributions. Both QCA and pair-correlations have now been extensively used in multiple scattering theory. The most commonly used pair-correlation is hole-correction [13]. Both QCA and hole-correction are examples of statistical closure approximations [2,3], which are techniques widely used in statistical physics. For multiple scattering, the accuracy of these approximations has been supported by theoretical [33,34], numerical [9], and experimental [59,67] evidence. These approximations also make no explicit assumptions on the frequency range, material properties, or particle volume fraction. We note, however, that to the best of our knowledge, there are no rigorous bounds for the error introduced by these approximations. For a brief discussion on these approximations see [21].
For an overview of the literature on multiple scattering in particulate materials, making use of closure approximations, see the books [55,31,38]. We now briefly summarize how calculating effective wavenumbers has evolved since the early work of Foldy and Lax. Over the last 60 years, corrections to the dilute limit have been sought, mainly by expanding in the volume fraction \phi and then attempting to determine the O(\phi 2 ) contribution to k \ast . Twersky [56] obtained an expression for this contribution as a function of f (\pi -2\theta inc ) and f (0), where \theta inc is the angle of incidence of an exciting plane wave (see Figure 1), and f (\theta ) is the far field scattering pattern from one particle [28]. The dependence on \theta inc implies that k \ast depends on the angle of incidence, which is counterintuitive. Waterman and Truell [63] obtained the same expression as Twersky but with \theta inc = 0. However, [63] used a``slab pair-correlation function"" that (theoretically) limits the validity of their approach to dilute dispersions (small \phi ); see [27] for comparisons with experiments, and see [5] for a discussion in two dimensions. Extensions that incorporate the hole-correction pair-correlation function were described by Fikioris and Waterman [13]. The Waterman and Truell expressions for three-dimensional elasticity are written down in [68,45]. Work in two-dimensional elasticity using QCA was reported by [8]. Lloyd and Berry [30] calculated the O(\phi 2 ) contribution by including both QCA and hole-correction for the scalar wave equation, although the language used stemmed from nuclear physics. More recently, [28,29] rederived the Lloyd and Berry formula for the effective wavenumber without appealing to the so-called extinction theorem used in many previous papers, such as [60], and without recourse to``resumming series."" The work was then extended in order to calculate effective reflection and transmission in [32]. Gower et al. [21] subsequently extended this result to model multispecies materials, i.e., to account for polydisperse distributions.
Other related work on effective wavenumbers and attenuation include comparing the properties of single realizations to those of effective waves [49,6,7,40], and effective waves in polycrystals [50,64] such as steel and ceramics. The polycrystal papers use a similar framework to waves in particulate materials, except they assume weak scattering which excludes multiple scattering.
1.2. Overview of this paper. A common assumption used across the field of random particulate materials, including those mentioned above, is to assume that there exists a single, unique, complex effective wavenumber k \ast that characterizes the material. For example, for an incident wave e ikx - i\omega t , of fixed frequency \omega , exciting a half-space (see Figure 1) filled with particles, the tacit assumption is that the Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license ensemble-averaged wave \langle u(x)\rangle travelling inside the particulate material takes the form See [32] for a brief derivation. This assumption has been widely used in acoustics [28,29,31,12], elasticity [61,42,43,46,10] (including thermo-viscous effects), electromagnetism [58,59,52], and even quantum waves [51]. For example, it is a key step in deducing radiative transfer equations from first principles [35,37].
In this work, we show, however, that there does not exist a single, unique effective wavenumber. Instead an infinite number of effective wavenumbers k 1 , k 2 , k 3 , . . . exist, so that the average field inside the particulate material takes the form In many scenarios, the majority of these waves are highly attenuating, i.e., k p has a large imaginary part for p > 1. In these cases, the least attenuating wavenumber k 1 will dominate the transmitted field in \langle u(x)\rangle , and k 1 will often be given by classical multiple scattering theory, as discussed in subsection 1.1. However, these other effective waves can still have a significant contribution to the reflected (backscattered) wave from a random particulate material, especially at higher frequencies and beyond the low volume fraction limit. Furthermore, there are scenarios where there are at least two effective wavenumbers, say k 1 and k 2 , with the same order of attenuation. In these cases using only one effective wavenumber, k 1 or k 2 , is insufficient to accurately calculate \langle u(x)\rangle , even for x far away from the interface between the homogeneous and particular materials.
We examine the simplest case that exhibits these multiple waves---two spatial dimensions (x, y) for the scalar wave equation---and consider particles placed in the half-space x > 0, which reflects incoming waves. We not only demonstrate that there are multiple effective wavenumbers, but we also use them to develop a highly accurate method to calculate \langle u(x)\rangle and the reflection coefficient. This method agrees extremely well with numerical solutions, calculated using a finite-difference method, but is more efficient. We provide software [17] that implements the methods presented and reproduce the results of this paper.
In a separate paper [18], we develop a proof that (1.2) is the analytical solution for the ensemble-averaged wave. However, the proof in [18] is not constructive, in contrast to the work presented here, where we present a method that determines all effective waves (1.2).
We begin by deducing the governing equation for ensemble-averaged waves in section 2. In section 3 we then show that multiple effective wavenumbers exist. To calculate these effective wavenumbers, we need to match them to the field near the interface at x = 0, which leads us to develop a discrete solution in section 4. The discrete solution also serves as the basis for a numerical method, which we use later as a benchmark. In section 5 we develop the matching method, which incorporates all of the effective waves. In section 6 we summarize the fields and reflection coefficients calculated by the matching method, the numerical method, and existing methods used in the literature. We subsequently compare their results in section 7. In section 8 we summarize the main results of the paper and discuss future work.
2. Ensemble-averaged multiple scattering. Consider a region \scrR filled with N particles or inclusions that are uniformly distributed. The field u is governed by Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license where k and k o are the real wavenumbers of the background and inclusion materials, respectively. For the sake of simplicity, we assume all particles are the same, except for their position and rotation about their center. For a distribution of particles, or multi-species, see [21].
In two dimensions, any incident wave 1 v j and scattered wave u j can be written in the form with (r j , \theta j ) being the polar coordinates of xx j , where x j = (x j , y j ) is a vector pointing to the center of the jth particle, from some suitable origin, and x is any vector; see Figure 1. The J n and H n are, respectively, Bessel and Hankel functions of the first kind. The representation (2.4) is valid when r j is large enough for (r j , \theta j ) to be outside of the jth particle for all \theta j . For example, in Figure 1 this distance would be r j > a o . The T-matrix is a linear operator, in the form of an infinite matrix, such that where we recall that N is the number of particles. The angle \tau j gives the particles' rotation about their center x j . Allowing particles to have different rotations, and assuming all \tau j \in [0, 2\pi ] to be equally likely, will lead to ensemble-averaged equations that are equivalent to the equations for circular particles [39]. This matrix T exists when scattering is a linear operation (elastic scattering), and can accommodate particles with a large variety of shapes and properties [15,16,62]; it is especially useful for multiple scattering [61,39,19].
The T-matrix also accounts for the particles' boundary conditions. For instance, if u represents pressure, \rho and c are the background density and wave speed, and the particles are circular with density \rho o , sound-speed c o , and radius a o , then continuity of pressure and displacement across the particles' boundary [28, section IV A] yields , In this case the T-matrix is independent of the rotation \tau j . In this paper we consider the incident plane wave (2.7) u inc (x, y) = e ik(x cos \theta inc +y sin \theta inc ) for \theta inc \in ( - \pi /2, \pi /2), which excites N particles, resulting in scattered waves of the form (2.4). See Figure 1 for an illustration. The total wave u, measured outside of all particles at x = (x, y), is the sum of all scattered waves plus the incident wave: To reach an equation for the coefficients u nj we write the total wave field incident on the jth particle v j (2.3) as a combination of the incident wave plus the waves scattered by the other particles: v j (r j , \theta j ) = u inc (x, y) + \sum i\not =j u i (r i , \theta i ). By then applying the Jacobi--Anger expansion to u inc (x, y), using Graf's addition theorem [31,1], multiplying both sides by T qn , summing over n, and then using (2.5), we obtain (2.9) u qj = u inc (x j , y j ) \infty \sum n= - \infty T qn (\tau j )e in(\pi /2 - \theta inc ) for all integers q and j = 1, 2, . . . , N , where (2.10) F n (X) = ( - 1) n e inΘ H n (R), and (R, \Theta ) are the polar coordinates of X.
2.1. Ensemble averaging. In practice the exact position of the particles is unknown, so rather than determine the scattering from an exact configuration of particles, we ensemble average the field u over all possible particle rotations and positions in \scrR . Sensing devices naturally perform ensemble averaging due to their size or from time averaging [36]. See [14,45,21] for an introduction to ensemble averaging of multiple scattering. For simplicity, we assume the particle positions are independent of particle rotations, so that the probability of the particles being centered at x 1 , x 2 , . . . , x N is given by the probability density function p(x 1 , x 2 , . . . , x N ). Hence, it follows that where each integral is taken over \scrR . Further, we have is the conditional probability density of having particles centered at x 1 , . . . , ✚ ✚ x j , . . . , x N (not including x j ), given that the jth particle is centered at x j . Given some function F (x 1 , . . . , x N ), we denote its ensemble average (over particle positions) by If we fix the location of the jth particle, x j , and average over the positions of the other particles, we obtain a conditional average of F given by We assume that one particle is equally likely to be centered anywhere in \scrR , when the position of the other particles is unknown: where we define the number density n = N/| \scrR | and the area of \scrR as | \scrR | .
Using the above we can express \langle u(x, y)\rangle , for (x, y) outside of the region \scrR , by taking the ensemble average of both sides of (2.8) to obtain where we assumed that all particles are identical (apart from their position and rotation). We also used equations (2.12), (2.15) and averaged both sides over particle rotations. Using the scattered field (2.4), we then reach The simplest scenario is the limit when the particles occupy the half-space x 1 > 0 [28], that is, \scrR = \{ (x, y)| x > 0\} . We focus on this case, although the method we present can be adapted to any region \scrR . In the limit of \scrR tending to a half-space, we let N \rightar \infty while n remains fixed. Due to the symmetry between the incident wave (2.7) and the half-space x 1 > 0, the field \langle u n1 \rangle x1 has a translational symmetry along y 1 , which allows us to write [21] (2.19) \langle u n1 \rangle x1 = \scrA n (kx 1 )e iky1 sin \theta inc . Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license For step-by-step details on deriving a governing equation for \scrA n (kx 1 ), see [28,10,21]. Here we only give an overview. First multiply by p(x 2 , . . . , x N | x 1 ) on both sides of (2.9), set j = 1, ensemble average over all particle rotations 2 and particle positions in x > 0, and then use the statistical assumptions hole-correction 3 and the quasicrystalline approximation to reach the system: where T m \delta mq = (2\pi ) - 1 \int 2\pi 0 T mq (\tau )d\tau , \delta mq = 1 if m = q and 0 otherwise, and a 12 is the minimum allowed distance between the center of any two particles. That is, a 12 is at least twice the radius for circular particles. For the case shown in Figure 1 we could choose a 12 = 2a o . This minimum distance a 12 guarantees that particles do not overlap.
When the T-matrix T is known, we can determine the field \scrA n from the system (2.20). The aim of this paper is to efficiently solve for \scrA n and in the process reveal that \scrA n is composed of a series of effective waves.
For the rest of the paper we employ the nondimensional variables where R o is the particles' nondimensional maximum radius (in Figure 1, R o = ka o ), \gamma \geq 2 is a chosen closeness constant, with \gamma = 2 implying that particles can touch, and \phi is the particle volume fraction. 4 Using nondimensional parameters helps to formulate robust numerical methods and to explore the parameter space.
3. Effective waves. An elegant way to approximate \scrA n is to assume it is a plane wave of the form [31] (3.1) \scrA n (X) = i n e - in\varphi A n e iXK cos \varphi for X > \= X, where K is the nondimensional effective wavenumber (kK is the dimensional effective wavenumber), with Im K \geq 0 to be physically reasonable, the factor i n e - in\varphi is for later convenience, and \= X is a length-scale we will determine later. We also restrict the complex angle \varphi by imposing that - \pi /2 < Re \varphi < \pi /2 and using which is due to the translational symmetry of equation (2.20) in y 1 ; see [21, equation (4.4)]. This relation is often called Snell's law.
As the material has been homogenized, it is tempting to make assumptions that are valid for homogeneous materials, such as assuming that only one plane wave (3.1) A. L. GOWER, W. J. PARNELL, AND I. D. ABRAHAMS is transmitted into the material. When the particles are very small in comparison to the wavelength, this is asymptotically correct [44], but in all other regimes this is not valid, especially close to the edge \= X = 0, as we show below. By substituting the ansatz (3.1) into (2.20), using (2.21) (see section SM1 in the supplementary material for details), and by restricting and (3.4) is often called the extinction theorem, though we will refer to it as the extinction equation. Using (3.3) we can calculate K by solving then the standard approach to calculate A n is to use (3.3) 1 and (3.4) and take \= X = 0, which avoids the need to know \scrA n or to calculate g( \= X). It is commonly assumed that there is only one viable K when fixing all the material parameters, including the incident wavenumber k. However, in general (3.7) admits many solutions, which we denote as K = K 1 , K 2 , . . . ; see Figure 2 for some examples. We order these wavenumbers so that Im K p increases with p. There is no reason why these wavenumbers are not physically viable. Therefore we write \scrA n as a sum of effective waves: where there are an infinite number of these effective wavenumbers [18], but to reach an approximate method we need only a finite number P . Technically, (3.8) is a solution to (2.20) for X > 0, that is, we could take \= X = 0. However, in this case, we found that close to X = 0 a very large number of effective waves P would be required to achieve an accurate solution. This is why we only use the sum of plane waves (3.8) for X > \= X > 0. One of these effective wavenumbers, in most cases the lowest attenuating K 1 , can be calculated using an asymptotic expansion for low \phi [28] and assuming it is a perturbation away from 1 (the background wavenumber).
Substituting  where the a p are determined from (3.9) 1 . However, only (3.10) remains to determine the vector α. As there is more than one effective wave, P > 1, (3.10) is not sufficient to determine α. This is because satisfying (3.9) and (3.10) only implies that the effective field (3.8) solves (2.20) for X 1 > \= X + \gamma R o . The missing information, needed to determine α, will come from solving (2.20) for 0 \leq X 1 < \= X + \gamma R o . We choose to do this by calculating a discrete solution for \scrA n within 0 \leq X 1 < \= X + \gamma R o , and then matching the \scrA n with the effective waves (3.8). The final result will be a (small) linear system (5.5).
Kristensson [23] deduced a similar one-dimensional integral equation for electromagnetism, and in [18] we showed that the analytic solution to (4.1) is a sum of effective plane waves.
We will use (4.1) to determine the effective waves (3.8) and to formulate a completely numerical solution to (2.20), which we use as a benchmark.
4.1. The discrete form. The simplest discrete solution of (4.1) is to use a regular spaced finite-difference method and a finite-section approximation. 5 A similar finite-difference solution was used in [22].
Let \scrA j n = \scrA n (X j ) for X j = jh and j = 0, . . . , J, with analogous notation for the other fields. We also define the vectors For implementation purposes, we consider all vectors to be column vectors. We also use the block matrix A with components A n1 = A n , that is, To discretize the integrals in (4.1), we use \int f (X)dX \approx \sum j f (X j )\sigma j , which in the simplest form is \sigma j = h for every j. Discretizing the integrals in (4.1), then substituting (4.2), and letting X j 1 = X j 2 = jh for j = 0, 1, . . . , J leads to where q = \lfloor R o \gamma /h\rfloor , \times \int X \ell +Ro\gamma The \sigma \ell j depend on \ell because the discrete domain of integration | j -\ell | \leq q changes with \ell , though the simplest choice would still be \sigma \ell j = h. If we did not include \scrE \ell nm and \scrR \ell nm , then the solution of (4.4) would represent the average wave in the layer 0 \leq X \leq X J . One method to calculate the solution for the whole half-space X \geq 0 is to extend X J until \scrA n (X J ) tends to zero. However, it is more computationally efficient to calculate \scrE \ell nm and \scrR \ell nm by approximating \scrA n (X) as a sum of plane waves, as shown below.
5. Matching the discrete form and effective waves. Here we formulate a system to solve for the unknown effective wave amplitudes \alpha p (3.11) and A (4.3). To do this, we substitute \scrA n (X 2 ) for the effective waves (3.8) in \scrE \ell nm and \scrR \ell nm , (4.6) and (4.7), and we calculate the integral g( \= X) in (3.10) by substituting \scrA n (X 2 ) for the discrete solution \scrA j n (4.2). Finally, to determine the \alpha p , and therefore the effective waves (3.8), we impose that (3.8) matches the discrete solution (4.2) in a thin layer near the boundary \= X. For an illustration see Figure 3. Imposing this match acts like a boundary condition for the effective waves. From here onwards we assume that \= X = X L .
X 0 X 1 . . . X L . . . X J - 1 X J \scrA 0 n \scrA 1 n \scrA L n i n \sum p e in\theta p A p n e iXKp cos \theta p Fig. 3. An illustration of the discrete solution A j n (4.2) (blue circles) and the effective waves (3.8) (black line). We restrict the coefficients A p n of the effective waves by imposing that the black line passes close to the A j n (i.e., satisfying the matching condition (5.10)) for X = X L , . . . , X J , where we choseX = X L . Increasing the number of effective waves will lead to a closer match between the discrete solution and effective waves. 5.1. Using the effective waves to calculate (4.6), (4.7). Substituting the effective waves (3.8) into (4.6), then integrating and using (A.2), we arrive at where we used X J -X \ell = X J - \ell \geq 0 for J \geq \ell , when substituting S n - m (X J -X \ell ) with (A.2). Employing (3.11), we write (5.1) in matrix form To calculate (4.7), we first discretize the integral then substitute the effective Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license where \sigma \ell j represents the discrete integral in the domain [X J , X \ell +q ]. Using (3.11), just as we did in (5.2), we can write the above in matrix form We can now rewrite the integral equation (4.4), using the above equations, in the compact form which is valid for all m. If α were known, then we could calculate the discrete solution A n from the above. However, the α also depends on the A n , as we show below.

5.2.
The effective waves in terms of the discrete form. The equations to determine the effective waves, so far, are (3.9) and (3.10). To calculate the integral in (3.10), we discretize and substitute (4.2), which leads to the discrete form of the extinction equation (3.10): where \cdot T denotes the transpose, we used (4.3), g( \= X) = G T A = \sum n G T n A n , w p = 2\phi \infty \sum n= - \infty e in\theta inc e - in\varphi p e i(Kp cos \varphi p - cos \theta inc )X L K p cos \varphi pcos \theta inc a p n , (5.7) (G n ) j = 2\phi e in\theta inc ( - i) n - 1 e - iX j cos \theta inc \sigma j , (5.8) and as the domain of the integral in (3.10) is only up to X L = \= X \leq X J , we set (G n ) j = 0 for j > L.
When using P effective wavenumbers, there are P unknowns \alpha 1 , . . . , \alpha P , with, so far, only one scalar equation (5.6) to determine them. To determine the \alpha p , we match the sum of effective waves (3.8) with the discrete form \scrA j n in the interval: X L < X < X J , as shown in Figure 3. To do this we could enforce (5.9) \scrA j n = i n \sum p e - in\varphi p e iX j Kp cos \varphi p a p n \alpha p = α T v j n for j = L, L + 1, . . . , J.
However, for n \not = 0 the coefficients \scrA j n and a p n can be very small, and the above would not enforce the extinction equation (5.6). So rather than use (5.9) for every n, it is more robust to minimize the difference: o cos \theta inc , Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license where the constraint enforces (5.6). For details on how to solve (5.10) see section SM2. The solution to the above is where w is the conjugate of w, and the block matrix L = [. . . , L - n , L 1 - n , . . . , L n , . . .], with Finally, substituting α (5.11) into (5.5), we reach an equation which we can solve for A: where E and R have components E m and R m , given by (5.2) and (5.4), respectively, while the components of the block matrices B and M are To summarize, the terms w, V, and L are given by (5.7), (5.13), and (5.12), Q mn is given by (4.5), and both A and b m are given by (4.2). The angle \theta inc is the angle of the incident plane wave (2.7), and R o is a nondimensional particle radius (2.21) which increases with the frequency. The block matrices G, B, A, E, R, L, and Z all have only one column. The elements of these columns are either column vectors (G m , B m , A m ) or matrices (E m , R m , L m , and Z m ).

The matching algorithm.
We can now understand how to truncate the effective wave series (3.8): assume the wavenumbers K p are ordered so that Im K p increases with p = 1, . . . , P . Then note that the larger Im (X J K p cos \varphi p ) is, the less the contribution this effective wave will make to the matching (5.10), w (5.6), \scrR \ell nm (5.3), and \scrE \ell nm (5.2). That is, we can choose P such that Im (X J K P cos \varphi P ) is large enough so that this wave will not affect the solution A.
To aid reproducibility, we explain how to solve (5.14), and determine A, by using an algorithm in section SM3.
6. The resulting methods. Here we summarize the matching method and other methods for solving (4.1). To differentiate between results for the different methods we use the superscripts M , D, and O. That is, we denote the field \scrA n (X) as \scrA M n (X) (matching method), \scrA D n (X) (discrete method), (6.1) \scrA O n (X) (one-effective-wave method). (6.2) For the matching method, we solve (5.14) to obtain (6.3) \scrA M n (X) = \Biggl\{ \scrA j n = (A n ) j , X = X j , i n \sum P p=1 e - in\varphi p e iXKp cos \varphi p a p n \alpha p , X > X J (matching method), Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license where the \alpha p are given from (5.11), and \varphi p , K p , and a p n are solutions to (3.2, 3.9, 3.11). For details on the matching method, see Algorithm SM3.1 in the supplementary material.
The one-effective-wave method is the typical method used in the literature. It consists in using only one effective wavenumber K 1 , that is, (3.1) with p = 1. This one wavenumber K 1 is often given explicitly in terms of either a low-volume-fraction or a low-frequency expansion. However, as we explore both moderate-frequency and moderate volume fractions, we will instead numerically solve for K 1 , the least attenuating wavenumber. To solve for K 1 and A 1 n we take \= X = 0 and numerically solve (3.9) and (3.10) for P = p = 1. The Snell angle \varphi 1 is determined from (3.2), with K = K 1 and \varphi = \varphi 1 . The result is (6.4) \scrA O n (X) = i n e - in\varphi 1 e iXK1 cos \varphi 1 A 1 n (one-effective-wave method).
From subsection 4.1, we can devise a purely numerical method, which requires a much larger meshed domain for X. The resulting field is This discrete method gives a solution for a material occupying the layer 0 < X < X J and Y \in R. If the layer is deep enough and the wave decays fast enough, then this discrete method will be the solution for an infinite half-space. Algorithm SM3.1 in the supplementary material can be used to calculate this discrete method by taking P = 1, J = L instead of step 7, as there is no matching region, and replacing steps 9--15 with: solve for A by using MA D = B instead of (5.14).
6.1. Reflection coefficient. The reflection coefficient R is the key information required for many measurement techniques. We can compare the different methods for calculating the average wave by comparing their resulting reflection coefficient, which is much simpler than comparing the resulting fields \scrA n (X).
Consider a particulate material occupying the region x > 0, and choose a point (x, y) to measure the reflection, with x < 0; then the ensemble average reflection coefficient R is such that (6.6) \langle u(x, y)\rangle = u inc (x, y) + Re ik( - x cos \theta inc+y sin \theta inc ) .
For the discrete method, we discretize (6.8), which leads to \sigma j \scrA j n e iX j cos \theta inc - in\theta inc (discrete method).
Alternatively, to obtain the reflection coefficient for one effective wave (6.4), we set J = 0 and P = 1 in (6.9) to reach which agrees with equations (41) and (42) from [32], when expanding for low volume fraction \phi .

Numerical experiments.
For simplicity, we consider circular particles (2.6) for all numerical experiments, in which case, the nondimensional radius (2.21) R o = a o k, where a o is the particle radius.
For the material properties we use a background material filled with particles which either strongly or weakly scatter the incident wave. These are given, respectively, by noting that \rho o \gg \rho leads to weaker scattering than \rho o \ll \rho . We will use a range of angles of incidence \theta inc , particle volume fractions \phi , and particle radiuses R o , which is equivalent to varying the incident wavenumbers k. Figure 4 shows several examples of \scrA M n from (6.3). As a comparison we have shown the one-effective-wave field \scrA O n (6.4) as well. To not clutter the figure, we have not shown the discrete field \scrA D n (6.5), which would lie exactly on top of \scrA M n . Figure 4 reveals how the discrete and effective wave parts of \scrA M n very closely overlap in the matching region X L \leq X \leq X J . This close overlap is not due to overfitting, as there are more than double the number of equations than unknowns.

Comparing the fields.
We now look closely at a specific case: particle volume fraction \phi = 20\% and nondimensional particle radius R o = 0.4 for the strong scatterers (7.1).  3) and the one-effective-wave field (6.4) for a material with circular particles, incident wave angle θ inc = 0, and properties (7.1). The nondimensional radius Ro = kao and volume fraction φ are shown on each graph. We used six effective wavenumbers (P = 6) for the bottom two graphs, and four effective wavenumbers (P = 4) for the top right graph. Note that the discrete and effective parts of the matching fields overlap in the match region. The one-effective-wave field in general loses accuracy close to the interface X = 0, which is why it gives inaccurate predictions for the reflection coefficient R O (6.11).
shows the effective wavenumbers used and how the greater the attenuation Im K p is, the lower the resulting amplitude | \alpha P | of the effective wave is, and therefore the less it contributes to the total transmitted wave. We also see in Figure 5(c) how increasing the number of effective waves (while fixing everything else) results in a smaller difference between the fields of the matching and discrete methods. This clearly confirms that the field \scrA n is composed of these multiple effective waves. Figure 6 shows how the matching method (6.3) and the discrete method (6.5) closely overlap with max X,n \| \scrA M n (X) -\scrA D n (X)\| = 4.5 \times 10 - 4 , which is similar to the matching error 4.7\times 10 - 5 given by the sum (5.10) 1 . The dotted and dashed curves in Figure 6 demonstrate how the matching method is only accurate when using the effective wavenumbers that satisfy (3.9). This agreement between the matching and discrete methods is not isolated to specific material properties and frequencies; we have yet to find a case where the two methods do not show excellent agreement. 6 Further, when increasing the number of effective wavenumbers P and lowering the tolerance tol in Algorithm SM3.1, the two methods converge to the same solution, as indicated by Figure 5(c). In this paper we will not explore this convergence in detail, but we will show that the two methods produce the same reflection coefficient for a large parameter range.  Figure (6). (a) shows the effective wavenumbers, with each marker corresponding to one wavenumber K P whose color darkens as the amplitude of its wave field α P increases. Clearly the larger the attenuation Im Kp, the lower the amplitude α P . (b) reveals how the amplitude α P decreases when the effective phase speed increases in magnitude.
(c) shows how the maximum error between the fields of the matching and discrete methods decreases when increasing the number of effective waves used by the matching method. Note: if we had not included the three lowest attenuating wavenumbers, the maximum error would be larger than 0.17. . This graph shows that the matching method (6.3) overlaps with the discrete method (6.5) (a purely numerical method). The effective wavenumbers used are shown in Figure 5, and the material properties are given by (7.1). The dashed and dotted curves also result from the matching method but use the wrong effective waves: the dotted curve wrong match A M 0 and A M 1 use the effective wavenumbers (3.9) multiplied by 1.2. The zero match fields zero all the effective wave amplitudes a p n = 0 and A M n (X) = 0 for X > 1.

Comparing reflection coefficients.
The reflection coefficient R is a simple way to compare the different methods in section 6. Many scattering experiments aim to estimate R [67,66]. The accuracy of estimating R is also directly related to the accuracy of calculating the transmitted waves.
In Figure 7 we compare the reflection coefficient for the discrete method R D (6.10), matching method R M (6.9), and two methods that use only one effective wavenumber (6.11): the one effective R O uses a numerical solution for K 1 (the wavenumber with the smallest imaginary part), while the low vol. frac R O uses a low-volumefraction expansion for the wavenumber [32].
In Figure 7(a) we compare the reflection coefficients for strong scatterers (7.1) when varying the particle radius R o (or likewise varying the wavenumber k) with a fixed volume fraction \phi = 20\%. We use at most 1600 points for the X mesh and fewer Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license than 100 points for the X mesh of the matching method, and aim for a tolerance of 10 - 5 for the fields.
We clearly see that R D and R M (6.9) overlap. For R o > 0.03 the maximum difference max Ro | R M -R D | < 0.0014. For R o < 0.03 we have not shown R D because the numerical truncation error became too large (compared to our tolerance). This occurs when the fields \scrA D (X) decay slowly, which occurs for small particles (or low frequency). However, for low frequency the one effective R O is asymptotically accurate [44], and we see that R M does converge to R O as R o \rightar 0. However, for larger R o the error of the one effective R O is as much as 20\%, while the low vol. frac. R O commits even larger errors. These larger errors are not unexpected, because the accuracy of the low-volume-fraction expansion depends on the type of scatterers and frequency [44], and can diverge in the limit R o \rightar 0 [20]. Figure 7(b) compares the reflection coefficients for weak scatterers (7.2). We use at most 2200 points for the X mesh and fewer than 100 points for the X mesh of the matching method, and aim for a tolerance of 10 - 5 for the fields.
Again, as before, we do not show R D for values of R o where the numerical truncation error become large (relative to our tolerance). For this case of weak scatterers we Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license see that the difference between the methods is less, though the reflection coefficient is also smaller with mean | R M | = 0.058. Still, the relative error of Im (R M - R O ) \approx 10\%. The imaginary part of the reflection coefficient, and where it changes sign, can be key for characterizing random microstructure [48]. The real part of the reflection coefficients is not shown, as the relative errors for the real part are even smaller.

Conclusions.
Our overriding message is that there is not one, but a series of waves, with different effective wavenumbers, that propagate (with attenuation) in an ensemble-averaged random particulate material. These waves must be included to accurately calculate reflection and transmission. Figure 2 shows examples of these effective wavenumbers.
Although there is an analytic proof [18] that there exist a series of effective waves, which solve (2.20), this current paper shows how to calculate these by using a matching method (6.3). In our numerical experiments in section 7, we show that the matching method converges to a numerical solution (the discrete method) for a broad range of wavenumbers k (or equivalently the nondimensional radius R o ), particle volume fractions, and two sets of material properties. For example, Figure 6 compares the average fields \scrA n (X), and Figure 7 the reflection coefficients R of the matching and discrete methods. The drawback of the discrete method (6.5) is that it is computationally intensive, especially for low wave attenuation, requiring a spatial mesh between 1600 and 2000 elements to reach the same tolerance as the matching method, which used only 100 elements.
For small incident wavenumbers k, the matching method converges to a result which assumes there exists only one effective wave for both strong and weak scatterers. Qualitatively, the fields \scrA n (X) from the one-effective-wave (6.4) and matching (6.3) methods agreed well when moving away from the material's interface; for example, see Figure 4. However, as the fields are not the same near the interface, the resulting reflection coefficients can significantly differ, as shown in Figure 7.
8.1. The next steps. Here we comment on a few directions for future work. One important limit, which we did not investigate here, is the low volume fraction limit: \phi \ll 1. In numerical experiments, not reported here, we found that the matching method converges to the one-effective-wave method in the limit for low \phi . It appears that as \phi decreases the Im K p , for p > 2, tends to +\infty , implying that the boundary layer \= X shrinks and makes all but K 1 insignificant. This limit deserves a detailed analytic investigation in a separate paper.
The consequences of this work have a direct impact on effective wave methods used for acoustic, elastic, electromagnetic, and even quantum wave scattering. That said, many of these fields use vector wave equations and require the average intensity. So one challenge is to translate the results of this paper to vector wave equations and the average intensity. Note that for electromagnetic waves, much of the groundwork for the average fields has already been done [23,24].
The radiative transfer equations are one outcome of properly deducing the averaged intensity for waves in particulate materials. For example, for electromagnetic waves, radiative transfer equations have been deduced under assumptions such as weak scattering, sparse particle volume fractions, and one effective wavenumber K 1 [37]. Within the confines of the assumptions used, radiative transfer methods (and modifications) are leading to accurate predictions of the reflected intensity [41,65,47]. We speculate that this work will eventually lead to accurate predictions for reflected intensity for a broad range of frequencies and particle properties. Downloaded 12/01/20 to 176.253. 25.198. Redistribution subject to CCBY license A. L. GOWER, W. J. PARNELL, AND I. D. ABRAHAMS Data and reproducibility. All results can be reproduced with the publicly available software [17], which has examples on how to calculate the effective wavenumbers and the matching method, as well as the finite-difference method that we present. i n e - in\theta inc e iX cos \theta inc , X \geq 0, ( - i) n e in\theta inc e - iX cos \theta inc , X < 0.
Then as X is bounded by | X| < R o \gamma , we can choose Y 1 such that X/Y 1 is below a prescribed tolerance.