## Abstract

We formally deduce closed-form expressions for the transmitted effective wavenumber of a material comprising multiple types of inclusions or particles (multi-species), dispersed in a uniform background medium. The expressions, derived here for the first time, are valid for moderate volume fractions and without restriction on the frequency. We show that the multi-species effective wavenumber is not a straightforward extension of expressions for a single species. Comparisons are drawn with state-of-the-art models in acoustics by presenting numerical results for a concrete and a water–oil emulsion in two dimensions. The limit of when one species is much smaller than the other is also discussed and we determine the background medium felt by the larger species in this limit. Surprisingly, we show that the answer is not the intuitive result predicted by self-consistent multiple scattering theories. The derivation presented here applies to the scalar wave equation with cylindrical or spherical inclusions, with any distribution of sizes, densities and wave speeds. The reflection coefficient associated with a halfspace of multi-species cylindrical inclusions is also formally derived.

## 1. Summary

Materials comprising mixtures of diverse particles, inclusions, defects or inhomogeneities dispersed inside a background medium arise in a wide range of applications, including composite materials, emulsions, gases, polymers, foods and paints. We will refer to these as multi-species materials. Of great importance is the ability to characterize these materials and their *microstructure*, such as particle size distribution and volume fractions. One approach to do this is to employ waves, including electromagnetic, acoustic and elastodynamic waves. If either the receivers are much larger than the inclusions, or the wavelength is much longer than the inclusions, then the receivers will measure the *ensemble-averaged* properties of the wave [1]. This includes the wave speed, attenuation and reflection. Even methods that estimate fluctuations of the wave on smaller scales, such as the averaged intensity, often require the ensemble-averaged wave properties as a first step [2–4]. So in order to improve material characterization, or to design materials with tailored properties, a crucial step is to rigorously calculate the sound speed and attenuation for multi-species materials.

In this paper, we present and formally deduce the effective wavenumber and reflected field of a plane wave scattered by a material comprising different families, or *species*, of particles with distributions of sizes and properties. The work here differs from the existing literature as our results are not limited to low frequencies and are valid for moderate number density. This is achieved by extending the methods introduced in [5] for calculating the effective transmission into a halfspace of a single-species material.

Our approach does not rely on an extinction theorem or the manipulation of divergent integrals or series. The one assumption that is employed is the *quasi-crystalline approximation* [6]. For a single species, this approximation is supported by theoretical [7,8], numerical [9] and experimental [10] evidence; however, the authors are unaware of any rigorous bounds for the error introduced by this approximation. We remark that the quasi-crystalline approximation (QCA) makes no assumption on the material properties, so in principal it is consistent for weak scattering, low or high frequency, or in dense or sparse mixtures. It is also exact when there is only one possible configuration for the particles, for example when the particle centres lie on the coordinates of a crystal lattice (Lax [11], where QCA is discussed under (4.3)). For simplicity, we also restrict attention to the case of circular cylindrical or spherical particles, although our methods can be extended to the case of general-shaped particles by using Waterman’s T-matrix approach, e.g. [12–14].

In the context of electromagnetic wave scattering, methods for predicting wave propagation and reflection for multi-species material have previously been developed [3,4,15]. These models have been useful for interpreting data from remote sensing, although it appears that such models cannot systematically reproduce experimental results [16]. In numerous contexts, but particularly in the context of electromagnetics, the standard approach is to employ the Lippman–Schwinger formulation [17,18]. However, such a formulation is restrictive as it is not valid for magnetic media in the electromagnetism context or for scatterers with varying density in acoustics, as identified in [19]. Although it is possible to extend the Lippman–Schwinger formulation to account for this effect [19], we found it simpler to extend the multiple scattering theory [2,20,21].

Our approach is also in contrast to coupled-phase theory where the first step is to estimate the ensemble average of the governing equations [22], without explicitly considering multiple scattering. Although this method can accommodate hydrodynamic interactions and has been extended to polydisperse inclusions (multi-species) it does not completely capture multiple scattering [23,24].

A suggestion for calculating the multi-species effective wavenumber came from Waterman & Truell, eqn (3.25a) in the conclusion of [25]. Their suggested formula has been extensively employed in acoustics, e.g. [26–28]. However, their formula is only valid for low frequency and dilute distributions of particles [21], so it does not properly account for multiple scattering. The approach in [25] combined with [29] led eventually to the state-of-the-art models for the effective acoustic wavenumber in colloidal dispersions [28]. We numerically compare our results with these authors.

