Effective dynamic constitutive parameters of acoustic metamaterials with random microstructure

A multiple scattering analysis in a non-viscous fluid is developed in order to predict the effective constitutive parameters of certain suspensions of disordered particles or bubbles. The analysis is based on an effective field approach, and uses suitable pair-correlation functions to account for the essential features of densely distributed particles. The effective medium that is equivalent to the original suspension of particles is a medium with space and time dispersion, and hence, its parameters are functions of the frequency of the incident acoustic wave. Under the quasi-crystalline approximation, novel expressions are presented for effective constitutive parameters, which are valid at any frequency and wavelength. The emerging possibility of designing fluid–particle mixtures to form acoustic metamaterials is discussed. Our theory provides a convenient tool for testing ideas in silico in the search for new metamaterials with specific desired properties. An important conclusion of the proposed approach is that negative constitutive parameters can also be achieved by using suspensions of particles with random microstructures with properties similar to those shown in periodic arrays of microstructures.


Introduction
The recent interest in effective material properties, especially for the effective mass density ρ eff , comes from the study of metamaterials [1]. Metamaterials are artificially fabricated structures (often periodic) which are designed so that they display unusual macroscopic properties. For example, various engineered materials with negative effective density and/or negative bulk modulus have been demonstrated [2][3][4]. Applications of these metamaterials include improved acoustic imaging [5,6] and sound wave manipulation for cloaking [7,8]. Most of these approaches rely on resonant inclusions and the resulting acoustic metamaterial parameters vary strongly with frequency. The prospect of modelling effective material properties of composites with microstructure brings about the need for fast and accurate computational schemes to test ideas in silico.
Here, our main concern is the dynamic acoustic properties of a mixture of spherical particles suspended in a non-viscous fluid. In order to estimate the effective mass density of fluid-particle mixtures, various different approaches are possible. The simplest starting point is the Archimedes principle in which the effective mass density of a two-phase composite (fluid matrix of density ρ 0 and inclusions of density ρ) is given by the volume average as 3 where φ is the volume fraction occupied by the inclusions. A similar formula can be used to calculate the effective bulk modulus κ eff . The inherent assumption in definition (1) is that ρ eff is insensitive to the relative motion of the phases of the composite. In some circumstances however, a different estimate for ρ eff is obtained. The reason is that regions of different phases do not necessarily move in unison, even in the long-wavelength limit, and therefore the local average momentum is not simply the product of the local averaged mass and sound speed. Ament [10] presented another expression, known as the 'Ament estimate', for the effective mass density that differs significantly from the intuitive mixture law. In the absence of viscosity, his estimate becomes Fikioris and Waterman [11] have provided an independent analytical verification of equation (2) in the long-wavelength limit using a multiple-scattering formulation. The same result was also obtained later by Kuster and Toksöz [12] and Gaunaurd and Wertman [13] using different methodologies. The effective mass density would satisfy the mixture law (1) if the matrix were an elastic medium. The expression (2) results from the inertia effects [14] which occur when the matrix is a fluid: the fluid can slide past the particle and induce an effect of added mass. Sound speed estimates based on equation (2) have also been shown to agree with experiments [15]. For small φD, formula (2) is approximated as Berryman [17] considered a 'self-consistent' version of equation (2) resulting in a quadratic equation for ρ eff , which can be compared with equation (2). For small φD, equation (4) is approximated as ρ eff ρ 0 1 + 3φD + 6φ 2 D 2 3ρ ρ 0 + 2ρ which differs from equation (3) in the quadratic order. It is of interest to extend the Ament estimate (2) to finite frequencies and wavelength. In this paper, a comprehensive description of acoustic scattering by multiple spheres using a multiple-scattering formulation is presented. It is based on the approach developed by Fikioris and Waterman [11]. This formulation describes a mixture of spherical particles suspended in a non-viscous fluid as a linearly viscoelastic material, whose elastic parameters with respect to the coherent acoustic wave motion are spatially constant but frequency dispersive and complex valued. Based on this approach (see section 5), the effective mass density depends on the volume fraction φ and on the scattering dispersion parameterω = ka, where k is the wavenumber in the fluid matrix and a is the particle radius. This dependence is later shown in section 5 to be as follows: where i is the imaginary unit, and r (ω) is a certain complex-valued function, which tends to zero along with its first derivatives atω = 0. The real part of equation (6) is recognized as the Ament estimate (2).

