Reflection from a multi-species material and its transmitted effective wavenumber.

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.

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 and background wavenumber k, our main results for a multi-species material comprising circular cylinders are the effective transmitted wavenumber: and for an incident wave uin=ei, with , the averaged reflected wave from the inhomogeneous halfspace, where and the functions 〈 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.

Figure 7.

The far-field reflected angle θref, where k=(α,β) and k* is the effective transmitted wavenumber defined in §4. The wavenumber k* results from ensemble averaging all scattered waves originating from x>0 (to the right of the dashed vertical line). The reflected field measured at x can be understood as the scattering (the grey circles) of the transmitted wave by an effective particle (grey particle). In the figure, θref equals π−θ*−θin, but for small number density , which is why θref=π−2θin appears in (1.3) and (1.4).

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 where 〈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.

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, the pressure u outside all the cylinders satisfies the scalar Helmholtz equation and inside the jth cylinder the pressure u satisfies where ∇2 is the two-dimensional Laplacian and We consider an incident plane wave and use for each cylinder the polar coordinates where x is the centre of the jth 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 as the scattered pressure field from the jth cylinder, where H are Hankel functions of the first kind, are arbitrary coefficients and characterizes the type of scatterer: with q=(ρc)/(ρc). In the limits q→0 or , the coefficients for Dirichlet or Neumann boundary conditions are recovered, respectively.

Table 1.

Summary of material properties and notation. The index j refers to properties of the jth species, see figure 1 for an illustration.

background properties:		density ρ	sound speed c
species properties:	number density n_j	density ρj	sound speed cj	radius aj

Figure 1.

Represents a multi-species material comprising different species of cylinders to the right of the origin O=(0,0). The vector x points to the centre of the jth cylinder, with a local polar coordinate system (R,Θ). Each cylinder has a radius a, density ρ and wave speed c, while the background has density ρ and wave speed c. The vector k is the direction of the incident plane wave. (Online version in colour.)

The pressure outside all cylinders is the sum of the incident wave uin and all scattered waves, and the total field inside the jth cylinder is In the above, J are Bessel functions of the first kind. The arbitrary constants and in (2.4) and (2.7) will be determined from the boundary conditions of the jth cylinder R=a. The boundary conditions of continuity of pressure and normal velocity on the cylinder boundaries are given, respectively, by recalling that ρ and ρ denote the material densities of the background and of the jth cylinder respectively. To impose the boundary conditions, we now express the relevant fields in terms of the (R,Θ) coordinate system. For the incident wave

where I=uin(x,y) following the Jacobi–Anger expansion [37]. For the scattered waves (2.4), we use Graf’s addition theorem (9.1.79) in [38], where (Rℓ,Θℓ) is the polar form of the vector x−xℓ. Using the above and (2.9) we can impose the boundary conditions (2.8) to arrive at the following system of equations for j=1,…,N and all integers m. Furthermore, the coefficients associated with the pressure inside the cylinder (2.7) are then given by and subsequently, the field 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 , denoted by , from which we can calculate an effective wavenumber and reflection. Note that the 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.

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 x1,x2,…,x with the scattering properties s1,s2,…,s, where s denotes the properties of the jth cylinder, i.e. here these are s=(a,ρ,c). Each x is in the region , where is the total number density and is the area of . The properties s are taken from the set . For example, we could have , so that a∈[0,1], ρ∈[1,2] and c∈[100,200].

The probability of the cylinders being in a specific configuration is given by the probability density function p({x1,s1},{x2,s2},…,{x,s}). Using the compact notation ={x,s} to denote the properties of the i-th cylinder, it follows that where each integral is taken over both (for x) and (for s) with d=dx ds. 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 and where p(1,…, | ) is the conditional probability of having cylinders with the properties 1,…, (not including ), given that the jth cylinder has the properties . Likewise, p(1,…, | ℓ,) is the conditional probability of having cylinders with the properties 1,…, (not including ℓ and ) given that there are already two cylinders present, with properties ℓ and .

Given some function F(1,…,), we denote its average, or expected value, by If we fix the location and properties of the jth cylinder, and average over all the properties of the other cylinders, we obtain a conditional average of F given by where we do not integrate over . The average and conditional averages are related by where 〈F〉 is the conditional average when fixing both and ℓ, and 〈F〉=〈F〉.