Given an overall particle number density *k*, our main results for a multi-species material comprising *circular cylinders* are the effective transmitted wavenumber:
*u*_{in}=*e*^{iαx+iβy}, with *f*_{°}〉 and 〈 *f*_{°°}〉 are defined in (4.12) and (4.13). The formula (1.2) is briefly deduced in §7a, and in figure 7 we give a pictorial representation, although we stress that the choice *θ*_{ref}=*π*−2*θ*_{in} is not due to a simple geometric argument, but appears from rigorous derivations. From the reflection coefficient (1.2), it is possible to choose effective material properties [30]. However, because the reflection coefficient depends on the angle of incidence via 〈 *f*_{°}〉(*θ*_{ref}) and 〈 *f*_{°°}〉(*θ*_{ref}), it is likely that these effective material properties change with the angle of incidence.

In the electronic supplementary material, we provide a brief self-contained version of these formulae, and the corresponding result for spherical particles, both for a finite number of species. We also provide open source code that implements these formulae [31]. For spherical inclusions, the effective transmitted acoustic wavenumber becomes
*F*_{°}〉 and 〈*F*_{°°}〉 are functions associated with scattering from the spherical particle and are defined in (A.2). Note that 〈*F*_{°°}〉 has no *θ* dependency. For a longer discussion of multiple scattering from spheres, see [30].

By developing *multi-species* formulae valid for higher number densities and frequencies, we open up the possibility of characterizing and designing a wide range of advanced materials. The effects of multiple scattering appear only for moderate number density, i.e. in the term 〈 *f*_{°°}〉(0) in (1.1) and 〈*F*_{°°}〉 in (1.5). One important consequence of this multiple scattering term is that a multi-species material can exhibit properties not exhibited by that of the background medium with only one constituent species. We stress that even for just two types of circular cylindrical particles, the effects of multiple scattering are neither intuitive nor easily deduced from the single-species case. This becomes apparent in the simple example of a multi-species material where one species is much smaller than the other. In this scenario, we compare our expression for the multi-species effective wavenumber with the state-of-the-art models from acoustics [28] and a *self-consistent* type approximation [32–34], which can be calculated from the single-species formula via an iterative approach: first one *homogenizes* the small particle and background mixture before considering the multiple scattering of the larger particles in the new (homogenized) background medium. We show analytically that this naive self-consistent methodology is not even correct in the low-frequency limit. This is then demonstrated numerically for the cases of an emulsion and concrete. Our results are in in-line with [35], who discuss an effective medium model of a three-phased material in the low-frequency limit.

The outline of this paper is as follows. In §2, we describe the exact theory of multiple scattering for *N* cylinders of any radius, density and sound speed. From there we calculate the effective (ensemble-averaged) equations and apply statistical approximations in §3. In §4, we deduce the governing system for the effective wavenumbers at arbitrary total number density and arbitrary frequency, before specializing the result to the case of moderate number fraction and low frequency. In §5, we investigate the specific, representative case of two types of circular cylindrical species and compare different approximations graphically. To calculate the reflected or transmitted wave we also need the effective amplitude, which we calculate in §7 followed by the effective reflected wave. We close in §8, where we discuss avenues for improvement of the techniques and more general further work.

## 2. Multipole method for cylinders

In this section, we describe the exact theory for scalar multiple wave scattering from a finite number *N* of circular cylinders possessing different densities, wave speeds and radii. Parameters associated with the medium are summarized in table 1. Naturally, the system of equations describing this problem bears strong similarities to that obtained by Záviška (see references in [36] and in [5] for the single species circular cylindrical particle context). Assuming time-harmonic dependence of the form *e*^{−iωt}, the pressure *u* outside all the cylinders satisfies the scalar Helmholtz equation
*j*th cylinder the pressure *u*_{j} satisfies
^{2} is the two-dimensional Laplacian and
**x**_{j} is the centre of the *j*th cylinder and **x**=(*x*,*y*) is an arbitrary point with origin *O*. (See figure 1 for a schematic of the material properties and coordinate systems.) Then we can define *u*_{j} as the scattered pressure field from the *j*th cylinder,
*H*_{m} are Hankel functions of the first kind, *q*_{j}=(*ρ*_{j}*c*_{j})/(*ρc*). In the limits *q*_{j}→0 or

The pressure outside all cylinders is the sum of the incident wave *u*_{in} and all scattered waves,
*j*th cylinder is
*J*_{m} are Bessel functions of the first kind. The arbitrary constants *j*th cylinder *R*_{j}=*a*_{j}. The boundary conditions of continuity of pressure and normal velocity on the cylinder boundaries are given, respectively, by
*ρ* and *ρ*_{j} denote the material densities of the background and of the *j*th cylinder respectively. To impose the boundary conditions, we now express the relevant fields in terms of the (*R*_{j},*Θ*_{j}) coordinate system. For the incident wave