4
Estimating the effective material properties of inhomogeneous composites has been a subject of scientific and engineering interest with a long history. So far, many approaches and predictive schemes have been proposed. Early works were limited to either static or very-longwavelength approximations, e.g. [10,12,17]. More recently, frequency-dependent descriptions have emerged, e.g. [22][23][24][25][26]28]. Some schemes include the self-consistent and effective medium methods for which a substantial body of literature may be found. Although variants exist, these methods often consider the scattering problem of a coated particle in an effective medium. The coated particle consists of a layered particle of inner radius a and a concentric layer of background material with an outer radius b. The actual boundary condition of the problem is imposed on r = a; within the coating, the wavenumber is k; outside the coating (r > b), the wavenumber is K (which is unknown). Transmission conditions are imposed across r = b, and the exciting field is a plane wave travelling with the wavenumber K. This leads to what Twersky [18] called the 'schizoid' or 'two-space' problem (see also [19]). Given b and K the coated-particle problem can be solved exactly. Common choices for the radius b are b = a (no coating) or such that (a/b) 3 = φ, the volume fraction occupied by the particles. The effective medium parameters are determined by applying the condition that total scattering vanishes in the limit when Kb 1. This approach allows determination of the effective mass density; however, the result is formally restricted to the low-frequency and long-wave ranges (see, e.g., [20] and references therein). Using only the leading-order terms, this estimate for the effective mass density is reduced to equation (6), as is shown in the remainder of the paper. It is necessary to note that a deficiency of all current versions of the effective medium methods is their failure to describe the influence of the spatial distribution of particles on the effective constitutive parameters. Such a description is possible in the framework of a self-consistent scheme called the effective field method (see, e.g., [37]) and our work is within the framework of this scheme.
The effective material properties of periodic acoustic structures can be easily calculated in the long-wavelength limit (see, e.g., [21,69]). The periodicity means that the equation of motion can be solved exactly. For random microstructures however, this equation can be solved exactly only for a small number of particles. In order to simplify the problem of solving the equation of motion in the presence of multiple particles, an approximation which neglects the finite size of the particle is customary. Several authors, for example, have treated a random suspension as a mixture of particles distributed in a fluid and the effective wavenumber K has been calculated explicitly [9,22,26,27]. This approach yields a simple formula for the effective mass density ρ eff , which is linear in volume fraction [28]. At higher concentrations, however, this formula is unjustified and leads to incorrect results. This point is further discussed in the present paper.
Acoustic metamaterials can be made by setting up a periodic distribution of resonant particles in a fluid [29]. When the particles are arranged in structures having a square or hexagonal symmetry, the resulting system effectively behaves like a homogeneous and isotropic fluid. This behaviour is obtained when the wavelength of the propagating acoustic wave is larger than about a few times the lattice separation between individual scatterers. Materials behaving with an effective anisotropic mass density can also be engineered by arranging particles on nonsymmetric lattices [30]. Acoustic parameters depending on the local coordinates can be achieved by simply changing the dimensions of the particles [31]. We show in the present paper that acoustic metamaterials can also be engineered by using suspensions of particles with random microstructures.
The paper proceeds as follows. The problem is formulated in section 2 where the multiplescattering theory is outlined. In section 3, an implicit form of the wavenumber equation is obtained using matrix notation. Although closed-form asymptotics of the effective wavenumber can be derived with the desired accuracy [33], the remainder of section 3 corroborates results for the effective wavenumber accurate to second order in volume fraction φ for finite-size or point-like particles. The main results of the paper are then derived in sections 4 and 5. Based on the above-mentioned matrix formulation and the asymptotic approximations of section 3, the reflection coefficient for the random array of spheres is obtained in section 4. Using the ansatz that the sought effective medium is described acoustically by the effective wavenumber and the reflection coefficient for the associated semi-infinite suspension of particles, the corresponding effective material properties are derived in section 5. A significant outcome of the proposed approach is that explicit expressions are obtained for the effective material parameters. Although the statement of the problem was completed by employing a specific pair-correlation function, the entire analysis is generalized to arbitrary pair-correlation functions in section 6. The possibility of designing random fluid-particle mixtures to form acoustic metamaterials is discussed in section 7. Low-frequency, long-wavelength expansions of the effective parameters are also derived in the appendices.
We consider a semi-infinite suspension of particles in a non-viscous fluid of density ρ 0 , bulk modulus κ 0 and sound speed c 0 . The spheres are of identical geometry and elastic properties and uniformly and randomly distributed in the half-space z > 0, where their number density is n 0 .
A plane wave is incident on the half-space and propagates in the z-direction, as shown in figure 1. Since viscosity is neglected, only pressure waves can propagate in the fluid. These waves cause the fluid particles to move in the z-direction. The dynamics of the current multiplescattering problem can be expressed in terms of appropriate scalar potentials. If the time 6 dependence is of the form exp(−iωt), the displacement potential φ in is given by where φ 0 is an amplitude factor, ω is the angular frequency, and the time dependence has been omitted for brevity. We assume that the scattering properties of an isolated sphere are fully described by a set of coefficients T n , which are known. Therefore, the transition operator T for every sphere is assumed to have translational invariance with where ψ n (r) = h n (kr )P n (cos θ ) is an outgoing spherical wavefunction, h n ≡ h (1) n is a spherical Hankel function, and P n is a Legendre polynomial. Observe in equation (8) that there is no dependence on the azimuthal angle ϕ, owing to the axial symmetry of the sphere.
The far-field scattering amplitude of an isolated sphere is then given by The angular shape function f (θ ) is an infinite series with coefficients equal to the scattering amplitudes T n , i.e.
Note that f (0) and f (π ) are measures, in the forward and backward directions, of the displacement potential due to scattering by a single sphere.
The incident wave (7) is subjected to multiple scattering between the spheres. For any given configuration (ensemble) of a finite number of spheres, the scattered motion may be computed exactly by solving the appropriate system of equations. However, for the case in hand the distribution of spheres is known in a probabilistic sense only. Thus, it is common to regard the volume containing the spheres as a random medium with certain average (homogenized) properties. We perform ensemble averaging [19,22,34,35] in order to calculate the average (coherent) motion. There exist a number of approaches to solving such problems; a survey can be found, for example, in [36]. The method of effective fields [37] is used in our work. This allows us to take into account the spatial distributions of particles via specific pair-correlation functions g(ρ ji ). For simplicity, the spatial distribution of particles is assumed homogeneous and isotropic, and the 'hole correction' (that prevents the overlapping of the particles in the averaging process) [11] is adopted in the following, for which, where H is the Heaviside unit function: H (x) = 1 for x > 0 and H (x) = 0 for x < 0. The hole radius b satisfies b 2a so that spheres are not allowed to overlap. More complex forms of pair-correlation functions are considered in [33]. Appendix A outlines the relevant multiplescattering theory.
Subject to the special case of the hole-correction (11), the subsequent analysis provides a comprehensive description of the coherent wave propagation in fluid-particle mixtures in terms of two distinct scattering processes: single and multiple scattering. A generalization of the main results to arbitrary pair-correlation functions is also provided in section 6.