Returning to the task of obtaining effective properties for a multi-species medium, we multiply the system (2.11) by p(2,…,|1) and average over 2,…,, to reach where, without loss of generality, we have chosen j=1, used the conditional average definition (3.2b) and defined to make the dependency on 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=, then , because equation (2.11) uniquely determines the coefficients from the position and scattering properties ℓ. We use this to obtain Our aim is to solve the system above for ; however, this requires that we make assumptions about p(2 | 1) and . These approximations are discussed in §3a; however, for the moment, we assume that an appropriate substitution has been imposed.

With , we can calculate the average total pressure (incident plus scattered), measured at some position x outside , by averaging (2.6) to obtain where 〈uin(x,y)〉=uin(x,y) because the incident field is independent of the scattering configuration. We can then rewrite the average outgoing wave u by fixing the properties of the jth cylinder and using equation (3.2a) to reach Likewise, for the conditionally averaged scattered field (2.4) measured at x, we obtain We use the above to calculate the reflection from a halfspace in §7a and to obtain (1.2). To proceed we need to solve the system (3.6) and, in line with existing approaches, we do this by making statistical approximations.

Statistical approximations

In order to solve (3.6) for , we need an approximation for and the pair distribution p(2 | 1). In this work, we adopt the standard closure approximation for single species, but extended to multi-species, the QCA [5,6]: This approximation still makes sense for multi-species because it replaces the dependence of in 1 by its expected value in 1. Note also that the expected difference in 2: for a large number of scatterers.

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 For clarity, we introduce the spatial random variables X1,X1,…,X and the scattering property random variables S1,S1,…,S, and write probability density functions in the form, e.g. In the first instance, we assume the random uniform distribution where P(S1=s1) is the probability density in that the particle will have scattering property s1. The above assumes that , i.e. that the position x1 of the cylinder is independent of the scattering property s1. This is not always the case, for example, depending on the size of the cylinder, some positions near the boundary of may be infeasible. However, these boundary effects are negligible when taking the limit .

For the remaining distribution in (3.12), we use followed by which assumes that the scattering properties are statistically independent. Next, we assume that the cylinders are equally likely to be anywhere but do not overlap (a hole correction correlation [20]), to write where R21:=∥x1−x2∥, a21=b1+b2 for some b1≥a1 and some b2≥a2, and b1 is the radius of exclusion around x1 which is usually chosen to be proportional to the radius a1. Note that when integrating (3.17) above in x1 and x2, we obtain for .

Ultimately, substituting (3.16) and (3.17) into (3.15) in tandem with (3.14) leads to the pair distribution where H(x) denotes the Heaviside function, under the assumption . In the next section, we will the approximations (3.11) and (3.18) to solve the system in (3.6) for .

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 where To calculate the effective wavenumber for the pair-correlation (3.19), a common choice is to assume that the scatterers are uniformly randomly distributed, which leads to used by Linton & Martin [5], Section D, [42], [43], eqn (27), where is some distance large enough for the scatterers at x1 and x2 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 R21>a21, the hole correction replaces p(2 | 1) with its expected value in 1: just as QCA (3.10) assumes that . Similarly for R21≤a21 we would set both p(2 | 1)=0 and for QCA and hole correction. Another reason to set χ=0 is because we are interested in the limit for small . In this limit, it is expected that χ→0 when for uniformly distributed scatterers, which in turn indicates that the contribution of χ to the effective wave is smaller than  [5,45].

Infinitely many cylinders in the halfspace

In preceding sections, we considered a finite number of scatterers in a bounded domain . Now we consider the limit and where the region tends to the halfspace 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 where for brevity, we write with p(s2)=P(S2=s2). By taking the limits and , while fixing the number density , equation (3.22) takes the form which represents the governing system for our semi-infinite multiple-species problem.

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

Effective wavenumber

To solve the system (3.24), first we use the symmetry of the problem to rewrite that is, if is a solution to (3.24), then so is defined by for every y′, then taking y′=y we see that (4.1) is also a solution, recalling that I1=ei and Sufficiently far away from the boundary, say , we assume a plane wave ansatz where the factor i e−i 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 , we obtain Snell’s Law noting that both θ* and k* are complex numbers. We also require that Im α*>0, so that the integral over x2 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 , a fact we use to calculate the reflected wave. The result is that (3.24) reduces to the system and in terms of the unknown parameters and k*, where and , as we have assumed hole correction (3.17). For a more general pair distribution (3.19), we obtain where further details may be found in appendix Ba. We also remark that equation (4.5) reduces to ([5], eqn (87),[46], eqn (33)) for a single particle species.

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 . On the other hand, if is a solution to (4.5), then so is for any constant c. To completely determine we need to use (4.5) and (4.6).