where *I*_{j}=*u*_{in}(*x*_{j},*y*_{j}) following the Jacobi–Anger expansion [37]. For the scattered waves (2.4), we use Graf’s addition theorem (9.1.79) in [38],
*R*_{ℓj},*Θ*_{ℓj}) is the polar form of the vector **x**^{j}−**x**^{ℓ}. Using the above and (2.9) we can impose the boundary conditions (2.8) to arrive at the following system of equations
*j*=1,…,*N* and all integers *m*. Furthermore, the coefficients associated with the pressure inside the cylinder (2.7) are then given by
*u*(*x*,*y*) is entirely prescribed.

In any given material, it is impossible to know the exact position and properties of all constituent particles. Our goal is, therefore, to solve (2.11) not for one particular configuration of scatterers, but instead to calculate the average value of the coefficients *effective* field describes the *ensemble-averaged* field that is usually measured in an acoustic experiment, as the receiver face is typically much larger than the particles and the distance between them [1,39]. In our case, we obtain an ensemble average by averaging over all particle configurations and all the material properties of the particles. This approach is general and can be tailored to different scenarios, e.g. when detailed information *is* known about the particle material properties.

## 3. Averaged multiple scattering

For an introduction to ensemble averaging of multiple scattering, see [2,40], where the result for a classical dilute isotropic mixture was determined. Here we present a brief self-contained explanation tailored to multi-species.

Consider a configuration of *N* circular cylinders centred at **x**_{1},**x**_{2},…,**x**_{N} with the scattering properties **s**_{1},**s**_{2},…,**s**_{N}, where **s**_{j} denotes the properties of the *j*th cylinder, i.e. here these are **s**_{j}=(*a*_{j},*ρ*_{j},*c*_{j}). Each **x**_{j} is in the region **s**_{j} are taken from the set *a*_{j}∈[0,1], *ρ*_{j}∈[1,2] and *c*_{j}∈[100,200].

The probability of the cylinders being in a specific configuration is given by the probability density function *p*({**x**_{1},**s**_{1}},{**x**_{2},**s**_{2}},…,{**x**_{N},**s**_{N}}). Using the compact notation *Λ*_{i}={**x**_{i},**s**_{i}} to denote the properties of the *i*-th cylinder, it follows that
*both* **x**_{j}) and **s**_{j}) with *d**Λ*_{j}=*d***x**_{i} *d***s**_{i}. Note that *p*(*Λ*_{1},*Λ*_{2}) is the probability of one cylinder having the properties *Λ*_{1} and another having the properties *Λ*_{2}, when the properties of all the remaining *N*−2 cylinders are unknown. And as the cylinders are indistinguishable: *p*(*Λ*_{1},*Λ*_{2})=*p*(*Λ*_{2},*Λ*_{1}). Furthermore, we have
*p*(*Λ*_{1},…,*Λ*_{N} | *Λ*_{j}) is the conditional probability of having cylinders with the properties *Λ*_{1},…,*Λ*_{N} (not including *Λ*_{j}), given that the *j*th cylinder has the properties *Λ*_{j}. Likewise, *p*(*Λ*_{1},…,*Λ*_{N} | *Λ*_{ℓ},*Λ*_{j}) is the conditional probability of having cylinders with the properties *Λ*_{1},…,*Λ*_{N} (not including *Λ*_{ℓ} and *Λ*_{j}) given that there are already two cylinders present, with properties *Λ*_{ℓ} and *Λ*_{j}.

Given some function *F*(*Λ*_{1},…,*Λ*_{N}), we denote its average, or *expected value*, by
*j*th cylinder, *Λ*_{j} and average over all the properties of the other cylinders, we obtain a *conditional average* of *F* given by
*Λ*_{j}. The average and conditional averages are related by
*F*〉_{ΛℓΛj} is the conditional average when fixing both *Λ*_{j} and *Λ*_{ℓ}, and 〈*F*〉_{ΛℓΛj}=〈*F*〉_{ΛjΛℓ}.