The wavenumber matrix and the associated eigenvector
The quasi-crystalline approximation [32] allows us to develop the dispersion equation for the effective wavenumber K of the coherent acoustic field in the fluid-particle mixture. This equation follows from equation (A.7a) and can be written in matrix form [33] where e = (1, 1, . . .) t is a constant vector, I is the identity matrix and the scalars ε and ξ are The infinite square matrices R and Q have elements where and In equation (14), G is a Gaunt coefficient [38] and the summation over p is finite, covering the range |n − ν| to (n + ν) in steps of 2, so that ( p + n + ν) is even.
A solution x = 0 can only exist for the particular value of ξ that makes the matrix premultiplying x in equation (12) singular. We then obtain an implicit equation for the effective wavenumber K, where the effective scattering amplitude F is given by [33] Equation (18) is an exact expression for the effective wavenumber. (The only assumption made to arrive at this result is the quasi-crystalline approximation). Note that the formula (19) separates the implicit form of the effective wavenumber into two distinct parts. One part is defined by the scattering matrix Q, which describes the response of a single particle to a plane incident harmonic wave. The other part, R, is defined by the spatial arrangements of the particles and accounts for multiple scattering. Note that for a regular array of particles, the quasi-crystalline approximation is exact, in which case the multiple-scattering matrix R can be reduced to a known lattice sum [39].
The infinite eigenvector x of amplitudes X n , associated with the wavenumber equation (18), follows from equation (A.7b) of appendix A and equation (12) as 8 The eigenvector x will be useful later in the analysis for calculating the reflection coefficient for the coherent wave reflected at the interface z = 0.
Equations (18) and (20) are the most convenient starting point for the calculation of the effective parameters of the fluid-particle mixture. We will now study some important specific cases to render the wavenumber equation (18) more tractable.

Asymptotic expansion
At low concentration (εa 3 1), we assume an asymptotic expansion in powers of ε as follows: Following the procedure detailed by Caleap et al [33], the coefficients in the expansion (21) are obtained as where R 0 = R(0). The elements of the matrix R 0 are then found by expanding the function R p (Kb) in equation (16) for small (Kb − kb). One finds that where In terms of the Fourier series, the coefficients in equations (22) may be written as The results in classical multiple scattering theories are usually defined in terms of the angular shape function f (θ ) for the scattering of a plane wave by a single particle. By using the definition of the Gaunt coefficients (A.10) and the angular shape function (10), equations (22) (or (25)) may be expressed as where and It is useful to describe yet another form of the coefficients (26). This is done in the next section, which considers point-particle approximation: in addition to εa 3 1, we also require kb 1.

