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.

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 characterize 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-8]. Reducing an 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-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 half-space filled with particles (figure 1).

Figure 1.

When an incident plane wave ei, with k = k(cosθinc, sinθinc), encounters an (ensemble-averaged) particulate material, it excites many transmitted plane waves and one reflected plane wave. The transmitted waves are of the form ei with wavenumbers s = S(cosθ, sinθ), where both S and θ are complex numbers. The larger Im s, 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. (Online version in colour.)

Figure 3.

Examples of effective wavenumbers S which satisfy the dispersion equation (5.1) with the properties (5.2). 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 S1 and S2, 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 S2, as it often has the lowest imaginary part; however, for this case, because S1 has a smaller attenuation it will have a significant contribution to both transmission and reflection. (Online version in colour.)

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 scatter 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, (ch. 5. Wiener–Hopf Technique)] and [20,21] for an introduction. It is especially useful for semi-infinite domains [22-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 factorization of a matrix function. To explicitly perform this factorization is difficult [28-31]. Indeed, this is often the hardest aspect of employing the Wiener–Hopf technique, although there exist approximate methods for this purpose [28,32-34]. We do not focus in this article on these analytic factorizations, 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.

Figure 1 shows the main set-up and result of this paper: an incident plane wave excites the half-space x > 0 filled with ensemble-averaged particles (the blue region), which generates a reflected wave and many effective transmitted waves. The 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 summarizing the equations that govern ensemble-averaged waves in two dimensions in §2. Following this, in §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 §4, we specialize 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 §3, admits an infinite number of solutions, the effective wavenumbers. In §5, we deduce asymptotic forms for the effective wavenumbers in both a low- and high-frequency limit. In §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 §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: and where k and k 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-averaged 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 uinc(x, y) and waves scattered by each particle, the jth scattered wave being u(x, y), we can write A simple and useful scenario to consider is when all particles are placed only within the half-space[1] x > 0, which are then excited by a plane wave, with implicit time dependence ei, 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 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: 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 a be the maximum distance from the particle's centre to its boundary, then we can set b = γa, where γ≥2 so as to avoid two particles overlapping. The volume fraction that does not include the exclusion zone ϕ′, as used in [18, (eqn (4.7))], is then ϕ′ = 4ϕ/γ2.

The T are the coefficients of a diagonal T-matrix [35-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 details of this procedure are omitted here for brevity. To represent the ensemble-averaged scattered wave from a particle, whose centre is fixed at (x1, y1), we use for , so that (X, Y ) is on the outside of this particle, with (R, Θ) being the polar coordinates of (X, Y ), H(1) are Hankel functions of the first kind, and A 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 eqn (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 ψ is given by Later we show that, as expected, is independent of x.

The system governing A(x) is given by eqn (4.7) of [18]: for all integers m. Kristensson [40, (eqn (15))] presents an equivalent integral equation for electromagnetism and particles in a slab.

Our main aim is to reach an exact solution for A(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 transform can lead to elegant and simple solutions. To facilitate, we must analytically extend (2.14) for all by defining for integers m, where if the A(x) were known for x > 0, then the D(x) would be given from the left-hand side. Note that the kernel ψ defined in (2.13) is already analytic in the domain .

The field D0(x) is not just an abstract construct, it is closely related to the reflected wave: by directly comparing (1) with the reflection coefficient (2.12), for x <  − b, we find that

To solve (1), we employ the Fourier transform and its inverse, which we define as for any smooth function f. We then define

We can determine where and are analytic by assuming[6] that and for some (possibly small) positive constant c. This leads to being analytic for Im s >  − c, while is analytic for Im s < c. In other words, both and are analytic in the overlapping strip

To apply the Wiener–Hopf technique, we also need to specify the large s behaviour for both and . To achieve this, we assume, on physical grounds, that A(x) is bounded when x → 0+, and D(x) is bounded when x → 0−. Then, it can be shown [21,41] that in their respective half-planes of analyticity.

Applying a Fourier transform to both sides of equation (1), the left-hand side becomes in which is well defined (i.e. analytic) for s in the strip: see appendix A for details. The right-hand side of (1) becomes 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

If (3.12) is satisfied, then we can combine (3.9), (3.11) and (A 6), to obtain the Fourier transform of (1) in matrix form: where and are vectors with components and , respectively, and where, for reference, and is given by (A 6). In the above, θ and S are chosen to satisfy Later we identify S and θ 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 θ, we take sgn(Re s) = sgn(Re S) which together with (3.19) leads to where both S and θ, when considered as functions[7] of s, contain branch-points at s =  ± iksinθinc with finite branch-cut running between −iksinθinc and iksinθ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 will be a key step in solving (3.13), and so the following identities will be useful: and where (3.21) results from (3.18) and (A 10). Equation (3.22) then follows from using (3.21)1, (3.15) and appendix C.

Multiple waves solution

To solve (3.13), we use a matrix product factorization [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 factorization exists [42], as this will lead to a proof that (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 s >  − ϵ) and another analytic in the lower half-plane. This is achieved below where we define so that −(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 (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 (s) for all s, i.e. (s) is entire.

To determine (s), we need to estimate its behaviour as |s| → ∞. From (3.8), we have that as |s| → ∞ in the upper half-plane, and from (3.15) to (3.17): for s in the strip (3.12). From this, we know that the factors +(s) and −(s) must be as |s| → ∞, in their respective half-planes of analyticity [28]. So, the left-hand side of (3.26) behaves as as |s| → ∞ in Im s >  − ϵ. We can therefore use Liouville's theorem to conclude that (s)≡0, which means the Wiener–Hopf equation (3.26) is formally equivalent to and Let C+(s) be the cofactor matrix of +(s), so that From the property (3.22)1, we can write for some function f. Then, for every root s = s of , with Im s > 0, we have that −s is also a root, and vice-versa. From here onwards, we assume: For any truncated matrix (s), i.e. evaluating m, n =  − M, …, M in (3.15), the roots s are discrete. In §5, we demonstrate asymptotically that they are indeed discrete for the limits of low and high wavenumber k. For the numerical results presented in this paper, we numerically solve the above dispersion relation for the truncating the matrix (s), and then increase M until the roots converge (typically no more than M = 4 was required).

Given , every root of must either be a root of or a root of . For [+(s)]−1 to be analytic in the upper half-plane, must only have roots s =  − s. As a consequence, only has roots s = s.

To use the residue theorem below, we need to calculate for s close to the root −s, in the form where we use ( instead of (, because it is more difficult to numerically evaluate .

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 For x > 0, the integral over s∈[ − ∞, ∞] in (3.32) is, by Jordan's lemma, the same as a clockwise integral over the closed contour which surrounds the poles −s1, −s2, …, 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.

Figure 2.

An illustration of the contour integral over , used to calculate (3.34) for x < 0, and the contour integral over , used to calculate (3.32) for x > 0. The −s (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). (Online version in colour.)

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 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)−(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 (x) is a sum of plane waves with different effective wavenumbers s, each satisfying (3.30). In §5, we discuss these roots in more detail, and in §4, we see that usually only a few effective wavenumbers are required to obtain accurate results.

Reflection coefficient

Substituting (3.34) in (3.2) leads to Alternatively, the reflection coefficient can be calculated from (2.12) by employing the form of (x) from (3.32), which is the more common approach. To simplify, we use which then implies that ψ(x1 − x) = (2/α)ie−iei for x1≥x. The above is shown in [43, (eqn (37))] and [11, (eqn (65))]. This result together with (3.32) substituted into (2.12) leads to the form where we used that Im s > 0. The above agrees with [13, (eqn (39))] and[8] [18, (eqn (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 factorization (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 factorization (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 §6.

Wiener–Hopf factorization

For scalar problems, there are well-known techniques to factorize Ψ00(s) = Ψ−00(s)Ψ+00(s), such as Cauchy's integral formulation, for details see [19, (Section 5. Wiener–Hopf Technique)] and [21].

For monopole scatterers we use S2 − k2 = s2 − α2 and rewrite with N0(bS) given by (3.17). Then, because q(s) → 1 as |s| → ∞, we can factorize q(s) = q−(s)q+(s) using and 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 = s1, s2, …, s where as p increases, the s 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.

Explicit solution for monopole scatterers

For monopole scatterers A(x) = D(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, (3.32) reduces to for x > 0, where we used (4.3), C+(s) = 1,  = iT0, 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 then indicates how many effective waves are needed to accurately describe the field near the boundary x = 0.

Multiple effective wavenumbers

Equation (3.32) clearly shows that (x) is a sum of attenuating plane waves, each with a different effective wavenumber s. These s satisfy the dispersion equation (3.30): with given by (3.16) and the first identity follows from (3.22).

An important conclusion from is that the wavenumbers S are independent of the angle of incidence θinc. We focus on showing the results for S, rather than s, 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 where a 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 S, i.e. the roots of (3.30). Numerically calculating the wavenumbers S 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 S, where |S| → ∞ as p → ∞. To confirm their existence, and to find their locations as |p| → ∞, we develop asymptotic formulae in appendix B. The results of the asymptotics are summarized below.

For monopole scatterers, where n = m = 0 in (3.15), equations (B 7) give the effective wavenumbers S at leading order: and 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 remarkable agreement with numerically calculated wavenumbers, except for the two lowest attenuating wavenumbers, as shown in figure 4. In the figure, we denoted S* as the effective wavenumber that can be calculated by low volume fraction expansions [11,14].

Figure 4.

Comparison of the asymptotic formula (5.3), which predicts an infinite number of effective wavenumbers, with numerical solutions for the effective wavenumbers (5.1). The parameters used are given by (5.2), with their definitions explained in (2.4)–(2.9). Here we chose b = 1.0, so the non-dimensional wavenumbers bS are the same as shown. The asymptotic formula is surprisingly accurate except for the two lowest attenuating wavenumbers. The wavenumber S* can be calculated by using low volume fraction expansions [11]. (Online version in colour.)

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 (B 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 formulae (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 formulae, nor the low volume fraction expansions of the wavenumber [11], are able to accurately estimate the low attenuating backward travelling effective wavenumber such as S1 shown in figure 3 (in this case not related to the S1 and S1 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 s1 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 per unit area. When we choose to fix the wavenumber, as we do for figures 5 and 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 100 times smaller than the wavelength, then only one effective wavenumber in the sum (4.4) would be necessary to accurately calculate A0(X).

Figure 5.

Compares the absolute value of the average field A0(x) calculated by different methods. The field A0(x) is closely related to the average transmitted wave [13]. The non-dimensional wavenumber kb = 1.001, the other parameters are given by (6.1), with their definitions explained by (2.4)–(2.9). Using the Wiener–Hopf solution (4.4), we approximate A0(x) by using either 352 effective wavenumbers s1,  s2,  …, s352, or just 1 effective wavenumber s1. The Matching 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 A10ei and the low volume fraction method are inaccurate, which would potentially lead to inaccurate predictions for transmission and reflection. (Online version in colour.)

Figure 6.

Demonstrates, with a log–log graph, how increasing the number of effective waves P leads to a more accurate reflection coefficient , 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 is the reflection coefficient given by (4.5). The error continuously drops as P increases because of the rapid transition that occurs to A0(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. (Online version in colour.)

To start we compare the average scattering coefficient A0(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 A10ei 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 factorization (3.23) for any truncated square matrices. However, there are methods [13,18,20,29,30,34] to calculate A(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?

In figure 6, we show how increasing the number of effective waves P reduces the error between (4.6) and (4.5). To calculate a highly accurate reflection coefficient , we could use either (4.5) or the Matching Method [18,46], as both give approximately the same .

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 is a low volume fraction expansion[11] of just one effective wavenumber [13]. The reflection coefficient is calculated from the Matching Method [18,46]. The general trend is clear: becomes more inaccurate as we increase the background wavenumber kb. On the other hand, both and agree closely over all k.

Figure 7.

Compares different methods for calculating the reflection coefficient when varying the non-dimensional wavenumber kb. The other parameters used are given by (6.1), with their definitions explained in (2.4)–(2.9). Here is given by the Wiener–Hopf solution (4.5), uses a low volume fraction expansion of just one effective wavenumber [13], and is calculated from the Matching Method [18]. (Online version in colour.)

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 A0(x) as a superposition of truncated waves, which is only accurate as long as the discarded terms are small. So, for a truncation number P, we can expect the instability to occur when ei is not small, i.e. x ≈ 1/Im s. 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 §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 characterize microstructure [47,48].

To numerically implement the Wiener–Hopf technique, we considered particles that scatter only in their monopole mode in §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 factorization 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], establishes 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 quasi-crystalline 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.