Returning to the task of obtaining effective properties for a multi-species medium, we multiply the system (2.11) by *p*(*Λ*_{2},…,*Λ*_{N}|*Λ*_{1}) and average over *Λ*_{2},…,*Λ*_{N}, to reach
*j*=1, used the conditional average definition (3.2b) and defined **s**_{ℓ} explicit. To further simplify the above, note that all terms in the sum over ℓ give the same value. That is, the terms in the integrand depend on ℓ solely through the dummy variable *Λ*_{ℓ}. In particular, the probability distribution is the same for each cylinder, and if *Λ*_{2}=*Λ*_{l}, then *Λ*_{ℓ}. We use this to obtain
*p*(*Λ*_{2} | *Λ*_{1}) and

With **x** outside *u*_{in}(*x*,*y*)〉=*u*_{in}(*x*,*y*) because the incident field is independent of the scattering configuration. We can then rewrite the average outgoing wave *u*_{j} by fixing the properties of the *j*th cylinder *Λ*_{j} and using equation (3.2a) to reach
**x**, we obtain

### (a) Statistical approximations

In order to solve (3.6) for *p*(*Λ*_{2} | *Λ*_{1}). In this work, we adopt the standard closure approximation for single species, but extended to multi-species, the QCA [5,6]:
*Λ*_{1} by its expected value in *Λ*_{1}. Note also that the expected difference in *Λ*_{2}:

Using QCA, we introduce the notation

Next, we determine a suitable approximation for the pair distribution *p*(*Λ*_{2}|*Λ*_{1}), beginning with (3.2a) to write
**X**_{1},**X**_{1},…,**X**_{N} and the scattering property random variables **S**_{1},**S**_{1},…,**S**_{N}, and write probability density functions in the form, e.g.
*P*(**S**_{1}=**s**_{1}) is the probability density in **s**_{1}. The above assumes that **x**_{1} of the cylinder is independent of the scattering property **s**_{1}. This is not always the case, for example, depending on the size of the cylinder, some positions near the boundary of

For the remaining distribution in (3.12), we use
*R*_{21}:=∥**x**_{1}−**x**_{2}∥, *a*_{21}=*b*_{1}+*b*_{2} for some *b*_{1}≥*a*_{1} and some *b*_{2}≥*a*_{2}, and *b*_{1} is the radius of exclusion around **x**_{1} which is usually chosen to be proportional to the radius *a*_{1}. Note that when integrating (3.17) above in **x**_{1} and **x**_{2}, we obtain

Ultimately, substituting (3.16) and (3.17) into (3.15) in tandem with (3.14) leads to the pair distribution
*H*(*x*) denotes the Heaviside function, under the assumption

We now include a discussion of other commonly used pair distributions. We remark that for densely packed scatterers, other pair distributions [41] are preferred and take the form
**x**_{1} and **x**_{2} to no longer effect each other. One popular choice for *χ* is the Percus–Yevick function, which assumes all scatterers are uniformly randomly distributed [44], though *χ* can also be used to specify if some species are more likely to be closer or further apart.

In this work, we set *χ*=0 for simplicity (unless otherwise stated), but also because it is not clear that the error introduced by using *χ*=0 is in any way greater than the error committed due to QCA (3.10). Both the hole correction (3.17) and QCA (3.10) make similar assumptions: for *R*_{21}>*a*_{21}, the hole correction replaces *p*(*Λ*_{2} | *Λ*_{1}) with its expected value in *Λ*_{1}:
*R*_{21}≤*a*_{21} we would set both *p*(*Λ*_{2} | *Λ*_{1})=0 and *χ*=0 is because we are interested in the limit for small *χ*→0 when *χ* to the effective wave is smaller than

### (b) Infinitely many cylinders in the halfspace

In preceding sections, we considered a finite number of scatterers in a bounded domain *x*>0. We will follow [5] and limit the cylinders to the halfspace *x*>0, as it allows us to avoid divergent integrals, such as those in [46], e.g. between their eqn (32) and (33).

Substituting the approximations (3.18) and (3.11) into the governing system (3.6) leads to
*p*(**s**_{2})=*P*(**S**_{2}=**s**_{2}). By taking the limits

Incidentally, when all cylinders are identical this system reduces to eqn (54) in [5], that is when *p*(**s**_{2})=*δ*(*a*_{2}−*a*)*δ*(*c*_{2}−*c*)*δ*(*ρ*_{2}−*ρ*) in (3.23), where *δ*(*x*) represents Dirac’s delta function.

## 4. Effective wavenumber

To solve the system (3.24), first we use the symmetry of the problem to rewrite
*y*′, then taking *y*′=*y* we see that (4.1) is also a solution, recalling that *I*_{1}=*e*^{iαx+iβy} and
*i*^{m} *e*^{−imθ*} is introduced for later convenience. We could have for generality considered a sum of plane waves, but for low number density this is unnecessary, as we would find a unique **k**_{*} for a halfspace.