Next, we determine closed-form estimates for k* from (4.5) and determine the corresponding coefficients for low number density in §7.

Explicit expressions for k* via expansions in the number density

We now consider the expansions where we use ∼ to denote an asymptotic expansion in , which is formally equivalent to an expansion in volume fraction. We show in appendix Bb how substituting the above into equation (4.5) leads to where we assumed K*0=k2 , though we deduce this in §7. Here we have and we introduce the multiple scattering pattern where and for convenience we define We remark that 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 and from (2.4) gives We can interpret the terms on the right-hand side of (4.11) in the following way: the first k2 corresponds to the incident wave, the second is the contribution from the incident wave scattered once from every cylinder (so-called ‘single scattering’) and the last is the contribution of this scattered wave being re-scattered by every cylinder (so-called ‘multiple scattering’).

We can further specialize the wavenumber (4.11) by considering wavelengths 2π/k larger than the largest cylinder radius, or more precisely which leads to where the integral converges because 〈 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 nth species has volume fraction ϕ, bulk modulus K and density ρ, the effective bulk modulus K* and density ρ* take the form where and D=(ρ−ρ)/(ρ+ρ).

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 and expressions for the reflection coefficient from the inhomogeneous halfspace.

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.

Figure 2.

(a) Two-species material comprising large (blue) and small (green) inclusions in a background material (yellow) with incident wavenumber k, (b) one-species material comprising the small inclusions alone, and (c) one-species material with large cylinders alone in a background with incident wavenumber k*, which is the effective wavenumber of (b). (Online version in colour.)

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 and c, densities ρ and ρ, and number fractions and , respectively. We assume that both types of cylinders have low volume fractions and and are proportional to one another ϕ∝ϕ, so we will discard terms, where ϕ=ϕ+ϕ 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* 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 s1 and p(s1)=1. Assuming c≠0 and ρ≠0, we use ([8], eqn (24)) for small cylinder radius, which in our notation (recall , where is given in (2.5)) where Substituting the above into the simplified (4.17) leads to after taking a Taylor series for small for the square root. We also calculate the effective density ([8], eqn (1)) or refer to (4.18) with N=1, given by which is appropriate for the approximation (5.3), see [8] for more details.

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