Point particles: the small kb limit
The leading order term in kb of equation (A.8) is M p (Kb) (K/k) p . Then, it follows from equation (16) that By means of these approximations, equation (27) simplifies. Then, by using the generating function for the Legendre polynomials, we obtain where The result in equation (31) is accurate to second order in terms of the concentration expansion term. Note that the first two terms of equation (31) give the thin slab result of Fermi [42] who considered the propagation of neutrons through a thin sheet of scatterers. Lloyd and Berry [43] have applied the multiple scattering treatment of Lloyd [44,45] (developed for the treatment of scattering in electronic structures of liquids and disordered alloys) and obtained precisely the result in equation (31). Equation (31) defines what is called the 'Lloyd-Berry estimate' in the remainder.

The reflection coefficient
The interaction of the incident wave (7) with the scatterers causes it to be reflected. The reflection coefficient represents the average back-scattered amplitude in the domain {z < 0}. As shown in [11], equation (A.1) (with n(r j ) = n 0 ) yields for the reflected field where the reflection coefficient R is defined as As expected, R depends on the infinite eigenvector of amplitudes X n .
Using the notation from section 3, equation (A.7b) may be written as Once the eigenvector x is determined, the effective forward-and back-scattered amplitudes, F 0 and F π , follow from equations (A.7b) and (34), respectively, which can be recast in matrix form as where the additional infinite matrix J has elements J nν = 0, J νν = cos νπ . Then, the reflection coefficient is obtained from Equation (37) is an exact expression for the reflection coefficient. As in section 3.1, we consider next the low-concentration asymptotic approximation. We seek an asymptotic expansion of R in powers of the parameter ε: Using the form (20) in equations (36) and (37), the individual terms in (38) may be obtained as Subsequent terms become more complicated but the procedure for finding them is straightforward. An alternative expression to equation (38) can be obtained in terms of Fourier series, i.e.
By using the definition of the Gaunt coefficients (A.10) and the angular shape function (10), equation (40) may also be expressed where the function S is defined in equation (27), and To conclude this section, we give the point-particle approximation of the low-concentration expansion (41). In addition to a 3 1, we also require kb 1. We use the same procedure as that described in section 3.2. In terms of n 0 , we obtain Based on the knowledge of the effective wavenumber and the reflection coefficient of the associated semi-infinite suspension of particles in a fluid, we calculate next the effective dynamic properties of the fluid-particle mixture.

Effective dynamic constitutive parameters
In this section, we extend the effective mass density expressions (2) and (3) to finite frequencies and wavelength. Specifically, we show that a frequency-dependent effective mass density can be obtained from the reflection coefficient of the previous section, which reduces to equations (2) and (3) in the long-wavelength limit and which can yield complex values. The bulk modulus κ eff of the fluid composite is also found to be complex-valued and frequency-dependent. Note that energy dissipation is induced only by scattering between multiple scatterers. Therefore, the imaginary part of the effective parameters is only due to the diffusive scattering loss.

Effective interface conditions
The effective medium corresponding to a suspension of particles in a fluid can be described as a linearly viscoelastic fluid from the standpoint of coherent wave propagation (see, e.g., [28,46]). The viscoelastic fluid has mass density ρ eff and coherent acoustic waves propagate with the effective wavenumber K of equation (18). The reflection coefficient at the interface between an inviscid fluid and a viscoelastic fluid may then be written as Observe that ρ eff k/ρ 0 K is the ratio of the acoustic impedance of the viscoelastic fluid to that of the lossless surrounding fluid. Equation (45) is obtained by imposing continuity of pressure and normal velocity at the interface. Fikioris and Waterman [11] showed that these continuity conditions are justified, at least in the low-frequency limit, and effective parameters can be determined. (Note that pressure is given by the product of a scalar potential with density). A similar approach has been used in [28,46,47] in the framework of Waterman and Truell [26]. Effective interface conditions at the boundary of the half-space can also be derived in the analytical framework of Fikioris-Waterman (see [48] for two-dimensional (2D) problems).