Equating (4.1) and (4.3), for *θ*_{*} and *k*_{*} are complex numbers. We also require that Im *α*_{*}>0, so that the integral over **x**_{2} in (3.24) converges.

In appendix B, we present the derivation for the system below, which is obtained by substituting (4.1) and (4.3) into (3.24). In the process we establish that *k*_{*}≠*k*, and find that there is no restriction on the length *k*_{*}, where
*hole correction* (3.17). For a more general pair distribution (3.19), we obtain

To determine the effective wavenumber *k*_{*} we need only use (4.5). That is, the solution *k*_{*} is the one that leads to non-trivial solutions for the function *c*. To completely determine

Next, we determine closed-form estimates for *k*_{*} from (4.5) and determine the corresponding coefficients

### (a) Explicit expressions for *k*_{*} via expansions in the number density

We now consider the expansions
*K*_{*0}=*k*^{2} , though we deduce this in §7. Here we have
*f*_{°} corresponds to the far-field scattering pattern for a single circular cylinder. This is evident by taking *N*=1 in (2.11), which leads to
*k*^{2} corresponds to the incident wave, the second

We can further specialize the wavenumber (4.11) by considering wavelengths 2*π*/*k* larger than the largest cylinder radius, or more precisely
*f*_{°}〉′(0)=0. This expression and the derivation is analogous to that given in ([5], eqn (86)) for a single species. Although the sums in (4.11) converge quickly for small *ka*_{*}, the form (4.17) is convenient as it is written in terms of the far-field scattering pattern 〈 *f*_{°}〉.

An alternative approach [47] that is very useful in the context of low-frequency propagation is to take the quasi-static limit of the system (4.5). For small *ka*, the monopole and dipole scattering coefficients are both *O*((*ka*)^{2}), which are the only contributions to the effective bulk modulus and density, respectively. Following the approach in [47], it is straightforward to show that for the *N*-species case, where the *n*th species has volume fraction *ϕ*_{n}, bulk modulus *K*_{n} and density *ρ*_{n}, the effective bulk modulus *K*_{*} and density *ρ*_{*} take the form
*D*_{n}=(*ρ*_{n}−*ρ*)/(*ρ*_{n}+*ρ*).

Next, we explore how the expression (4.11) compares with other approaches by evaluating it numerically. In §7, we develop analytical expressions for the average scattering coefficient

## 5. Two species of cylinders

In this section, we analytically compare two approaches to calculating the effective wavenumber of a multi-species material. The first self-consistent type method homogenizes the small cylinder distribution and then determines effective properties for a large cylinder distribution embedded in the homogenized background, as shown in figure 2. The second determines the multi-species result using the approach outlined in previous sections.

### (a) One small and one large species

We begin by assuming that there are only two species, *S* and *L*, that have constant wave speeds *c*_{S} and *c*_{L}, densities *ρ*_{S} and *ρ*_{L}, and number fractions *ϕ*_{S}∝*ϕ*_{L}, so we will discard *ϕ*=*ϕ*_{S}+*ϕ*_{L} denotes the total volume filling fraction. Note that it is more precise to assume small *ϕ*, rather than a small number density, since *ϕ* is non-dimensional.

First, the effective wavenumber *k*_{*S} of a material at long wavelengths with only a single species of *S*-cylinders is obtained by simplifying (4.17), where the far-field pattern (4.12) is, therefore, just the *S* cylinder species, i.e. there is no integral over **s**_{1} and *p*(**s**_{1})=1. Assuming *c*_{S}≠0 and *ρ*_{S}≠0, we use ([8], eqn (24)) for small cylinder radius, which in our notation (recall *N*=1, given by

Next, we determine the effective wavenumber for large scatterers embedded in a background described by *k*_{*S} and *ρ*_{*S}. For this step, we introduce the notation *f*_{°}(0,**s**_{1})=*f*_{°}(0,**s**_{1},*ρ*,*k*), which expresses the problem in terms of density and wavenumber in place of density and wave speed. Consequently, from (5.4), we have
*f*_{°L}(0):=*f*_{°}(0,**s**_{L},*ρ*,*k*).

To calculate the wavenumber *k*_{*LS} for the *L*-cylinders in a material with a wavenumber *k*_{*S}, we use the formula (4.11) with *k* replaced by *k*_{*S}, 〈 *f*_{°}〉 replaced with *f*_{°}(0,**s**_{L},*ρ**_{S},*k*_{*S}) above and keeping only the integrands, that is, removing the integrals over the multi-species **s**_{1} and **s**_{2}, to arrive at
*O*(*ϕ*_{L}*ϕ*_{S}) in the above does not agree with (4.11), even in the limit *a*_{S}→0, as we show next.

For only two species of cylinders, and assuming the cylinders are uniformly distributed, the probability density function for the scattering properties becomes
*a*_{LS}=*a*_{SL}. Assuming that *a*_{S}≪1 and using *a*_{LS}=*b*_{S}+*b*_{L}, with *b*_{S}≥*a*_{S} and *b*_{L}≥*a*_{L}, we expand *d*_{m} (4.14) as
*b*_{S}∝*a*_{S}. Substituting the above, (5.1), (5.3) and (5.8), into (5.10) we obtain
*Z*^{p}(**s**_{S})/*a*^{2}_{S} converges when *a*_{S}→0, see (5.1).

The terms in the brackets in (5.12) account for the interaction between the two types of cylinders, which is where the wavenumbers *k*_{*LS} (5.8) and *k*_{*} (5.12) differ. The leading-order error is non-vanishing even as the radius of the small species vanishes and is given by
*δf*_{LS} is the change in the far-field scattering pattern of the *L*-cylinders due to changing the background wavenumber from *k* to *k*_{*S}, while *G*_{0} accounts for the multiple scattering between the *L* and *S*-cylinders, which becomes significant when both *ϕ*_{S} are *ϕ*_{L} are large. Ultimately, this means that *k*_{*}≈*k*_{*LS} if either the *S* or *L*-cylinders are very dilute.

For comparison, we also give the two-dimensional version of eqn (23) of [28] given by

## 6. Numerical examples

In this section, we consider a selection of numerical examples to demonstrate the efficacy of (5.10) and other expressions. For *k*_{*LS}, we use the exact formula (4.11) for one species and then equation (5.7). This way, *k*_{*LS} is valid for *S*-cylinders with approximately zero density such as air. In the graphs that follow we use
*k*_{*} will be replaced with *k*_{*LS} and *k*_{*0} depending on the context.

For reference, we provide Julia [49] code to calculate the effective wavenumbers.

### (a) Two-dimensional emulsion

Here we consider an emulsion composed of hexadecane (oil) and glycerol in water [50], table 2, where the glycerol forms very small inclusions. The graphs of figure 3*a* show how *k*_{*}, *k*_{*LS} and *k*_{*C} differ when varying only *a*_{S} the radius of the glycerol inclusions, for a fixed angular frequency *ω*=*c*/*k*_{0}≈3×10^{6} Hz. We observe that the difference between *k*_{*} and *k*_{*LS} persists even as *a*_{S}→0, as expected according to (5.12). Meaning that, no matter how small the *S*-cylinders become, the larger cylinders *L* do not perceive the *S*-cylinders as a homogeneous material, in the naive way described in §5a.

