A proof that multiple waves propagate in ensemble-averaged particulate materials

Effective medium theory aims to describe a complex inhomogeneous material in terms of a few important macroscopic parameters. To characterize wave propagation through an inhomogeneous material, the most crucial parameter is the effective wavenumber. For this reason, there are many published studies on how to calculate a single effective wavenumber. Here, we present a proof that there does not exist a unique effective wavenumber; instead, there are an infinite number of such (complex) wavenumbers. We show that in most parameter regimes only a small number of these effective wavenumbers make a significant contribution to the wave field. However, to accurately calculate the reflection and transmission coefficients, a large number of the (highly attenuating) effective waves is required. For clarity, we present results for scalar (acoustic) waves for a two-dimensional material filled (over a half-space) with randomly distributed circular cylindrical inclusions. We calculate the effective medium by ensemble averaging over all possible inhomogeneities. The proof is based on the application of the Wiener–Hopf technique and makes no assumption on the wavelength, particle boundary conditions/size or volume fraction. This technique provides a simple formula for the reflection coefficient, which can be explicitly evaluated for monopole scatterers. We compare results with an alternative numerical matching method.

2. The Matching Method should be briefly mentioned in the abstract as it forms the basis for the work presented. 3. Page 4. The spatial extent of the incident plane wave is not mentioned as an assumption that simplifies the analysis. A plane wave with spatial extent in y would encounter different boundary conditions along the boundary as the relative acoustic impedance changes, due to local changes in the particulate material. Also near and far field effects would need to be considered in the analysis. 4. Page 4, lines [47][48][49]. The statement about the agreement between the methods is not clearly referenced. It could be read as a conclusion of the present work. On behalf of the Editor, I am pleased to inform you that your Manuscript RSPA-2019-0344 entitled "A Proof that Multiple Waves Propagate in Ensemble-Averaged Particulate Materials" has been accepted for publication subject to minor revisions in Proceedings A. Please find the referees' comments below.
The reviewer(s) have recommended publication, but also suggest some minor revisions to your manuscript. Therefore, I invite you to respond to the reviewer(s)' comments and revise your manuscript. Please note that we have a strict upper limit of 28 pages for each paper. Please endeavour to incorporate any revisions while keeping the paper within journal limits. Please note that page charges are made on all papers longer than 20 pages. If you cannot pay these charges you must reduce your paper to 20 pages before submitting your revision. Your paper has been ESTIMATED to be 21 pages pages. We cannot proceed with typesetting your paper without your agreement to meet page charges in full should the paper exceed 20 pages when typeset. If you have any questions, please do get in touch.
It is a condition of publication that you submit the revised version of your manuscript within 7 days. If you do not think you will be able to meet this date please let me know in advance of the due date.
To revise your manuscript, log into https://mc.manuscriptcentral.com/prsa and enter your Author Centre, where you will find your manuscript title listed under "Manuscripts with Decisions." Under "Actions," click on "Create a Revision." Your manuscript number has been appended to denote a revision.
You will be unable to make your revisions on the originally submitted version of the manuscript. Instead, revise your manuscript and upload a new version through your Author Centre.
When submitting your revised manuscript, you will be able to respond to the comments made by the referee(s) and upload a file "Response to Referees" in "Section 6 -File Upload". You can use this to document any changes you make to the original manuscript. In order to expedite the processing of the revised manuscript, please be as specific as possible in your response to the referee(s).
IMPORTANT: Your original files are available to you when you upload your revised manuscript. Please delete any redundant files before completing the submission process.
In addition to addressing all of the reviewers' and editor's comments, your revised manuscript MUST contain the following sections before the reference list (for any heading that does not apply to your work, please include a comment to this effect): • Acknowledgements • Funding statement See https://royalsociety.org/journals/authors/author-guidelines/ for further details.
When uploading your revised files, please make sure that you include the following as we cannot proceed without these: 1) A text file of the manuscript (doc, txt, rtf or tex), including the references, tables (including captions) and figure captions. Please remove any tracked changes from the text before submission. PDF files are not an accepted format for the "Main Document".
2) A separate electronic file of each figure (tif, eps or print-quality pdf preferred). The format should be produced directly from original creation package, or original software format.
3) Electronic Supplementary Material (ESM): all supplementary materials accompanying an accepted article will be treated as in their final form. Note that the Royal Society will not edit or typeset supplementary material and it will be hosted as provided. Please ensure that the supplementary material includes the paper details where possible (authors, article title, journal name). Supplementary files will be published alongside the paper on the journal website and posted on the online figshare repository (https://figshare.com). The heading and legend provided for each supplementary file during the submission process will be used to create the figshare page, so please ensure these are accurate and informative so that your files can be found in searches. Files on figshare will be made available approximately one week before the accompanying article so that the supplementary material can be attributed a unique DOI.
Alternatively you may upload a zip folder containing all source files for your manuscript as described above with a PDF as your "Main Document". This should be the full paper as it appears when compiled from the individual files supplied in the zip folder.

Article Funder
Please ensure you fill in the Article Funder question on page 2 to ensure the correct data is collected for FundRef (http://www.crossref.org/fundref/).

Media summary
Please ensure you include a short non-technical summary (up to 100 words) of the key findings/importance of your paper. This will be used for to promote your work and marketing purposes (e.g. press releases). The summary should be prepared using the following guidelines: *Write simple English: this is intended for the general public. Please explain any essential technical terms in a short and simple manner. *Describe (a) the study (b) its key findings and (c) its implications. *State why this work is newsworthy, be concise and do not overstate (true 'breakthroughs' are a rarity). *Ensure that you include valid contact details for the lead author (institutional address, email address, telephone number).

Cover images
We welcome submissions of images for possible use on the cover of Proceedings A. Images should be square in dimension and please ensure that you obtain all relevant copyright permissions before submitting the image to us. If you would like to submit an image for consideration please send your image to proceedingsa@royalsociety.org Once again, thank you for submitting your manuscript to Proceedings A and I look forward to receiving your revision. If you have any questions at all, please do not hesitate to get in touch. Comments to the Author(s) The manuscript demonstrates that one may have more than one effective wavenumber when describing wave propagation in an inhomogeneous medium via the effective medium (ensemble average) approach. The result is established for a half-infinite 2D space via the Wiener-Holf technique. Overall, the paper is interesting and well written. However, I do have a few concerns/suggestions that I hope the authors can address in the revised manuscript.
The assumption that the free-space (vacuum) wavenumber has a small imaginary part seems critical to all of the derived results, e.g., see eq. (3.10). Physically, $k_0$ (wavenumber inside a particle) may and should contain a small imaginary part but this is not the case for the freespace. I understand the mathematical necessity to add a small imaginary part to $k$ but I wonder if one could obtain similar results for strickly real $k$. I suspect the answer is no and thus, the relevance of this manuscript to real physical systems is in question. It would be great if authors could elaborate on this point beyond "[The] dissipation will facilitate the use of the Wiener-Hopf technique, and after reaching the solution we can take $k$ to be real." Decision letter (RSPA-2019-0344.R1)

14-Aug-2019
Dear Mr Gower I am pleased to inform you that your manuscript entitled "A Proof that Multiple Waves Propagate in Ensemble-Averaged Particulate Materials" has been accepted in its final form for publication in Proceedings A.
Our Production Office will be in contact with you in due course. You can expect to receive a proof of your article soon. Please contact the office to let us know if you are likely to be away from email in the near future. If you do not notify us and comments are not received within 5 days of sending the proof, we may publish the paper as it stands. Note that if you have opted for open access then payment will be required before the article is published? payment instructions will follow shortly.
If you wish to opt for open access then please inform the editorial office (proceedingsa@royalsociety.org) as soon as possible.
Your article has been estimated as being 21 pages long. Our Production Office will inform you of the exact length at the proof stage.
Proceedings A levies charges for articles which exceed 20 printed pages. (based upon approximately 540 words or 2 figures per page). Articles exceeding this limit will incur page charges of ?150 per page or part page, plus VAT (where applicable).
Under the terms of our licence to publish you may post the author generated postprint (ie. your accepted version not the final typeset version) of your manuscript at any time and this can be made freely available. Postprints can be deposited on a personal or institutional website, or a recognised server/repository. Please note however, that the reporting of postprints is subject to a media embargo, and that the status the manuscript should be made clear. Upon publication of the definitive version on the publisher?s site, full details and a link should be added.
You can cite the article in advance of publication using its DOI. The DOI will take the form: For tips on promoting your accepted paper see our blog post: https://blogs.royalsociety.org/publishing/promoting-your-latest-paper-and-tracking-yourresults/ Effective medium theory aims to describe a complex inhomogeneous material in terms of a few important macroscopic parameters. To characterise wave propagation through an inhomogeneous material, the most crucial parameter is the effective wavenumber. For this reason, there are many published studies on how to calculate a single effective wavenumber. Here we present a proof that there does not exist a unique effective wavenumber; instead, there are an infinite number of such (complex) wavenumbers. We show that in most parameter regimes only a small number of these effective wavenumbers make a significant contribution to the wave field. However, to accurately calculate the reflection and transmission coefficients, a large number of the (highly attenuating) effective waves is required. For clarity, we present results for scalar (acoustic) waves for a two-dimensional material filled (over a half space) with randomly distributed circular cylindrical inclusions. We calculate the effective medium by ensemble averaging over all possible inhomogeneities. The proof is based on the application of the Wiener-Hopf technique and makes no assumption on the wavelength, particle boundary conditions/size, or volume fraction. This technique provides a simple formula for the reflection coefficient, which can be explicitly evaluated for monopole scatterers. We compare results with an alternative numerical matching method.

Introduction
Materials comprising particles or inclusions that are randomly distributed inside a uniform host medium occur frequently in the world around us. They occur as synthetically fabricated media and also in nature. Common examples include composites, emulsions, suspensions, complex gases, and polymers. Understanding how electromagnetic, elastic, or acoustic waves propagate through these materials is necessary in order to characterise the properties of these materials, and also to design new materials that can control wave propagation.
The wave scattered from a particulate material will be influenced by the positions and properties of all particles, which are usually unknown. However, this scattered field, averaged over space or over time, depends only on the average particle properties. Many measurement systems perform averaging over space, if the receivers or incident wavelength are large enough [1], or over time [2]. In most cases, this averaging process is the same as averaging over all possible particle configurations. Such systems are sometimes called ergodic [2,3]. In this paper, we focus on ensemble averaged waves, satisfying the scalar wave equation in two-dimensions, reflecting from, and propagating in, a half-space particulate material. In certain scenarios, such as light scattering, it is easier to measure the average intensity of the wave. However, even in these cases, the ensemble-averaged field is often needed as a first step [4,5].
One driving principle, often used in the literature, is that the ensemble-averaged wave itself satisfies a wave equation with a single effective wavenumber [6][7][8]. Reducing a inhomogeneous material, with many unknowns, down to one effective wavenumber is attractive as it greatly reduces the complexity of the problem. For this reason many papers have attempted to deduce this unique effective wavenumber from first principles in electromagnetism [3,9,10], acoustics [11][12][13][14][15] and elasticity [16,17]. See [18] for a short overview of the history of this topic, including typical statistical assumptions employed within the methods, such as hole-correction and the quasi-crystalline approximation, which we also adopt here.
The assumption that the ensemble averaged wave field satisfies a wave equation, with an effective wavenumber, has never been fully justified. Here we prove that there does not exist a unique effective wavenumber but instead there are an infinite number of them. Gower et al. [18] first showed that there exist many effective wavenumbers, and provided a technique, the Matching Method, to efficiently calculate the effective wave field. In the present paper and [18], we show that for some parameter regimes, at least two effective wavenumbers are needed to obtain accurate results, when compared with numerical simulations. We also provide examples of how a single effective wave approximation leads to inaccurate results for both transmission and reflection for a halfspace filled with particles, see Figure 1.
Although the Matching Method developed in [18] gave accurate results, when compared to numerical methods and known asymptotic limits, the limitations of the method were not immediately clear. Here however we illustrate that the Matching Method is robust, because combining many effective wavenumbers is not just a good approximation, it is an analytical solution to the integral equation governing the ensemble averaged wave field. We prove this by employing the Wiener-Hopf technique and then, for clarity, illustrate the solution for particles that scatterer only in their monopole mode. The Wiener-Hopf technique also gives a simple and elegant expression for the reflection coefficient.
The Wiener-Hopf technique is a powerful tool to solve a diverse range of wave scattering problems, see [19,Chapter 5. Wiener-Hopf Technique] and [20,21] for an introduction. It is especially useful for semi-infinite domains [22][23][24][25][26][27] and boundary value problems of mixed type. In this work, the Wiener-Hopf technique clearly reveals the form of the analytic solution, but to compute the solution would require an analytic factorisation of a matrix-function. To explicitly perform this factorisation is difficult [28][29][30][31]. Indeed this is often the hardest aspect of employing the Wiener-Hopf technique, although there exist approximate methods for this purpose [28,[32][33][34]. We do not focus in this article on these analytic factorisations, as there already exists a method to compute the required solution [18]. Instead, the present work acts as proof that the Matching Method [18] faithfully reproduces the form of the analytic solution.
Ensembled averaged particulate material Im sp, the more quickly the wave attenuates as it propagates into the half-space and the smaller the drawn vector for that wave above. The results shown here represent the effective wavenumbers for parameters (5.2), which are shown in Figure 3. Figure 1 shows the main setup and result of this paper: an incident plane wave excites the halfspace x > 0 filled with ensemble-averaged particles (the blue region), which generates a reflected wave and many effective transmitted waves. The sp are the transmitted wavevectors, and the smaller the length of the vector, the faster that effective wave attenuates as it propagates further into the material.
The paper begins by summarising the equations that govern ensemble averaged waves in two-dimensions in section 2. Following this, in section 3 we apply the Wiener-Hopf technique to the governing integral equation and deduce that the solution is a superposition of plane waves, each with a different effective wavenumber. A simple expression for the reflection coefficient is also derived. In section 4 we specialise the results for particles that scatter only in the monopole mode, which leads to a closed form analytic solution.
The dispersion relation (3.30), derived in section 3, admits an infinite number of solutions, the effective wavenumbers. In section 5, we deduce asymptotic forms for the effective wavenumbers in both a low and high frequency limit. In section 6 we compare numerical results for monopole scatterers, using the Wiener-Hopf technique, with classical methods that assume only one effective wavenumber [11,13], and the Matching Method introduced in [18]. In general, when comparing predicted reflection coefficients, the Wiener-Hopf and Matching Method agree well, whereas the classical single-effective-wavenumber method can disagree by anywhere up to 20%. These results are discussed in section 7 together with anticipated future steps.

Waves in ensemble averaged particles
Consider a region filled with particles or inclusions that are uniformly distributed. The field u is governed by the scalar wave equations: where k and ko are the real wavenumbers of the background and inclusion materials, respectively. We assume all particles are identical, except for their position and orientation, for simplicity. For a distribution of particles, or multi-species, see [15].
Our goal is to calculate the ensemble average field u(x, y) , that is, the field averaged over all possible particle positions and orientations. For clarity, and ease of exposition, we consider that the particles are equally likely to be located anywhere except that they cannot overlap (this is often called the hole correction assumption). We also assume the quasi-crystalline approximation; for details on this, and for further details on deducing the results in this section, see [11,15,18].
By splitting the total steady wave field u(x, y) into a sum of the incident wave u inc (x, y) and waves scattered by each particle, the jth scattered wave being u j (x, y), we can write: A simple and useful scenario to consider is when all particles are placed only within the halfspace 1 x > 0, which are then excited by a plane wave, with implicit time dependence e iωt , incident from a homogeneous region: where we restrict the incident angle − π 2 < θ inc < π 2 , as shown in Figure 1, and consider a slightly dissipative medium with Re k > 0 and Im k > 0. (2.5) This dissipation will facilitate the use of the Wiener-Hopf technique, and after reaching the solution we can take k to be real 2 .
To describe the particulate medium we employ the following notation: b = the minimum distance between particle centres, (2.6) n = number of particles per unit area, (2.7) Tn = the coefficients of the particle's T-matrix, Although the area fraction φ, normally called the volume fraction, is a combination of other parameters, it is useful because it is non-dimensional. If we let ao be the maximum distance from the particle's centre to its boundary, then we can set b = γao, where γ ≥ 2 so as to avoid two particles overlapping. The volume fraction that does not include the exclusion zone φ , as used in [18, equation (4.7)], is then φ = 4φ/γ 2 . The Tn are the coefficients of a diagonal T-matrix [35][36][37][38][39]. The T-matrix determines how the particle scatters waves, and so depends on the particle's shape and boundary conditions. A diagonal T-matrix can be used to represent either a radially symmetric particle, or particles averaged over their orientation, assuming the orientations have a random uniform distribution.
The results of ensemble averaging (2.3) from first principals are deduced in a number of references [15,18] and so deails of this procedure are omitted here for brevity. To represent the ensemble averaged scattered wave from a particle, whose centre is fixed at (x 1 , y 1 ), we use is on the outside of this particle, with (R, Θ) being the polar coordinates of (X, Y ), H n are Hankel functions of the first kind, and An is some field we want to determine 3 .
By choosing 4 x < −b, which is outside of the region filled with particles, then taking the ensemble average on both sides of (2.3) results in equation (6.7) of [18], given by: which is the incident wave plus an effective reflected wave with reflection coefficient: where we assumed particles are distributed according to a uniform distribution, and the kernel ψn is given by Later we show that, as expected, R is independent of x.
The system governing Am(x) is given by equation (4.7) of [18]: for all integers m. Kristensson [40, equation (15)] presents an equivalent integral equation for electromagnetism and particles in a slab. Our main aim is to reach an exact solution for An(x) by employing the Wiener-Hopf technique to (2.14). We show how this also leads to simple solutions for the reflection coefficient by using (2.11). We acknowledge the authors of [11], as they noticed that (2.14) is a Wiener-Hopf integral equation; but apparently did not follow the steps indicated in the following sections 5 .

Applying the Wiener-Hopf technique
Equation (2.14) is convolution integral equation with a difference kernel. This means applying a Fourier transforms can lead to elegant and simple solutions. To facilitate, we must analytically extend (2.14) for all x 1 ∈ R by defining for integers m, where if the An(x) were known for x > 0, then the Dn(x) would be given from the left hand-side. Note that the kernel ψn defined in (2.13) is already analytic in the domain R. 3 The factor e iβy 1 appears due to the translational invariance of u1(x1, y1) in y1, which is a result of the material being statistically homogeneous, see [15] for details. 4 We define the reflection coefficient only for x < −b, instead of x < −b/2, so that we can use ψn in the formula for R, which will in turn facilitate calculating R. 5 However, they were unable to solve it because, it seems, of a mistake in the integrand of equation (37)  The field D 0 (x) is not just an abstract construct, it is closely related to the reflected wave: by directly comparing (3.1) with the reflection coefficient (2.12), for x < −b, we find that To solve (3.1) we employ the Fourier transform and its inverse, which we define aŝ for any smooth function f . We then definê We can determine whereÂ + n andD − n are analytic by assuming 6 that for some (possibly small) positive constant c. This leads toÂ + n (s) being analytic for Im s > −c, whileD − n (s) is analytic for Im s < c. In other words, bothÂ + n (s) andD − n (s) are analytic in the overlapping strip |Im s| < c. (3.7) To apply the Wiener-Hopf technique we also need to specify the large s behaviour for botĥ A + n (s) andD + n (s). To achieve this, we assume, on physical grounds, that An(x) is bounded when x → 0 + , and Dn(x) is bounded when x → 0 − . Then, it can be shown [21,41] that An in whichψn(s) is well defined (i.e. analytic) for s in the strip: see appendix A for details. The right-hand side of (3.1) becomeŝ where for the last step we assumed Im (s + α) > 0, which is why we use the superscript + on (s + α) + . This assumption, together with (3.7) and (3.10), is satisfied if |Im s| < , where = min{c, (1 − | sin θ inc |)Im k, Im α}. If (3.12) is satisfied then we can combine (3.9),(3.11) and (A 6), to obtain the Fourier transform of (3.1) in matrix form: where, for reference, andψ n−m (s) is given by (A 6). In the above θ S and S are chosen to satisfy Later we identify S and θ S as the effective wavenumber and transmission angle. The above does not determine the sign of S for any given complex s. To fully determine S and θ S , we take sgn(Re s) = sgn(Re S) which together with (3.19) leads to where both S and θ S , when considered as functions 7 of s, contain branch-points at s = ±ik sin θ inc with finite branch-cut running between −ik sin θ inc and ik sin θ inc . However, Ψ (s) is an entire matrix function having only zeros in s and no branch-points; see the end of appendix A for details. Determining the roots of det Ψ (s) = 0 will be a key step in solving (3.13), and so the following identities will be useful

(a) Multiple waves solution
To solve (3.13), we use a matrix product factorisation [42] of the form: where Ψ − (s), and its inverse, are analytic in Im s < , and Ψ + (s), and its inverse, are analytic for Im s > − . See (3.12) for the definition of . For our purposes, it is enough to know that such a factorisation exists [42], as this will lead to a proof that A(x) is a sum of attenuating plane waves.
Multiplying both sides of (3.13) by [Ψ − (s)] −1 and by (s − α) − leads to where (s + α) + is analytic for Im s > −Im α, while (s − α) − is analytic for Im s < Im α. We need to rewrite the last term above as a sum of a function which is analytic in the upper half-plane (Im 7 With our choice of branch-cut, the simplest way to compute S, with most software packages, is to use S = sgn(Re s) s 2 + (k sin θinc) 2 with the default cut location for the square root function, i.e. along the real negative line. s > − ) and another analytic in the lower half-plane. This is achieved below: so that g − (s) does not have a pole at s = −α and is therefore analytic for Im s < .

Substituting (3.25) into (3.24) leads to
Because both sides are analytic in the strip |Im s| < , we can equate each side to E(s), some analytic function in the strip. Further, as the left-hand side (right-hand side) of (3.26) is analytic for Im s > ( Im s < − ), we can analytically continue E(s) for all s, i.e. E(s) is entire.
To determine E(s) we need to estimate its behaviour as |s| → ∞. From (3.8) we have that A + (s) = (|s| −1 ) as |s| → ∞ in the upper half-plane, and from (3.15 -3.17): for s in the strip (3.12). From this we know that the factors Ψ + (s) and Ψ − (s) must be O(|s|) as |s| → ∞, in their respective half-planes of analyticity [28]. So, the left hand-side of (3.26) behaves as O(|s| −1 ) as |s| → ∞ in Im s > − . We can therefore use Liouville's theorem to conclude that E(s) ≡ 0, which means the Wiener-Hopf equation (3.26) is formally equivalent tô Let C + (s) be the cofactor matrix of Ψ + (s), so that From the property (3.22) 1 we can write det Ψ (s) = f (s 2 ) for some function f . Then, for every root s = sp of det Ψ (s), with Im sp > 0, we have that −sp is also a root, and vice-versa. From here onwards we assume: det Ψ (sp) = det Ψ (−sp) = 0 with Im sp > 0 and p = 1, 2, · · · , ∞.  To use the residue theorem below, we need to calculate det Ψ + (s) for s close to the root −sp, in the form

Im s
Re s α −s p Figure 2. An illustration of the contour integral over C D , used to calculate (3.34) for x < 0, and the contour integral over C A , used to calculate (3.32) for x > 0. The −sp (the red points) are roots of (3.30), and also the poles of (3.28). The single blue point α is the only pole of (3.29).
Using the above, and that C + (S) is analytic for Im s > − , we can apply an inverse Fourier transform (3.3) 2 to both sides of (3.28) and using residue calculus we find (3.32) is, by Jordan's lemma, the same as a clockwise integral over the closed contour C A which surrounds the poles −s 1 , −s 2 , . . ., i.e. roots of (3.30), as shown by Figure 2. Note that the cofactor matrix C + (s) contains no poles and so does not contribute additional residual terms. The yellow striped region in Figure 2 is the domain where Ψ is analytic. On the other hand, for x < 0, the integral (3.32) is the same as an integral over the counter-clockwise closed contour within the region Im s > 0 (not shown in Figure 2). The integrand has no poles in this domain and hence evaluates to zero. Likewise, by applying an inverse Fourier transform to (3.29), we obtain: For x < 0 the above integral is the same as a counter-clockwise closed integral over C D which surrounds the pole s = α (recalling that Im α > 0), as shown in Figure 2. The result is just the residue at this pole. That is, the function Ψ − (s)g − (s) contains no other singularities within Im s > 0. On the other hand, for x > 0 the integral is the same as a closed clockwise integral around the region Im s < 0 which evaluates to zero, as there are no singularities in this region (not shown in Figure 2). Clearly (3.32) shows that A(x) is a sum of plane waves with different effective wavenumbers sp, each satisfying (3.30). In section 5 we discuss these roots in more detail, and in section 6, we see that usually only a few effective wavenumbers are required to obtain accurate results.

(b) Reflection coefficient
By substituting (3.34) in (3.2) leads to Alternatively, the reflection coefficient can be calculated from (2.12) by employing the form of A(x) from (3.32), which is the more common approach. To simplify, we use i n e −inθinc e iαX for X > 0, (3.36) which then implies that ψn( The above is shown in [43, equation (37)] and [11, equation (65)]. This result together with (3.32) substituted into (2.12) leads to the form where we used that Im sp > 0. The above agrees with [13, equation (39)] and 8 [18, equation (6.9)].

Monopole scatterers
For particles that scatter only in their monopole mode, i.e. the scattered waves are angularly symmetric about each particle, we can easily calculate the factorisation (3.23). This type of scattered wave tends to dominate in the long wavelength limit for scatterers with Dirichlet boundary conditions. In acoustics, these correspond to particles with low density or low sound speed.
Once we know the factorisation (3.23), we can then calculate the average scattering coefficient (3.32) and average reflection coefficient (3.35). We will compare both of these against predictions from other methods in section 6.
For monopole scatterers we use S 2 − k 2 = s 2 − α 2 and rewrite Ψ 00 (s) = (s 2 − α 2 )q(s), with q(s) = 1 + 2πn with N 0 (bS) given by (3.17). Then, because q(s) → 1 as |s| → ∞, we can factorise q(s) = q − (s)q + (s) using where the integral path for q + (s) (q − (s)) has to be in the strip where q(s) is analytic, with the path for q + (s) (q − (s)) passing below (above) z. We then have 9 that where (4.3) 3 holds if −s is below the integration path of (4.2) and s is above the integration path of (4.1). From (3.32) we see that we need only evaluate Ψ + 00 (s), and therefore q + (s), for s = s 1 , s 2 , . . . , sp where as p increases, the sp become more distant from the real line. Then for large z, by inspection of (3.17), we have that and therefore we can accurately approximate the integral (4.1) by truncating the integration domain for large z.

(b) Explicit solution for monopole scatterers.
For monopole scatterers An(x) = Dn(x) = 0 for |n| > 0. Using this in (3.14 -3.17) leads to all vectors and matrices having only one component, given by setting n = m = 0. In this case A (3.32) reduces to for x > 0, where we used (4.3), C + (s) = 1, B = iT 0 , and dΨ00 ds (−s) = − dΨ00 ds (s) for every s. Likewise for (3.35) we arrive at Alternatively, using (3.37), we can calculate the contribution of P effective waves to the reflection coefficient where the error |R P − R| then indicates how many effective waves are needed to accurately describe the field near the boundary x = 0. The blue points represent waves travelling forwards (i.e. deeper into the material), while the red represent waves travelling backwards. All these waves are excited in a reflection experiment. Two wavenumbers in particular stand out as having the lowest attenuation S 1 and S 2 , both inside the grey dashed circle. The graph on the right is a magnification of the region close to these two wavenumbers. Out of these two, most efforts in the literature have focused on calculating S 2 , as it often has the lowest imaginary part; however for this case, because S 1 has a smaller attenuation it will have a significant contribution to both transmission and reflection.

Multiple effective wavenumbers
An important conclusion from det G(Sp) = 0 is that the wavenumbers Sp are independent of the angle of incidence θ inc . We focus on showing the results for Sp, rather than sp, because then we do not need to specify θ inc .
As a specific example, let us consider circular particles with Dirichlet boundary conditions (i.e. particles with zero density or soundspeed), and the parameters , kb = 1.001, kao = 0.5, φ = 30%, where ao is the radius of the particle. With the above parameters, we found that truncating the matrix Ψ (s), with |n| ≤ 3 and |m| ≤ 3 in (3.15-3.17), led to accurate results when calculating the effective wavenumbers Sp, i.e. the roots of (3.30). Numerically calculating the wavenumbers Sp then leads to Figure 3.
The effective wavenumbers with the lowest attenuation (smallest imaginary part) contribute the most to the transmitted wave. In Figure 3 we see two wavenumbers have lower attenuation then the rest, both within the dashed grey circle. The blue point represents the wavenumber that most of the literature focuses on calculating: it has a positive real part and therefore propagates forwards along the x−axis (into the material) as is expected for a transmitted wave. However, the other wavenumber, with negative real part, is equally as important because it actually has lower attenuation. Figure 1 illustrates several effective wavenumbers, some travelling forward into the material, while others have negative phase direction (travel backwards).
In Figure 3 we see what appears to be an infinite sequence of effective wavenumbers Sp, where |Sp| → ∞ as p → ∞. To confirm their existence, and to find their locations as |p| → ∞, we develop asymptotic formulas in appendix B. The results of the asymptotics are summarised below.
For monopole scatterers, where n = m = 0 in (3.15), equations (A 7) give the effective wavenumbers S o p at leading order: 3) and for any integer p. We use the superscript "o" to distinguish these wavenumbers for monopole scatterers from others. Even though (5.3) was deduced for large integer p, it gives remarkably agreement with numerically calculated wavenumbers, except for the two lowest attenuating wavenumbers, as shown in Figure 4. In the figure we denoted S o± * as the effective wavenumber that can be calculated by low volume fraction expansions [11,14]. The asymptotic formula is surprisingly accurate except for the two lowest attenuating wavenumbers. The wavenumber S o * can be calculated by using low volume fraction expansions [11].
For multipole scatterers, where both n and m could potentially range from −∞ to ∞ in (3.15), we can also calculate an infinite sequence of effective wavenumbers. To show this explicitly, we consider the limit of large bk, with |k| ∼ |S|. In the opposite limit bk 1, the Rayleigh limit, only one effective wavenumber is required [44,45].
At leading order, the asymptotic solution of (A 11) leads to the effective wavenumbers: for integer p. This confirms that there are an infinite number of effective wavenumbers for large scatterers, i.e. bk 1. The distribution of these wavenumbers is similar to the monopole wavenumbers shown in Figure 4. These asymptotic formulas (5.3) and (5.5) demonstrate the existence of multiple effective waves in the limit of small (monopole and Dirichlet) scatterers (5.3) and large scatterers (5.5). However, neither of these formulas, nor the low volume fraction expansions of the wavenumber [11], are able to accurately estimate the low attenuating backward travelling effective wavenumber such as S 1 shown in Figure 3  given above). There is currently no way to analytically estimate these types of wavenumbers, even though they are necessary to accurately calculate transmission due to their small attenuation. The only approach it seems is to numerically solve (3.30).

Numerical results
Here we present numerical results for monopole scatterers, as these have explicit expressions for reflection (4.5) and the transmitted wave (4.4) (or more accurately the average scattering coefficients). We compare our analytic solution with a classical method that assumes only one effective wavenumber [11,13], and the Matching Method [18], recently proposed by the authors.
It should be noted that all of these approaches aim to solve the same equation (2.14).
Note that for monopole scatterers, using only one effective wavenumber s 1 can, in some cases, lead to accurate results. However, for multipole scatterers (a more common scenario practically) this is rarely the case because, as shown by Figure 3, there can be at least two effective wavenumbers with low attenuation, and therefore both are needed to obtain accurate results.
For the numerical examples we use the parameters which implies that the number fraction n ≈ 0.38 per unit area. When we choose to fix the wavenumber, as we do for Figure 5 and Figure 6, we use bk = 1.001. This leads to a wavelength (2π/k) which is roughly six times larger than the particle diameter. If the particle was, say, more than a hundred times smaller than the wavelength, then only one effective wavenumber in the sum (4.4) would be necessary to accurately calculate A 0 (X). Method also accounts for multiple effective wavenumbers, and is described in [18]. The low volume fraction method assumes a low volume fraction expansion for just one effective wavenumber [11]. The small graph on the right is a magnification of the region around x = 0. Close to the boundary x = 0, both A 1 0 e is 1 x and the low volume fraction method are inaccurate, which would potentially lead to inaccurate predictions for transmission and reflection. To start we compare the average scattering coefficient A 0 (x) calculated by the Wiener-Hopf solution (4.4) with other methods in Figure 5. The most accurate of these other methods is the Matching Method [18,46], and it closely agrees with the Wiener-Hopf solution when using 352 effective wavenumbers. The exception is the region close to the boundary x = 0, where the Wiener-Hopf solution experiences a rapid transition. The low volume fraction method is the most commonly used in the literature: it assumes a small particle volume fraction 10 and just one effective wavenumber [11,13]. One significant conclusion we can draw from Figure 5 is that both the low volume fraction method and A 1 0 e is1x are inaccurate near the boundary x = 0. This means that both of these methods lead to inaccurate reflection coefficients.
In general, the Wiener-Hopf method does not lead to an explicit formula for the reflection coefficient (3.35), because we do not have an exact factorisation (3.23) for any truncated square matrices. However, there are methods [13,18,20,29,30,34] to calculate An(x), from which we can obtain the reflection coefficient (2.12). The method [18] also accounts for multiple effective wavenumbers. So one important question is: when using (2.12), how many effective wavenumbers do we need to obtain an accurate reflection coefficient?  Figure 6. Demonstrates, with a log-log graph, how increasing the number of effective waves P leads to a more accurate reflection coefficient R P , when using (4.6). The non-dimensional wavenumber is kb = 1.001, and the other parameters used are given by (6.1), with their definitions explained by (2.4-2.9). Here R is the reflection coefficient given by (4.5).
The error |R P − R| continuously drops as P increases because of the rapid transition that occurs to A 0 (x) near the boundary x = 0, see Figure 5. However, methods such as the Matching Method [18] are able to accurately calculate the reflection coefficient without taking into account this rapid transition.
In Figure 6 we show how increasing the number of effective waves P reduces the error between R P (4.6) and R (4.5). To calculate a highly accurate reflection coefficient R, we could use either (4.5) or the Matching Method [18,46], as both give approximately the same R. Now we ask: how does the reflection coefficient (4.6), deduced via the Wiener-Hopf technique, compare with other methods across a broader range of wavenumbers. The result is shown in Figure 7, where R O is a low volume fraction expansion 11 of just one effective wavenumber [13]. The reflection coefficient R M is calculated from the Matching Method [18,46]. The general trend is clear: R O becomes more inaccurate as we increase the background wavenumber kb. On the other hand both R M and R agree closely over all k.
One result to note is the "instability" exhibited by the Wiener-Hopf solution near the boundary x = 0, see Figure 5. This instability occurs because we represented A 0 (x) as a superposition of truncated waves, which is only accurate as long as the discarded terms are small. So, for 10 For the low volume fraction method we used a small volume fraction expansion for the wavenumber, but we numerically evaluated the wave amplitude. This is because the alternative, a small volume fraction expansion of the wave amplitude, led to poor results. 11 We use the reflection coefficient [13, equation (39)], rather than the explicit low volume fraction expansion [13, equation (40 -41) ]. This is because using equations (40 -41)  ). Here R is given by the Wiener-Hopf solution (4.5), R O uses a low volume fraction expansion of just one effective wavenumber [13], and R M is calculated from the Matching Method [18].
a truncation number P , we can expect the instability to occur when e is P x is not small, i.e.
x ≈ 1/Im sp. However, this instability does not affect the accuracy of the reflection coefficient (4.5) deduced by the Wiener-Hopf technique, as demonstrated by close agreement with the Matching Method in Figure 7.

Conclusion and Next Steps
The major result of this paper is to prove that the ensemble-averaged field in random particulate materials consists of a superposition of waves, with complex effective wavenumbers, for one fixed incident wavenumber. These effective wavenumbers are governed by the dispersion equation (5.1) and are independent of the angle of incidence θ inc . We showed asymptotically in section 5 that this has an infinite number of solutions, and hence there are an infinite number of effective wavenumbers. The Wiener-Hopf technique also provides a simple and elegant expression for the reflection coefficient (3.35), whose form can be used to guide and assess methods to characterise microstructure [47,48].
To numerically implement the Wiener-Hopf technique, we considered particles that scatter only in their monopole mode in section 6. There we saw that when close to the interface of the half-space, a large number of effective wavenumbers were necessary to reach accurate agreement with an alternative method from the literature, the Matching Method as introduced by the authors in [18]. To obtain a constructive method via the Wiener-Hopf technique for general scatterers, and not just monopole scatterers, will require the factorisation of a matrix-function [31], which is challenging. For these reasons the Matching Method [18] is presently more effective than using the Wiener-Hopf technique. However, there is ongoing work to use approximate methods [33,34,49] which exploit the symmetry and properties of the matrix (3.15).
Moving forwards, this paper together with [18], establish accurate and robust solutions to the governing equation (2.14). These same methods can now be translated to three spatial dimensions and vectorial waves (e.g. elasticity and electromagnetics), with much of the groundwork already available [12,16,40]. Some clear challenges, that can now be addressed, are to verify the accuracy of the statistical assumptions used to deduce (2.14). These include the hole-correction and the quasicrystalline approximations. As these are now the only assumptions used, we could compare the solution of (2.14) with multipole methods [50,51] in order to investigate their accuracy and limits of validity.
We have now rewritten the abstract to accommodate these points, though are limited due to space; i.e. to state the problem and conclusions more clearly, mentioned the matching method, and that the reflection coefficient, deduced by the Wiener-Hopf technique, can be explicitly evaluated for monopole scatterers.
3. Page 4. The spatial extent of the incident plane wave is not mentioned as an assumption that simplifies the analysis. A plane wave with spatial extent in y would encounter different boundary conditions along the boundary as the relative acoustic impedance changes, due to local changes in the particulate material. Also near and far field effects would need to be considered in the analysis.
We agree that having an incident plane wave of the form e i(αx+βy) for x < 0 and y ∈ R does indeed simplify the analysis as long as the statistics of the distribution of particles is invariant in the orthogonal y-direction. Clearly there will be local changes along the boundary for any specific configuration of particles, but the act of ensemble averaging removes these local changes, and leads to the field u(x, y) having a translational symmetry along y, see [1] for details. We have added a footnote just below equation (2.10) on this.
We are not sure which specific near and far field effects the referee refers to. There are a few different fields which could be considered near-fields: when two particles are close to each other we could describe their near field effects, which are included in the governing equation (2.14). We could also consider the fields near the boundary x = 0 to be near-fields, which is where these multiple effective waves have the greatest affect. These types of near-field effects, and fields far from the boundary, are all considered in our analysis. Only if the region containing the inclusions was finite in y would we need to consider global far-field effects in our work (i.e. far from the multiple scattering region); but as already mentioned, we preclude this in our analysis.