Effective mass density and bulk modulus
Equating the reflection coefficients of equations (37) and (45), and using equation (A.7b), one finds that Equation (46) is the generalization of the Ament estimate for the effective mass density. It reduces exactly to equation (2) in the long-wavelength limit (see equation (B.13) in appendix B).
If one defines the effective bulk modulus κ eff by where the effective wavenumber K is given by equation (18) and the effective mass density ρ eff by equation (46), then κ eff is given by Once the effective scattering amplitudes F 0 and F π are determined from equations (36), the effective mass density and bulk modulus follow from equations (46) and (48).
Note that the low-frequency approximations of the effective mass density and bulk modulus are derived in appendix B, equations (B.5) and (B.6), respectively. Li and Chan [54] have derived the same formulae for ρ eff and κ eff by using an effective medium method (i.e. the coherent potential approximation). They showed that it becomes possible to achieve negative effective bulk modulus and negative effective mass density through resonance behaviour of the scattering coefficients T 0 (ω) and T 1 (ω). Structures are designed in which the resonances occur at very low frequency using silicone rubber spherical particles in water. Our general results for ρ eff and κ eff , equations (46) and (48), extend this previous analysis to include as many diffraction orders (monopole, dipole, quadrupole and so on) as needed for convergence, and are not restricted to the low-frequency range.
Note that equations (46) and (48) are the basis of two different closed-form formulae for effective parameters to be derived in the remainder of this section, depending on the level of approximation (i.e. low concentration and/or point particles).

Explicit formulae for ρ eff (ω) and κ eff (ω)
Here we combine the results of the previous two sections to give explicit formulae for the effective parameters ρ eff and κ eff valid in two limits: low concentration and point particles, respectively. First the low-concentration asymptotic approximation is considered. From equation (46), one obtains for the effective mass density where the coefficientsρ 1 andρ 2 are given in matrix notation bỹ In terms of the angular shape function f (θ) the coefficients (50) may be expressed as 13 By combining the expansions in ε of equations (21) and (49) with equation (47), the effective bulk modulus becomes κ eff κ 0 = 1 2 1 + ε where ξ 1 and ξ 2 are given by equations (22) (or (26)). We now consider the point-particle limit. In addition to a 3 1, we also require kb 1. We use the same procedure as that described in the previous sections. Equations (49) and (52) become, in terms of n 0 , The closed-form formulae (53) and (54) are direct consequences of the Lloyd-Berry estimate (31) for the effective wavenumber and the reflection coefficient (43) for the associated semi-infinite suspension of particles, which are both correct to quadratic terms in the concentration. Note that equation (49) (or (53)) agrees precisely with the small-φ estimate of equation (3) in the long-wavelength limit (see equation (B.21) in appendix B). This agreement provides a further check on the calculations. A comparison with previous work is also made in appendix C.
Apart from their dependence on k and n 0 , the effective dynamic parameters given by (31), (43), (53) and (54) are all completely determined when the angular shape function f (θ ) for an isolated particle is known. If this scattering amplitude can be determined analytically, numerically or experimentally, then the effective medium equivalent to the fluid-particle mixture is fully described.
It is worth noting that if one wants to study the behaviour of effective parameters at high concentrations (where such expansions are no longer appropriate), the general implicit equations (46) and (48) have to be used and/or more accurate pair-correlation functions have to be considered. In the following section, the entire analysis is generalized to account for arbitrary pair-correlation functions.

Generalization to an arbitrary pair-correlation function
The advantage of the quasi-crystalline approximation is that statistics, describing the spatial distribution of scatterers, can be incorporated into the effective field method. All the above developments start with the general equations (18) and (20). They are both given explicitly in terms of the multiple-scattering matrix R. Statement of the problem was completed by using the hole correction (11), for which the matrix R has elements given by equation (14). We shall now render our analysis more general; that is, an arbitrary pair-correlation function g(r ) (that satisfies the conditions (A.4)) is used to derive expressions for K, R, ρ eff and κ eff .
Note first that the form of R p in equation (16) is the result of hole correction. More generally, this function has values The integral in equation (55) describes the statistics of particle locations in terms of an arbitrary pair-correlation function g(r ). Next, in the limit of K → k, equation (55) gives Finally, in the point-particle limit, equation (56) reduces to The above equations enable us to calculate the effective dynamic parameters accounting for pair-correlated particles. If the simple hole correction (11) is used in equations (55)-(57), then equations (16), (24) and (29) are recovered.
To generalize equations (31), (43), (53) and (54) to an arbitrary pair-correlation function g(r ), one only needs to calculate e t QR 0 Qe and e t QJR 0 Qe in terms of the angular shape function f (θ). Combining equations (27) and (57) and using the generating function for Legendre polynomials, equation (30) is replaced with where One obtains where The effective wavenumber and the reflection coefficient are then obtained as Equations (63) and (64) are accurate to second order in concentration for point-like particles with an arbitrary two-point correlation function g(r ). Similar expressions can be derived for the effective mass density and bulk modulus by using equations (60) and (61).