In fig. 21 of [28], they observed that experimentally measured wave speeds were shifted in comparison to the *k*_{*C} predictions, even for low-frequency. We can see this same discrepancy in figure 3*b*, where the angular frequency is varied between 1 KHz<*ω*<12 MHz while the radius *a*_{S}=25 μ*m* is fixed. This discrepancy is due to the terms of order *k*_{*C} (5.15) and *k*_{*LS} (5.10). Although all three wavenumbers are similar in figure 3*b*. The same is not true when we increase the frequency.

In figure 3*c*, we show how *k*_{*C}, valid only for low-frequency, strays from the more accurate *k*_{*} as the frequency increases,^{1} where we did not include *k*_{*LS} as it is only valid for low frequency. There we can see that *k*_{*C} performs well up to about *ka*_{S}=0.3, at which point *ka*_{L}=3.0.

All the approximations *k*_{*0} (5.16), *k*_{*LS} (5.10) and *k*_{*C} (5.15) are missing second-order terms in the number density. In figure 4, we see the effect of these missing terms by varying the volume fraction while fixing *ω*=3×10^{6}, or equivalently *ka*_{S}=0.5. In the limit of low volume fraction, all three effective wavenumbers agree, as expected. For the largest volume fraction 40%, the expected error^{2} of *k*_{*} is 6%. However, the relative difference between the attenuations of *k*_{*LS} and *k*_{*0} and the multi-species attenuation of *k*_{*} reaches 30%.

Summarizing figures 3 and 4, all the approximations are similar for either low frequency or low volume fraction. This is because the three phases in table 2 have similar properties. In our next example one of the phases, air, will be very different from the others, which will lead to more dramatic differences.

### (b) Two-dimensional concrete