To calculate the wavenumber k* for the L-cylinders in a material with a wavenumber k*, we use the formula (4.11) with k replaced by k*, 〈 f°〉 replaced with f°(0,s,ρ*,k*) above and keeping only the integrands, that is, removing the integrals over the multi-species s1 and s2, to arrive at where we used . The above is an attempt to calculate the multi-species wavenumber by using only the single-species formula. However, the term of order O(ϕϕ) in the above does not agree with (4.11), even in the limit a→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 which when substituted into (4.11) leads to the multi-species result where we used a=a. Assuming that a≪1 and using a=b+b, with b≥a and b≥a, we expand d (4.14) as where we use b∝a. Substituting the above, (5.1), (5.3) and (5.8), into (5.10) we obtain where and Note that Z(s)/a2 converges when a→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* (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 A numerical investigation of this limit is conducted in §6. The physical meaning of these two terms are quite different: δf is the change in the far-field scattering pattern of the L-cylinders due to changing the background wavenumber from k to k*, while G0 accounts for the multiple scattering between the L and S-cylinders, which becomes significant when both ϕ are ϕ are large. Ultimately, this means that k*≈k* 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 and another commonly used approximation [48] describing the effective wavenumber for two species. However, the expression (5.15) is missing the interaction between the and species (the term ) and is only valid for low frequency. Additionally, the estimate (5.16) ignores terms of the order .

Numerical examples

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

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

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 3a show how k*, k* and k* differ when varying only a the radius of the glycerol inclusions, for a fixed angular frequency ω=c/k0≈3×106 Hz. We observe that the difference between k* and k* persists even as a→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.

Table 2.

Material properties used to approximate an emulsion.

	density (kg m−3)	speed (m s−1)	radius (μm)	volume %
distilled water	ρ=998	c=1496	—	84%
hexadecane	ρL=773	cL=1338	aL=250	ϕL=11%
glycerol	ρS=1260	cS=1904	aS=25	ϕS=11%

Figure 3.

Comparison of sound speeds and attenuation using wavenumbers k* (5.15), k* (5.7), and k* (5.10) for the water and oil emulsion from table 2. Code to generate figure: https://github.com/arturgower/EffectiveWaves.jl/tree/master/examples/emulsion. (a) Fixed wavenumber k=k0 while changing the radius a, (b) fixed radius a while changing the wavenumber k and (c) fixed radius a while changing the wavenumber k. (Online version in colour.)

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

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

All the approximations k*0 (5.16), k* (5.10) and k* (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×106, or equivalently ka=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* and k*0 and the multi-species attenuation of k* reaches 30%.

Figure 4.

Comparison of sound speeds and attenuation calculated from the effective wavenumbers k* (5.15), k*0 (5.16), k* (5.7) and the more accurate k* (5.10), as the total volume fraction of the inclusions increases (for the emulsion shown in table 2 with ka=0.5). Code: https://github.com/arturgower/EffectiveWaves.jl/tree/master/examples/emulsion. (Online version in colour.)

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.

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.

Table 3.

Material properties for our concrete-like material. Note that we used compacted limestone with very low porosity [51].

	density (kg m−3)	speed (m s−1)	radius (mm)	volume %
limestone	ρ=2460	c=4855	—	84%
brick	ρL=1800	cL=3650	aL=2.0	ϕL=10%
dry air	ρS=1.293	cS=331.4	aS=0.2	ϕS=6%

Figure 5 shows that it is only in the low-frequency limit, ka<0.05, that the wavenumbers k* and k* agree with the more exact k*, which has a maximum expected relative error of only ϕ3=0.163≈0.4%. And in figure 5b, the wavenumber k* 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 5c shows the limitations of k* as the frequency increases. Even though k* is only valid for low frequencies, its results are quite close to k*, having a relative difference of around 25%.

Figure 5.

Sound speed and attenuation from the approximate wavenumbers k* (5.15) and k* (5.7) with the more accurate k* (5.10) for the concrete-like mixture shown in table 3. Code: https://github.com/arturgower/EffectiveWaves.jl/tree/master/examples/concrete. (a) Fixed wavenumber k=k0 while changing the radius a, (b) same as (a) but for smaller ak and (c) fixed radius a while changing the wavenumber k. (Online version in colour.)

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%.

Figure 6.

Sound speed and attenuation from the three effective wavenumbers k* (5.15), k*0 (5.16) and k* (5.7), with the more accurate k* (5.10), against the total volume fraction of the inclusions for the concrete shown in table 3. (Online version in colour.)

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 , for low number density shown in (4.10)2.

In order to calculate we first use (4.6) and (4.11), except here we can deduce that K*0=k2. As θ* appears in (4.6), we expand which combined with Snell’s equation (4.4) and the number density expansions (4.10) gives for the first two orders:

For , which appears in equation (4.6), we assume that, as , is a fixed width large enough for the effective wave ansatz (4.3) to hold, meaning that Using the above in (4.6) we conclude that K*0=k2, to ensure that appears on the left-hand side. Subsequently, using K*0=k2 in (7.2) leads to and where θ*2 is given for completeness. We use the above to expand then we substitute the leading order term in the above into (4.6) leading to: . However, from (B.17) we found that is independent of n and s2, therefore, . This means that tends, in the limit , to the scattering coefficient of one lone cylinder: where we used K*0=k2 and θ*0=θin, and the last equation is from (4.15).

To calculate the next order in of equation (4.6) we need to calculate . To do so, we assume that as for every x2>0. That is, in the limit where there are no cylinders, except one fixed at x2, the averaged scattering coefficient tends to the scattering coefficient of one lone cylinder, even for . As a result which when substituted in from (4.6), together with (4.12), leads to . Substituting this, (7.6), and (7.5) into (4.6), and then ignoring second-order terms we obtain where we also used which is a result of the property , see (2.5), implying that . In appendix Bb, we showed that the quantity F*, given by (B.20), is independent of n and s2. So if we substitute for F*, we can then take F* outside the sum and integral in (7.8), and then substitute back to arrive at where we used (4.12)2. The above reduces to the one-species case given in ([52], eqn (27)). With and , we can now calculate reflection from a halfspace.

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 , which was also used in [52], where they showed that this approach matches other homogenization results in the low-frequency limit. We note, however, that Felbacq 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 where we used uin(x,y)=e, , , substituted (3.14), used and took the limit . Then using (B.9) and (B.11), we obtain noting that x1−x>0. Using the above in (7.10), we reach where θref=π−θ*−θin. The reflected wave shown by (1.2) is calculated by expanding for small , including , and then substituting the results from §7. For a single-species this formula reduces to ([52], eqn (41) and (42)). Figure 7 gives a pictorial representation of the reflection coefficient in (7.12).

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.