Numerical illustrations and discussion
In order to illustrate the general behaviour of different approximations for effective dynamic parameters, some numerical examples are presented. Before starting, note that in a previous paper [33] we have rigorously compared the three solutions for the effective wavenumber (one implicit: equation (18), and two explicit: equations (21) and (31)). We have also highlighted differences in the effective wavenumber due to the choice of the pair-correlation function. It was noted that the hole correction yielded some unphysical results in some cases, whereas the Virial expansion [55] and the Percus-Yevick results [56] were, in the cases considered, physically acceptable. In the following, the pair-correlation function g(r ) satisfying the conditions (A.4) is chosen to be given by the Percus-Yevick relation (see [33]). Also, whenever the exclusion distance b (i.e. the distance of the closest approach between the centres of adjacent spheres) needs to be specified in the calculations, the value b = 2a is taken, which is thought to be physically reasonable.
We consider spherical particles of radius a, mass density ρ and longitudinal and shear sound speeds, c L and c T , respectively. Two different fluid-particle mixtures are considered: steel and polystyrene particles, both suspended in water 2 .

Reflection (transmission) from an effective layer of particles
Results are presented for the reflection coefficient R for a half-space of steel particles suspended in water. The reflection coefficient is calculated by using the general expression (37) and its two expansions in volume fraction φ : finite-size and point-like particles, equations (41) and (43), respectively. The comparison of these formulations is shown in figure 2, for two volume fractions φ = 0.15 and φ = 0.35. On the horizontal axis, the scattering dispersion parameter ω is defined asω = ka = 2πa/ , where is the wavelength of the plane incident wave in water. Thus, for a value ofω 1, the incident wavelength is approximately three times the diameter of the particle. It is therefore a fairly long wavelength. We observe from figure 2 that the values for the reflection coefficient predicted via the general expression (37) and the expansion (41) for finite-size particles are almost identical to each other, but differ slightly from the results based on the point-particle approximation (43), and this difference increases as the distribution becomes denser (φ = 0.35). The expansion (41) for finite-size particles is shown to be an excellent approximation of the true expression (37) even for quite dense distributions. As such, the expansion formulae of effective parameters for finite-size particles are exploited in the following. The explicit dependence of the reflection coefficient (45) on the effective mass density ρ eff means that, in principle, measurement of R via experiment can provide useful knowledge for its determination. Knowing the effective wavenumber K and mass density ρ eff , the effective acoustic impedance Z can easily be calculated, such as where z 0 = ρ 0 ω/k is the acoustic impedance of the inviscid fluid. With Z, one can describe acoustic wave propagation inside and outside a finite layer of particles suspended in a fluid by keeping track of the successive internal reflections from the boundaries of the equivalent homogeneous layer. There is an infinite number of reflections, which yields the reflection and transmission coefficients in the form of two infinite series. It is easy to interpret physically each term in the series. The two series converge, respectively, to the expressions [73] where h is the thickness of the layer and is the reflection coefficient of the associated semi-infinite layer.
It is important to note that the solution of the multiple-scattering problem derived in this paper does not require a priori knowledge of an effective mass density or of an effective bulk modulus for the mixture of particles inside the fluid. Furthermore, no assumptions concerning these quantities are needed. There exist works, however, where such assumptions are made in order to evaluate reflection and transmission on either side of a finite layer [57][58][59]. The resulting effective mass densities are real valued and frequency independent. Figure 3 illustrates the importance of using a complex-valued frequency-dependent mass density as in equation (49). This figure represents the moduli of the reflection and transmission coefficients corresponding toafluidlayercontainingsteelparticlesimmersedinwater(seefootnote2).Onthehorizontal axis, the effective wavelength eff of the plane coherent wave is defined as eff = 2π/Re K. The volume fraction of particles is φ = 0.35; two layer thicknesses are considered: h = a and h = 10a. Four results are compared, which correspond, respectively, to the four mass densities given by equations (1) transmission coefficients. This suggests that serious errors can result from using inappropriate expressions to model the effective mass density in experimental approaches. Only in the lowfrequency limit does the mass density formula reduce to that corresponding to the Ament estimate.