When there is a high contrast in the properties of the inclusions, multiple scattering can have a dramatic effect. To demonstrate this we consider a concrete-like material made from a limestone possessing cylinders of brick and air, given in table 3.

Figure 5 shows that it is only in the low-frequency limit, *ka*_{S}<0.05, that the wavenumbers *k*_{*C} and *k*_{*LS} agree with the more exact *k*_{*}, which has a maximum expected relative error of only *ϕ*^{3}=0.16^{3}≈0.4%. And in figure 5*b*, the wavenumber *k*_{*LS} appears to hit a resonance which should not be present. This, and the dramatic changes in attenuation at low frequency, are expected because for an inclusion with low density, the effective wavenumber diverges for fixed volume fraction when *k* tends to zero [8]. Figure 5*c* shows the limitations of *k*_{*C} as the frequency increases. Even though *k*_{*C} is only valid for low frequencies, its results are quite close to *k*_{*}, having a relative difference of around 25%.

Again as expected, all the wavenumbers converge as the volume fraction decreases (figure 6), yet the differences in the wave speed are significant, reaching 100% in this example, when the total volume fraction *ϕ*=40%.

## 7. The average field and reflection

In this section, we determine the reflected field from a halfspace, which is achieved by deducing the averaged scattering coefficient _{2}.

In order to calculate *K*_{*0}=*k*^{2}. As *θ*_{*} appears in (4.6), we expand

For *K*_{*0}=*k*^{2}, to ensure that *K*_{*0}=*k*^{2} in (7.2) leads to
*θ*_{*2} is given for completeness. We use the above to expand
*n* and **s**_{2}, therefore, *K*_{*0}=*k*^{2} and *θ*_{*0}=*θ*_{in}, and the last equation is from (4.15).

To calculate the next order in *x*_{2}>0. That is, in the limit where there are no cylinders, except one fixed at **x**_{2}, the averaged scattering coefficient *F*_{*}, given by (B.20), is independent of *n* and **s**_{2}. So if we substitute *F*_{*}, we can then take *F*_{*} outside the sum and integral in (7.8), and then substitute back _{2}. The above reduces to the one-species case given in ([52], eqn (27)). With

### (a) Reflection from a halfspace

Here we calculate the reflected wave measured at (*x*,*y*), where *x*<0. To achieve this, we assume that the boundary layer around *x*=0 has little effect on the reflected wave, that is, we assume most of the scatterers behave as if they are in an infinite medium. This is the same as taking *et al.* [53] discuss the possibility of a boundary layer effect even in the low-frequency limit.

Substituting (3.9) into the total effective wave (3.8), and using the form (4.3) reveals
*u*_{in}(*x*,*y*)=*e*^{ik⋅x}, *x*_{1}−*x*>0. Using the above in (7.10), we reach
*θ*_{ref}=*π*−*θ*_{*}−*θ*_{in}. The reflected wave shown by (1.2) is calculated by expanding for small

## 8. Conclusion

We have deduced the effective wavenumbers (1.1) and (1.5), and reflection coefficient (1.2), for a multi-species material up to moderate number density and over a broad range of frequencies. This will enable experimental researchers to extract more information about the makeup of inhomogeneous media (see the electronic supplementary material for self-contained expressions for the wavenumbers and reflection coefficients in the case of a finite number of species). We also remark that the results may be extended straightforwardly to multiple scattering from cylinders in a number of contexts, including two-dimensional electromagnetism.

Characterization is not the only application; this theory can also be employed to design novel materials. We have shown that multiple scattering between different species can lead to effective properties that are not exhibited by single-species media. That is, using our multi-species formulae it is now possible to choose species so as to design impedance matched, highly dispersive and broad band attenuating materials.

We also saw that the multi-species effective wavenumbers derived in the acoustics literature were accurate for low frequency and low volume fraction. But to go beyond these limitations, our more precise effective wavenumber was needed. We also illustrated that a ‘self-consistent’ approach to calculating the effective wavenumber is not even accurate at low frequencies.

Two main issues of our method deserve further investigation: the effects of the boundary layer near the surface of the halfspace and the QCA. To calculate the reflection coefficient up to second order in the number density, we neglected the effects of the boundary layer. It is not clear how to theoretically make progress without these two approximations, nor what errors they introduce. We believe that these issues represent important future work.

## Data accessibility

We provide the code to generate all the graphs in [31]. We also provide electronic supplementary material with self-contained formulae.

## Authors' contributions