Acoustic metamaterials with random microstructure
Acoustic metamaterials, in which either the mass density or the bulk modulus is negative, can be used, for example, to design acoustic lenses to overcome the diffraction limit [60,61] or for the design of acoustic panels for sound insulation [62]. 'Double negative' acoustic metamaterials, in which both the mass density and the bulk modulus are negative, can cause negative refraction [63] and, as is known from electromagnetic wave theory, they can also be used to increase the resolution of conventional lenses [64,65].
The existence of frequency ranges where the effective medium presents negative constitutive parameters is related to resonances of the individual scatterers that constitute the metamaterial, i.e. either soft-particle resonances (e.g. [54]) or Helmholtz-like resonances (e.g. [66]). It is known that the monopolar resonances in the individual particles are responsible for negative bulk modulus and that the dipolar resonances are responsible for negative mass density (see, e.g., [54,69]). However, the full effect of the ensemble of scatterers that constitute the effective medium has only been explained partially, as higher diffraction orders (monopole, dipole, quadrupole and so on) are not included in the models; therefore such arguments are restricted only to the low-frequency range. In the low-frequency range, the effective medium is independent of whether the microstructure is random or periodic [21]. Not surprisingly, previous studies on acoustic metamaterials treated the microstructure as a regular lattice of scatterers, since for a periodic system, as opposed to a disordered one, the normal modes of the system may be easily obtained in the low-frequency limit (e.g. [67][68][69]).
In the following, we show that negative constitutive parameters can also be achieved by using suspensions of particles with random microstructures and are not restricted to the lowfrequency, long-wave limit. Here we report only one example of particles to realize acoustic metamaterials, i.e. homogeneous spherical particles. Obviously, more complex particles (e.g. homogeneous particles but with spherical anisotropy, oblate/prolate ellipsoidal particles, etc) could be used to enhance the frequency response, but a full analysis of this type is beyond the scope of this work. Figures 4 and 5 illustrate the effective mass density ρ eff /ρ 0 (=ρ eff + iρ eff ) and bulk modulus κ eff /κ 0 (=κ eff − iκ eff ) corresponding to a suspension of polystyrene particles in water (see footnote2).Ascanbeseen,inparticularfrequencyregions,boththeeffectivemassdensityand bulk modulus yield negative values in the vicinity of local resonances of individual particles. Figure 6 displays a contour plot of these negative regions. Also shown in figure 6 are the resonant frequencies for an isolated polystyrene sphere in water. Each resonance may be labelled by two symbols (n, ) in the Regge trajectories, where the first index defines the ordinal number of the resonance and the second indicates the family of surface waves. The = 1 family is usually referred to as the Rayleigh waves and the others (i.e. = 2, 3, . . .) are the Whispering Gallery waves. It is interesting to note that the negative regions for bothρ eff andκ eff are mainly due to the Rayleigh wave family of modes ( = 1). The first negative region forκ eff is related to the (n, ) = (2, 1) resonance, while forρ eff the first negative region is associated with the (n, ) = (3, 1) resonance. Note from figure 6 that the monopolar (n = 0) and dipolar (n = 1) resonances occur in a higher-frequency range. Figure 7 shows the behaviour of ρ eff /ρ 0 and κ eff /κ 0 versus the dimensionless frequencyω, for values of volume fraction φ ∈ (0.345, 0.355). The effect of every diffraction order is clearly seen in this figure. In particular, note that in the frequency region whereω ∈ (2.16, 2.22) (which corresponds to the (3, 1) resonance), bothρ eff andκ eff exhibit negative values simultaneously. In figure 8, the behaviour of ρ eff /ρ 0 and κ eff /κ 0 versus the volume fraction φ is presented for values of the dimensionless frequencyω ∈ (2.16, 2.22). In summary, figures 4 and 5 display a few narrow regions in which the effective mass density and/or bulk modulus are negative. These regions depend on both the frequency and the concentration of the suspended particles. To enlarge the negative regions of both ρ eff and κ eff , one could design the particles to broaden their resonant behaviour. For instance, considering heavier core particles, e.g. lead-cored polystyrene particles, will enhance the field oscillation in the polystyrene shell, which, in turn, will widen the resonant regions of ρ eff . On the other hand, if one engineers the particles to make them easier to deform and compress [70], e.g. air-filled polystyrene particles, this will widen the resonant regions of κ eff . The radius of the inner core of the particles is chosen such that the negative regions of both ρ eff and κ eff overlap. This and other applications will be exploited elsewhere.