A.L.G. drafted the manuscript, carried out and verified all the calculations. A.L.G., W.J.P. and I.D.A. conceived of the study. M.J.A.S. contributed to §5 and W.J.P to §4. All authors edited the manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

A.L.G., W.J.P. and I.D.A. gratefully acknowledge funding provided by EPSRC (EP/M026205/1). W.J.P. and M.J.A.S. gratefully acknowledge funding provided by EPSRC (EP/L018039/1). I.D.A. undertook part of this work while in receipt of a Royal Society Wolfson Research Merit Award, and part was supported by the Isaac Newton Institute under EPSRC (EP/K032208/1).

## Acknowledgements

A.L.G. thank the Mathematics of Waves and Materials (MWM) research group at the University of Manchester for their support and discussions.

## Appendix A. Effective wavenumber for multi-species spherical inclusions

In this section, we apply our multi-species theory to the results in [21] for spherical inclusions to reach the effective wavenumber (1.5). Details are omitted when the results follow by direct analogy. For spheres we define the ensemble-average far-field pattern and multiple scattering pattern,
*q* takes the values
*D*_{m}(*x*)=*xj*_{m}′(*x*)(*xh*_{m}′(*x*)+*h*_{m}(*x*))+(*x*^{2}−*m*(*m*+1))*j*_{m}(*x*)*j*_{m}(*x*), *P*_{n} are Legendre polynomials, *j*_{m} are spherical Bessel functions, *h*_{m} are spherical Hankel functions of the first kind and
*q*_{j}=(*ρ*_{j}*c*_{j})/(*ρc*), where

## Appendix B. Calculating the effective equations (4.6) and (4.5)

In this appendix, we provide an in-depth outline of the derivation for (4.6) and (4.5), which are given in terms of the unknowns *k*_{*}. The approach here is similar to Section IV in [5].

We begin by substituting (4.1) and (4.3) into the integral in the governing equation (3.24) which for *X*=*x*_{2}−*x*_{1} and *Y* =*y*_{2}−*y*_{1}, so that (*R*,*Θ*) are the polar coordinates of **X**=(*X*,*Y*), where *R*=*R*_{21} and *Θ*=*Θ*_{21}−*π*, and we define

To calculate the last integral in (B 1), it is necessary that *k*_{*}≠*k*, as *k*_{*}=*k* leads to a divergent integrand. Assuming *k*_{*}≠*k*, we observe that
^{2}*Φ*_{n−m}(*kR*,*Θ*)=−*k*^{2}*Φ*_{n−m}(*kR*,*Θ*) and

where **n** is the outwards pointing unit normal. For *x*_{1}>*a*_{21}, the region *R*>*a*_{21}, so that the boundary *R*=*a*_{21} and the line *λ*↦*α*_{*}, and where we have included the case *X*>0 for future reference).

From (B 5), it follows that (B 1) now becomes
*I*_{1}=*e*^{i(α*x1+βy1)}. As the above must hold for all *x*_{1}, we can equate the terms in the brackets to zero, which leads to (4.6) and (4.5).

**(a) Effective wave for any pair distribution**

To calculate the effective wave for any pair distribution function *χ* satisfying (3.21), we substitute (3.19) into (3.24), which leads to an extra integral appearing on the left side of (B 1):
*X*=*x*_{2}−*x*_{1} and *Y* =*y*_{2}−*y*_{1}, so that (*R*,*Θ*) are the polar coordinates of **X**=(*X*,*Y*), where *R*=*R*_{21}, *Θ*=*Θ*_{21}−*π*. The second integral on the right is zero when **X** can not satisfy both

For

**(b) Low number fraction**

To calculate *k*^{2}_{*} we need only equation (4.5), which after substituting (4.10) and
*d*_{m} defined by (4.14), and then equating terms of order

Turning to (B 13), we see that *m* and **s**_{1}. Let

Turning to (B 14), we see that
*m* and **s**_{1}, which we use to write (B 14) as
*K*_{*2}=−4*i*〈 *f*_{°°}〉(0), which together with (B 15), (4.10) leads to the effective wavenumber (4.11).

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4052399.

↵1 In this case, we did not exactly use the two-dimensional version of eqn (23) from [28], but instead used a more accurate version where we summed enough terms for the far-field patterns to converge.

↵2 If we disregard the error due to the low-frequency assumptions, we can estimate the expected errors

*ϕ*^{2}=0.4^{2}=16% for*k*_{*0}and 2*ϕ*_{L}*ϕ*_{S}=2*0.2^{2}=8% for both*k*_{*LS}and*k*_{*C}.

- Received December 11, 2017.
- Accepted March 7, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.