Conclusion and perspectives
Using the multiple-scattering approach proposed by Fikioris and Waterman [11], we have determined the coherent wave motion in a non-viscous fluid containing disordered heterogeneities, such as particles, bubbles or contrast agents. Scatterers can be homogeneous spheres, layered, shell-like with encapsulated liquids or gas, non-absorbing, or absorbing. They are randomly distributed in a uniform half-space. Based on the quasi-crystalline approximation, the fluid-particle mixture is then shown to behave as an effective dissipative medium from the standpoint of the coherent wave motion. The effective medium is fully described once the effective wavenumber and the reflection coefficient for the associated semi-infinite suspension of particles are known. The effective mass density ρ eff and bulk modulus κ eff of the effective medium are determined by solving the general equations (46) and (48). As expected, ρ eff and κ eff are complex valued and frequency dependent. The imaginary part of the effective parameters is only due to the diffusive scattering loss. Equations (46) and (48) are the basis of two closed-form formulae for ρ eff and κ eff valid in two limits: low concentration and point particles. The lowconcentration expansions (49) and (52) are accurate to second order in the number of scatterers per unit volume, n 0 , for finite-size particles. The acoustic limiting case of small spheres is then considered, in the process taking kb 1. This limit corresponds to point particles with a spherical excluded volume (i.e. a hole of radius b). The resulting new formulae for ρ eff and κ eff , equations (53) and (54), depend only on the angular shape function f(θ) for an isolated particle, apart from their dependence on k and n 0 .
We have shown that the effective mass density and bulk modulus can be made negative near resonances by choosing appropriate resonant particles. The theory provides a technique for searching for acoustic metamaterials with specific desired properties. To obtain negative mass density ρ eff and bulk modulus κ eff , one must engineer the particles to coincide and enlarge the negative regions of both ρ eff and κ eff . One way to achieve this might be to simply change the dimensions of the particles or to consider core-shell particles. Another way is to consider a twophase suspension comprising a fluid matrix and a distribution of clustered particles (e.g. [72]). The application of this work to more complex acoustic metamaterials is the subject of ongoing research.
on the scatterer at r j . The average on disorder is taken over all configurations keeping the scatterer at r j fixed. Thus, it is a partial average. The integral of this scattered potential over all positions in the volume accessible to scatterers represents the total contribution of all scattered motions. The point at which the wave motion is evaluated is r.
When all the scatterers occupy deterministic positions, a scatterer at r i 'sees' an exciting wave motion that is the sum of the incident motion and the motions scattered by all the other scatterers. Under the 'quasi-crystalline approximation' [32], the integral equation for the exciting displacement potential is given by φ ex (r|r i ) = φ in (r) + n(r j |r i )T (r j ) φ ex (r|r j ) dv j , where n(r j |r i ) is the conditional number density of the scatterer at r j if a scatterer is fixed at r i . For a uniform and random array of spheres of radius a, the conditional number density must be specified in terms of a pair-correlation function, such that n(r j |r i ) = n 0 g(ρ ji ) (A.3) where the pair-correlation function g satisfies the conditions g(ρ ji ) = 0 for ρ ji < 2a, and lim The former condition holds for non-overlapping sets of spheres; the later condition is correct if the correlation of particle locations disappears when the distance between their centres tends to infinity.

A.2. The Fikioris-Waterman dispersion relation
Assuming that the exciting displacement potential (A.2) is a 'bounded' solution of the Helmholtz equation in the fluid and is a regular function at the point r i , one finds that the most general expression of φ ex is φ ex (r|r i ) = ∞ n=0 E n (r i )ψ n (ρ i ), (A.5) whereψ n (ρ i ) = j n (kρ i )P n (cos θ i ) is a regular spherical wavefunction and j n is a spherical Bessel function. The unknown modal coefficients E n (r i ) characterize the average displacement potential of the exciting motion on the scatterer at r i . In order to solve the ensemble average equation (A.2), the exciting coefficients E n (r) are decomposed into refracted plane waves in the forward direction (z > 0) with a propagation constant K. The effective wavenumber K defines the coherent motion according to E n (r) = X n exp(iKz). (A.6) Both K and X n are unknown. Substituting (A.5), with (A.6), into (A.2) and using the addition theorem for spherical wavefunctions [33] yields an infinite set of linear algebraic equations for the determination of K and the coefficients X n . The details of the calculation are not shown here since they merely reproduce results obtained in [11] and [33]. (B.14) The real parts of effective parameters in (B.12)-(B.14) are identified as the corresponding results due to Fikioris and Waterman [11]. The real part of equation (B.13) is identical to the Ament estimate (2), whereas the real part of equation (B.14) is recognized as the Reuss average. Note that the signs Im ρ eff > 0 and Im κ eff < 0 are in agreement with the sign of dissipation.

B.2. Second-order expansions
It is straightforward to expand the low-frequency formulae obtained in the previous section. Instead, the point-particle expressions (31), (43), (53) and (54) are approximated in the lowfrequency limit (ka 1). This provides an additional check of the correctness of the results obtained in this paper.
The leading order contribution to the angular shape function f(θ) clearly only involves T 0 and T 1 , i.e. Finally, using the results (B.16), the effective acoustic (K and R) and material (ρ eff and κ eff ) parameters are obtained as