A Proof that Multiple Waves Propagate in Ensemble-Averaged Particulate Materials
rrspa.royalsocietypublishing.org
Research
Article submitted to journal
Subject Areas:
Wave motion, acoustics,mathematical physics
Keywords: wave propagation, random media,composite materials, backscattering,multiple scattering, ensembleaveraging, Wiener-Hopf
Author for correspondence:
Artur L. Gowere-mail: arturgower @ gmail.comwebsite: arturgower.github.io A Proof that Multiple WavesPropagate inEnsemble-AveragedParticulate Materials
Artur L. Gower , I. David Abrahams , andWilliam J. Parnell Department of Mechanical Engineering, TheUniversity of Sheffield, UK, Isaac Newton Institute for Mathematical Sciences, 20Clarkson Road, Cambridge CB3 0EH, UK School of Mathematics, University of Manchester,Oxford Road, Manchester M13 9PL, UK
Effective medium theory aims to describe a complexinhomogeneous material in terms of a few importantmacroscopic parameters. To characterise wave propagationthrough an inhomogeneous material, the most crucialparameter is the effective wavenumber . For this reason,there are many published studies on how to calculatea single effective wavenumber. Here we present aproof that there does not exist a unique effectivewavenumber; instead, there are an infinite numberof such (complex) wavenumbers. We show thatin most parameter regimes only a small numberof these effective wavenumbers make a significantcontribution to the wave field. However, to accuratelycalculate the reflection and transmission coefficients,a large number of the (highly attenuating) effectivewaves is required. For clarity, we present results forscalar (acoustic) waves for a two-dimensional materialfilled (over a half space) with randomly distributedcircular cylindrical inclusions. We calculate theeffective medium by ensemble averaging over allpossible inhomogeneities. The proof is based onthe application of the Wiener-Hopf technique andmakes no assumption on the wavelength, particleboundary conditions/size, or volume fraction. Thistechnique provides a simple formula for the reflectioncoefficient, which can be explicitly evaluated formonopole scatterers. We compare results with analternative numerical matching method. c (cid:13) The Authors. Published by the Royal Society under the terms of theCreative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author andsource are credited. a r X i v : . [ phy s i c s . c l a ss - ph ] A ug r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
1. Introduction
Materials comprising particles or inclusions that are randomly distributed inside a uniform hostmedium occur frequently in the world around us. They occur as synthetically fabricated mediaand also in nature. Common examples include composites, emulsions, suspensions, complexgases, and polymers. Understanding how electromagnetic, elastic, or acoustic waves propagatethrough these materials is necessary in order to characterise the properties of these materials, andalso to design new materials that can control wave propagation.The wave scattered from a particulate material will be influenced by the positions andproperties of all particles, which are usually unknown. However, this scattered field, averagedover space or over time, depends only on the average particle properties. Many measurementsystems perform averaging over space, if the receivers or incident wavelength are largeenough [1], or over time [2]. In most cases, this averaging process is the same as averaging overall possible particle configurations. Such systems are sometimes called ergodic [2,3]. In this paper,we focus on ensemble averaged waves, satisfying the scalar wave equation in two-dimensions,reflecting from, and propagating in, a half-space particulate material. In certain scenarios, such aslight scattering, it is easier to measure the average intensity of the wave. However, even in thesecases, the ensemble-averaged field is often needed as a first step [4,5].One driving principle, often used in the literature, is that the ensemble-averaged wave itselfsatisfies a wave equation with a single effective wavenumber [6–8]. Reducing a inhomogeneousmaterial, with many unknowns, down to one effective wavenumber is attractive as it greatlyreduces the complexity of the problem. For this reason many papers have attempted todeduce this unique effective wavenumber from first principles in electromagnetism [3,9,10],acoustics [11–15] and elasticity [16,17]. See [18] for a short overview of the history of this topic,including typical statistical assumptions employed within the methods, such as hole-correctionand the quasi-crystalline approximation, which we also adopt here.The assumption that the ensemble averaged wave field satisfies a wave equation, with aneffective wavenumber, has never been fully justified. Here we prove that there does not exista unique effective wavenumber but instead there are an infinite number of them. Gower etal. [18] first showed that there exist many effective wavenumbers, and provided a technique,the
Matching Method , to efficiently calculate the effective wave field. In the present paper and [18],we show that for some parameter regimes, at least two effective wavenumbers are needed toobtain accurate results, when compared with numerical simulations. We also provide examplesof how a single effective wave approximation leads to inaccurate results for both transmissionand reflection for a halfspace filled with particles, see Figure 1.Although the
Matching Method developed in [18] gave accurate results, when comparedto numerical methods and known asymptotic limits, the limitations of the method were notimmediately clear. Here however we illustrate that the Matching Method is robust, becausecombining many effective wavenumbers is not just a good approximation, it is an analyticalsolution to the integral equation governing the ensemble averaged wave field. We prove thisby employing the Wiener-Hopf technique and then, for clarity, illustrate the solution for particlesthat scatterer only in their monopole mode. The Wiener-Hopf technique also gives a simple andelegant expression for the reflection coefficient.The Wiener-Hopf technique is a powerful tool to solve a diverse range of wave scatteringproblems, see [19, Chapter 5. Wiener-Hopf Technique] and [20,21] for an introduction. It isespecially useful for semi-infinite domains [22–27] and boundary value problems of mixedtype. In this work, the Wiener-Hopf technique clearly reveals the form of the analytic solution,but to compute the solution would require an analytic factorisation of a matrix-function. Toexplicitly perform this factorisation is difficult [28–31]. Indeed this is often the hardest aspectof employing the Wiener-Hopf technique, although there exist approximate methods for thispurpose [28,32–34]. We do not focus in this article on these analytic factorisations, as there already r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... exists a method to compute the required solution [18]. Instead, the present work acts as proof thatthe Matching Method [18] faithfully reproduces the form of the analytic solution. Ensembled averagedparticulate material Re s Re s Re s Re s k θ in θ in y x Figure 1.
When an incident plane wave e i k · ( x,y ) , with k = k (cos θ inc , sin θ inc ) , encounters an (ensemble-averaged)particulate material, it excites many transmitted plane waves and one reflected plane wave. The transmitted waves are ofthe form e i s p · ( x,y ) with wavenumbers s p = S p (cos θ p , sin θ p ) where both S p and θ p are complex numbers. The largerIm s p , the more quickly the wave attenuates as it propagates into the half-space and the smaller the drawn vector forthat wave above. The results shown here represent the effective wavenumbers for parameters (5.2) , which are shownin Figure 3. Figure 1 shows the main setup and result of this paper: an incident plane wave excites the half-space x > filled with ensemble-averaged particles (the blue region), which generates a reflectedwave and many effective transmitted waves. The s p are the transmitted wavevectors, and thesmaller the length of the vector, the faster that effective wave attenuates as it propagates furtherinto the material.The paper begins by summarising the equations that govern ensemble averaged waves intwo-dimensions in section 2. Following this, in section 3 we apply the Wiener-Hopf technique tothe governing integral equation and deduce that the solution is a superposition of plane waves,each with a different effective wavenumber. A simple expression for the reflection coefficient isalso derived. In section 4 we specialise the results for particles that scatter only in the monopolemode, which leads to a closed form analytic solution.The dispersion relation (3.30), derived in section 3, admits an infinite number of solutions, theeffective wavenumbers. In section 5, we deduce asymptotic forms for the effective wavenumbersin both a low and high frequency limit. In section 6 we compare numerical results for monopolescatterers, using the Wiener-Hopf technique, with classical methods that assume only oneeffective wavenumber [11,13], and the Matching Method introduced in [18]. In general, whencomparing predicted reflection coefficients, the Wiener-Hopf and Matching Method agree well,whereas the classical single-effective-wavenumber method can disagree by anywhere up to .These results are discussed in section 7 together with anticipated future steps. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
2. Waves in ensemble averaged particles
Consider a region filled with particles or inclusions that are uniformly distributed. The field u isgoverned by the scalar wave equations: ∇ u + k u = 0 , (in the background material) , (2.1) ∇ u + k o u = 0 , (inside a particle) , (2.2)where k and k o are the real wavenumbers of the background and inclusion materials, respectively.We assume all particles are identical, except for their position and orientation, for simplicity. Fora distribution of particles, or multi-species, see [15].Our goal is to calculate the ensemble average field (cid:104) u ( x, y ) (cid:105) , that is, the field averaged overall possible particle positions and orientations. For clarity, and ease of exposition, we considerthat the particles are equally likely to be located anywhere except that they cannot overlap (this isoften called the hole correction assumption). We also assume the quasi-crystalline approximation;for details on this, and for further details on deducing the results in this section, see [11,15,18].By splitting the total steady wave field u ( x, y ) into a sum of the incident wave u inc ( x, y ) andwaves scattered by each particle, the j th scattered wave being u j ( x, y ) , we can write: u ( x, y ) = u inc ( x, y ) + (cid:88) j u j ( x, y ) . (2.3)A simple and useful scenario to consider is when all particles are placed only within the half-space x > , which are then excited by a plane wave, with implicit time dependence e i ωt , incidentfrom a homogeneous region: u inc ( x, y ) = e i ( αx + βy ) , with ( α, β ) = ( k cos θ inc , k sin θ inc ) , (2.4)where we restrict the incident angle − π < θ inc < π , as shown in Figure 1, and consider a slightlydissipative medium with Re k > and Im k > . (2.5)This dissipation will facilitate the use of the Wiener-Hopf technique, and after reaching thesolution we can take k to be real .To describe the particulate medium we employ the following notation: b = the minimum distance between particle centres , (2.6) n = number of particles per unit area , (2.7) T n = the coefficients of the particle’s T-matrix , (2.8) φ = π n b particle area fraction . (2.9)Although the area fraction φ , normally called the volume fraction, is a combination of otherparameters, it is useful because it is non-dimensional. If we let a o be the maximum distancefrom the particle’s centre to its boundary, then we can set b = γa o , where γ ≥ so as to avoidtwo particles overlapping. The volume fraction that does not include the exclusion zone φ (cid:48) , asused in [18, equation (4.7)], is then φ (cid:48) = 4 φ/γ .The T n are the coefficients of a diagonal T-matrix [35–39]. The T-matrix determines howthe particle scatters waves, and so depends on the particle’s shape and boundary conditions.A diagonal T-matrix can be used to represent either a radially symmetric particle, or particlesaveraged over their orientation, assuming the orientations have a random uniform distribution.The results of ensemble averaging (2.3) from first principals are deduced in a number ofreferences [15,18] and so deails of this procedure are omitted here for brevity. To represent the The case where particles can be placed anywhere in the plane can lead to ill-defined integrals [11]. Assuming Im k = (cid:15) > , rather than (cid:15) ≥ , will facilitate calculating certain integrals that appear below. However, afterreaching a solution, we can take the limit (cid:15) → to recover the physically viable solution for Im k = 0 . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... ensemble averaged scattered wave from a particle, whose centre is fixed at ( x , y ) , we use (cid:104) u ( x + X, y + Y ) (cid:105) ( x ,y ) = ∞ (cid:88) n = −∞ A n ( x ) e i βy H (1) n ( kR ) e i nΘ , (2.10)for R := √ X + Y > b/ , so that ( X, Y ) is on the outside of this particle, with ( R, Θ ) being thepolar coordinates of ( X, Y ) , H (1) n are Hankel functions of the first kind, and A n is some field wewant to determine .By choosing x < − b , which is outside of the region filled with particles, then taking theensemble average on both sides of (2.3) results in equation (6.7) of [18], given by: (cid:104) u ( x, y ) (cid:105) = u inc ( x, y ) + R e − i αx + i βy for x < − b, (2.11)which is the incident wave plus an effective reflected wave with reflection coefficient: R = e i αx n ∞ (cid:88) n = −∞ ˆ ∞ A n ( x ) ψ n ( x − x ) dx , (2.12)where we assumed particles are distributed according to a uniform distribution, and the kernel ψ n is given by ψ n ( X ) = ˆ Y >b − X e i βY ( − n H (1) n ( kR ) e i nΘ dy. (2.13)Later we show that, as expected, R is independent of x .The system governing A m ( x ) is given by equation (4.7) of [18]: n T m ∞ (cid:88) n = −∞ ˆ ∞ A n ( x ) ψ n − m ( x − x )d x = A m ( x ) − e i αx T m e i m ( π/ − θ inc ) , for x ≥ , (2.14)for all integers m . Kristensson [40, equation (15)] presents an equivalent integral equation forelectromagnetism and particles in a slab.Our main aim is to reach an exact solution for A n ( x ) by employing the Wiener-Hopf techniqueto (2.14). We show how this also leads to simple solutions for the reflection coefficient byusing (2.11). We acknowledge the authors of [11], as they noticed that (2.14) is a Wiener-Hopfintegral equation; but apparently did not follow the steps indicated in the following sections .
3. Applying the Wiener-Hopf technique
Equation (2.14) is convolution integral equation with a difference kernel. This means applying aFourier transforms can lead to elegant and simple solutions. To facilitate, we must analyticallyextend (2.14) for all x ∈ R by defining n T m ∞ (cid:88) n = −∞ ˆ ∞ A n ( x ) ψ n − m ( x − x )d x = (cid:40) A m ( x ) − e i αx T m e i m ( π/ − θ inc ) , x ≥ ,D m ( x ) , x < , (3.1)for integers m , where if the A n ( x ) were known for x > , then the D n ( x ) would be given fromthe left hand-side. Note that the kernel ψ n defined in (2.13) is already analytic in the domain R . The factor e i βy appears due to the translational invariance of (cid:104) u ( x , y ) (cid:105) in y , which is a result of the material beingstatistically homogeneous, see [15] for details. We define the reflection coefficient only for x < − b , instead of x < − b/ , so that we can use ψ n in the formula for R , whichwill in turn facilitate calculating R . However, they were unable to solve it because, it seems, of a mistake in the integrand of equation (37) of [11]; they used e i βY where they should have used cos( βY ) . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... The field D ( x ) is not just an abstract construct, it is closely related to the reflected wave: bydirectly comparing (3.1) with the reflection coefficient (2.12), for x < − b , we find that D ( x ) = T R e − i αx . (3.2)To solve (3.1) we employ the Fourier transform and its inverse, which we define as ˆ f ( s ) = ˆ ∞−∞ f ( x ) e i sx dx with f ( x ) = 12 π ˆ ∞−∞ ˆ f ( s ) e − i sx dx, (3.3)for any smooth function f . We then define ˆ A + n ( s ) = ˆ ∞ A n ( x ) e i sx dx, ˆ D − n ( s ) = ˆ −∞ D n ( x ) e i sx dx. (3.4)We can determine where ˆ A + n and ˆ D − n are analytic by assuming that | A n ( x ) | < e − xc for x → ∞ , (3.5) | D n ( x ) | < e xc for x → −∞ , (3.6)for some (possibly small) positive constant c . This leads to ˆ A + n ( s ) being analytic for Im s > − c ,while ˆ D − n ( s ) is analytic for Im s < c . In other words, both ˆ A + n ( s ) and ˆ D − n ( s ) are analytic in theoverlapping strip | Im s | < c. (3.7)To apply the Wiener-Hopf technique we also need to specify the large s behaviour for both ˆ A + n ( s ) and ˆ D + n ( s ) . To achieve this, we assume, on physical grounds, that A n ( x ) is bounded when x → + , and D n ( x ) is bounded when x → − . Then, it can be shown [21,41] that ˆ A + n ( s ) = O ( | s | − ) and ˆ D − n ( s ) = O ( | s | − ) for | s | → ∞ , (3.8)in their respective half-planes of analyticity.Applying a Fourier transform to both sides of equation (3.1), the left-hand side becomes ˆ ∞ A n ( x ) ˆ ∞−∞ ψ n − m ( x − x ) e i sx d x d x = ˆ A + n ( s ) ˆ ψ n − m ( s ) , (3.9)in which ˆ ψ n ( s ) is well defined (i.e. analytic) for s in the strip: | Im s | < (1 − | sin θ inc | )Im k, (3.10)see appendix A for details. The right-hand side of (3.1) becomes ˆ −∞ D m ( x ) e i sx dx + ˆ ∞ A m ( x ) e i sx d x − e i m ( π/ − θ inc ) T m ˆ ∞ e i x ( s + α ) d x = ˆ D − m ( s ) + ˆ A + m ( s ) − T m ie i m ( π/ − θ inc ) ( s + α ) + , (3.11)where for the last step we assumed Im ( s + α ) > , which is why we use the superscript + on ( s + α ) + . This assumption, together with (3.7) and (3.10), is satisfied if | Im s | < (cid:15), where (cid:15) = min { c, (1 − | sin θ inc | )Im k, Im α } . (3.12)If (3.12) is satisfied then we can combine (3.9),(3.11) and (A 6), to obtain the Fourier transformof (3.1) in matrix form: Ψ ( s ) ˆ A + ( s ) s − α = − ˆ D − ( s ) + B s + α , (3.13) The solutions for A n ( x ) and D n ( x ) , in the next section, show that these assumptions do hold. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... where ˆ A + ( s ) and ˆ D − ( s ) are vectors with components ˆ A + n ( s ) and ˆ D − n ( s ) , respectively and B m = i T m e i m ( π/ − θ inc ) , (3.14) Ψ mn ( s ) = G mn ( S )( − i ) n − m e i ( n − m ) θ S , (3.15) G mn ( S ) = ( s − α ) δ mn + 2 π n T m N n − m ( bS ) , (3.16) N m ( bS ) = bk J m ( bS )H (1) (cid:48) m ( bk ) − bS J (cid:48) m ( bS )H (1) m ( bk ) . (3.17)where, for reference, Ψ mn ( s ) = ( s − α ) (cid:104) δ mn − n T m ˆ ψ n − m ( s ) (cid:105) , (3.18)and ˆ ψ n − m ( s ) is given by (A 6). In the above θ S and S are chosen to satisfy s = S cos θ S with S sin θ S = k sin θ inc . (3.19)Later we identify S and θ S as the effective wavenumber and transmission angle. The above doesnot determine the sign of S for any given complex s . To fully determine S and θ S , we take sgn(Re s ) = sgn(Re S ) which together with (3.19) leads to θ S = arctan (cid:18) k sin θ inc s (cid:19) , S = (cid:113) s + ( k sin θ inc ) , (3.20)where both S and θ S , when considered as functions of s , contain branch-points at s = ± i k sin θ inc with finite branch-cut running between − i k sin θ inc and i k sin θ inc . However, Ψ ( s ) is an entirematrix function having only zeros in s and no branch-points; see the end of appendix A for details.Determining the roots of det Ψ ( s ) = 0 will be a key step in solving (3.13), and so the followingidentities will be useful Ψ mn ( − s ) T n = Ψ mn ( s ) T n ( − m − n e i ( m − n ) θ s = Ψ nm ( s ) T m , (3.21) det Ψ ( − s ) = det Ψ ( s ) and det Ψ ( s ) = det G ( S ) . (3.22)where (3.21) results from (3.18) and (A 10). Equation (3.22) then follows from using (3.21) , (3.15),and appendix C. (a) Multiple waves solution To solve (3.13), we use a matrix product factorisation [42] of the form: Ψ ( s ) = Ψ − ( s ) Ψ + ( s ) , (3.23)where Ψ − ( s ) , and its inverse, are analytic in Im s < (cid:15) , and Ψ + ( s ) , and its inverse, are analytic forIm s > − (cid:15) . See (3.12) for the definition of (cid:15) .For our purposes, it is enough to know that such a factorisation exists [42], as this will lead toa proof that A ( x ) is a sum of attenuating plane waves.Multiplying both sides of (3.13) by [ Ψ − ( s )] − and by ( s − α ) − leads to Ψ + ( s ) ˆ A + ( s )( s + α ) + = − ( s − α ) − [ Ψ − ( s )] − ˆ D − ( s ) + [ Ψ − ( s )] − B ( s − α ) − ( s + α ) + , (3.24)where ( s + α ) + is analytic for Im s > − Im α , while ( s − α ) − is analytic for Im s < Im α . We needto rewrite the last term above as a sum of a function which is analytic in the upper half-plane (Im With our choice of branch-cut, the simplest way to compute S , with most software packages, is to use S =sgn(Re s ) (cid:112) s + ( k sin θ inc ) with the default cut location for the square root function, i.e. along the real negative line. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... s > − (cid:15) ) and another analytic in the lower half-plane. This is achieved below: [ Ψ − ( s )] − B ( s − α ) − ( s + α ) + = − α ( s + α ) + [ Ψ − ( − α )] − B (cid:124) (cid:123)(cid:122) (cid:125) g + ( s ) + [ Ψ − ( s )] − B ( s − α ) − ( s + α ) + + [ Ψ − ( − α )] − B α ( s + α ) + (cid:124) (cid:123)(cid:122) (cid:125) g − ( s ) , (3.25)where we define lim s →− α g − ( s ) = (cid:20) I + 2 α [ Ψ − ( − α )] − d Ψ − d s ( − α ) (cid:21) [ Ψ − ( − α )] − B , so that g − ( s ) does not have a pole at s = − α and is therefore analytic for Im s < (cid:15) .Substituting (3.25) into (3.24) leads to Ψ + ( s ) ˆ A + ( s )( s + α ) + + g + ( s ) = − ( s − α ) − [ Ψ − ( s )] − ˆ D − ( s ) + g − ( s ) . (3.26)Because both sides are analytic in the strip | Im s | < (cid:15) , we can equate each side to E ( s ) , someanalytic function in the strip. Further, as the left-hand side (right-hand side) of (3.26) is analyticfor Im s > (cid:15) ( Im s < − (cid:15) ), we can analytically continue E ( s ) for all s , i.e. E ( s ) is entire.To determine E ( s ) we need to estimate its behaviour as | s | → ∞ . From (3.8) we have that A + ( s ) = ( | s | − ) as | s | → ∞ in the upper half-plane, and from (3.15 - 3.17): Ψ ( s ) = ( s − α ) I + O ( | s | ) as | s | → ∞ , (3.27)for s in the strip (3.12). From this we know that the factors Ψ + ( s ) and Ψ − ( s ) must be O ( | s | ) as | s | → ∞ , in their respective half-planes of analyticity [28]. So, the left hand-side of (3.26) behavesas O ( | s | − ) as | s | → ∞ in Im s > − (cid:15) . We can therefore use Liouville’s theorem to conclude that E ( s ) ≡ , which means the Wiener-Hopf equation (3.26) is formally equivalent to ˆ A + ( s ) = − α [ Ψ + ( s )] − [ Ψ − ( − α )] − B , (3.28) ˆ D − ( s ) = Ψ − ( s ) g − ( s )( s − α ) − . (3.29)Let C + ( s ) be the cofactor matrix of Ψ + ( s ) , so that [ Ψ + ( s )] − = [ C + ( s )] T det( Ψ + ( s )) . From the property (3.22) we can write det Ψ ( s ) = f ( s ) for some function f . Then, for every root s = s p of det Ψ ( s ) , with Im s p > , we have that − s p is also a root, and vice-versa. From hereonwards we assume: det Ψ ( s p ) = det Ψ ( − s p ) = 0 with Im s p > and p = 1 , , · · · , ∞ . (3.30)For any truncated matrix Ψ ( s ) , i.e. evaluating m, n = − M, . . . , M in (3.15), the roots s p arediscrete. In section 5 we demonstrate asymptotically that they are indeed discrete for the limits oflow and high wavenumber k . For the numerical results presented in this paper, we numericallysolve the above dispersion relation for the truncating the matrix Ψ ( s ) , and then increase M untilthe roots converge (typically no more than M = 4 was required).Given det Ψ ( s ) = det Ψ − ( s ) det Ψ + ( s ) , every root of det Ψ ( s ) must either be a root of det Ψ − ( s ) or a root of det Ψ + ( s ) . For [ Ψ + ( s )] − to be analytic in the upper half-plane, det Ψ + ( s ) must only have roots s = − s p . As a consequence, det Ψ − ( s ) only has roots s = s p . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... To use the residue theorem below, we need to calculate det Ψ + ( s ) for s close to the root − s p ,in the form det Ψ + ( s ) = det Ψ + ( − s p ) + ( s + s p ) d det Ψ + d s ( − s p ) + O (( s + s p ) )= s + s p det Ψ − ( − s p ) d det Ψ d s ( − s p ) + O (( s + s p ) ) (3.31)where we use d det Ψ ds ( − s p ) instead of d det Ψ + ds ( − s p ) det Ψ − ( − s p ) , because it is more difficult tonumerically evaluate d det Ψ + ds ( − s p ) . C D C A Im s Re s α − s p Figure 2.
An illustration of the contour integral over C D , used to calculate (3.34) for x < , and the contour integral over C A , used to calculate (3.32) for x > . The − s p (the red points) are roots of (3.30) , and also the poles of (3.28) . Thesingle blue point α is the only pole of (3.29) . Using the above, and that C + ( S ) is analytic for Im s > − (cid:15) , we can apply an inverse Fouriertransform (3.3) to both sides of (3.28) and using residue calculus we find A ( x ) = − απ ˆ ∞−∞ [ C + ( s )] T [ Ψ − ( − α )] − B det Ψ + ( s ) e − i sx d s = (cid:40)(cid:80) ∞ p =1 A p e i s p x , x > , , x < , (3.32)with A p = 2 α i det Ψ − ( − s p ) d det Ψ d s ( − s p ) [ C + ( − s p )] T [ Ψ − ( − α )] − B . (3.33)For x > , the integral over s ∈ [ −∞ , ∞ ] in (3.32) is, by Jordan’s lemma, the same as a clockwiseintegral over the closed contour C A which surrounds the poles − s , − s , . . . , i.e. roots of (3.30),as shown by Figure 2. Note that the cofactor matrix C + ( s ) contains no poles and so does notcontribute additional residual terms. The yellow striped region in Figure 2 is the domain where Ψ is analytic. On the other hand, for x < , the integral (3.32) is the same as an integral overthe counter-clockwise closed contour within the region Im s > (not shown in Figure 2). Theintegrand has no poles in this domain and hence evaluates to zero.Likewise, by applying an inverse Fourier transform to (3.29), we obtain: D ( x ) = 12 π ˆ ∞−∞ Ψ − ( s ) g − ( s )( s − α ) − e − i sx d s = (cid:40) i Ψ − ( α )[ Ψ − ( − α )] − B e − i αx , x < , , x > , (3.34)For x < the above integral is the same as a counter-clockwise closed integral over C D whichsurrounds the pole s = α (recalling that Im α > ), as shown in Figure 2. The result is just the r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... residue at this pole. That is, the function Ψ − ( s ) g − ( s ) contains no other singularities within Im s > . On the other hand, for x > the integral is the same as a closed clockwise integral around theregion Im s < which evaluates to zero, as there are no singularities in this region (not shown inFigure 2).Clearly (3.32) shows that A ( x ) is a sum of plane waves with different effective wavenumbers s p , each satisfying (3.30). In section 5 we discuss these roots in more detail, and in section 6, wesee that usually only a few effective wavenumbers are required to obtain accurate results. (b) Reflection coefficient By substituting (3.34) in (3.2) leads to R = i T − ∞ (cid:88) n,m = −∞ Ψ − n ( α )[ Ψ − ( − α )] − nm B m . (3.35)Alternatively, the reflection coefficient can be calculated from (2.12) by employing the form of A ( x ) from (3.32), which is the more common approach. To simplify, we use ψ n ( X ) = ( − n ˆ ∞−∞ e i kY sin θ inc H (1) n ( kR ) e i nΘ dY = 2 α i n e − i nθ inc e i αX for X > , (3.36)which then implies that ψ n ( x − x ) = α i n e − i nθ inc e i α ( x − x ) for x ≥ x . The above is shown in [43,equation (37)] and [11, equation (65)]. This result together with (3.32) substituted into (2.12) leadsto the form R = 2 n α ∞ (cid:88) n = −∞ i n e − i nθ inc ˆ ∞ A n ( x ) e i αx dx = 2 i n α ∞ (cid:88) n = −∞ ∞ (cid:88) p =1 i n e − i nθ inc A pn s p + α , (3.37)where we used that Im s p > . The above agrees with [13, equation (39)] and [18, equation (6.9)].
4. Monopole scatterers
For particles that scatter only in their monopole mode, i.e. the scattered waves are angularlysymmetric about each particle, we can easily calculate the factorisation (3.23). This type ofscattered wave tends to dominate in the long wavelength limit for scatterers with Dirichletboundary conditions. In acoustics, these correspond to particles with low density or low soundspeed.Once we know the factorisation (3.23), we can then calculate the average scatteringcoefficient (3.32) and average reflection coefficient (3.35). We will compare both of these againstpredictions from other methods in section 6. (a) Wiener-Hopf factorisation
For scalar problems, there are well known techniques to factorise Ψ ( s ) = Ψ − ( s ) Ψ +00 ( s ) , such asCauchy’s integral formulation, for details see [19, Section 5. Wiener-Hopf Technique] and [21].For monopole scatterers we use S − k = s − α and rewrite Ψ ( s ) = ( s − α ) q ( s ) , with q ( s ) = 1 + 2 π n T N ( bS ) S − k , When taking a zero thickness boundary layer, i.e. J = 0 , and appropriate substitutions. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... with N ( bS ) given by (3.17). Then, because q ( s ) → as | s | → ∞ , we can factorise q ( s ) = q − ( s ) q + ( s ) using q + ( s ) = exp (cid:18) π i & ∞−∞ log q ( z ) z − s dz (cid:19) , (4.1) q − ( s ) = exp (cid:18) − π i $ ∞−∞ log q ( z ) z − s dz (cid:19) , (4.2)where the integral path for q + ( s ) ( q − ( s ) ) has to be in the strip where q ( s ) is analytic, with thepath for q + ( s ) ( q − ( s ) ) passing below (above) z . We then have that Ψ +00 ( s ) = ( s + α ) + q + ( s ) , Ψ − ( s ) = ( s − α ) − q − ( s ) , Ψ − ( − s ) = − Ψ +00 ( s ) , (4.3)where (4.3) holds if − s is below the integration path of (4.2) and s is above the integrationpath of (4.1). From (3.32) we see that we need only evaluate Ψ +00 ( s ) , and therefore q + ( s ) , for s = s , s , . . . , s p where as p increases, the s p become more distant from the real line. Then forlarge z , by inspection of (3.17), we have that (cid:12)(cid:12)(cid:12)(cid:12) log q ( z ) z − s (cid:12)(cid:12)(cid:12)(cid:12) ∼ z / | z − s | , and therefore we can accurately approximate the integral (4.1) by truncating the integrationdomain for large z . (b) Explicit solution for monopole scatterers. For monopole scatterers A n ( x ) = D n ( x ) = 0 for | n | > . Using this in (3.14 - 3.17) leads to allvectors and matrices having only one component, given by setting n = m = 0 . In this case A (3.32)reduces to A ( x ) = ∞ (cid:88) p =1 A p e i s p x with A p = 2 αT Ψ +00 ( α ) Ψ +00 ( s p ) d Ψ d s ( s p ) = T s p − α q + ( s p ) q + ( α ) q (cid:48) ( s p ) , (4.4)for x > , where we used (4.3), C + ( s ) = 1 , B = i T , and dΨ ds ( − s ) = − dΨ ds ( s ) for every s .Likewise for (3.35) we arrive at R = Ψ ( α )( Ψ +00 ( α )) = π n T N ( bα )2( αq + ( s )) . (4.5)Alternatively, using (3.37), we can calculate the contribution of P effective waves to thereflection coefficient R P = 2 i n α P (cid:88) p =1 A p s p + α = 2 i n T αq + ( α ) P (cid:88) p =1 s p − α q + ( s p ) q (cid:48) ( s p ) with R = lim P →∞ R P , (4.6)where the error | R P − R | then indicates how many effective waves are needed to accuratelydescribe the field near the boundary x = 0 .
5. Multiple effective wavenumbers
Equation (3.32) clearly shows that A ( x ) is a sum of attenuating plane waves, each with a differenteffective wavenumber s p . These s p satisfy the dispersion equation (3.30): det Ψ ( s p ) = det G ( S p ) = 0 , (5.1)with Ψ given by (3.16) and the first identity follows from (3.22). Note that the factors q + ( s ) and q − ( s ) are singularity and pole free in their respective regions of analyticity, and so theirinverses [ q + ( s )] − and [ q − ( s )] − have the same property. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... − −
30 0 30 600123456 backward forward Re S I m S − − . . . . S S Re S Figure 3.
Examples fo effective wavenumbers S p which satisfy the dispersion equation (5.1) with the properties (5.2) .The blue points represent waves travelling forwards (i.e. deeper into the material), while the red represent waves travellingbackwards. All these waves are excited in a reflection experiment. Two wavenumbers in particular stand out as having thelowest attenuation S and S , both inside the grey dashed circle. The graph on the right is a magnification of the regionclose to these two wavenumbers. Out of these two, most efforts in the literature have focused on calculating S , as itoften has the lowest imaginary part; however for this case, because S has a smaller attenuation it will have a significantcontribution to both transmission and reflection. An important conclusion from det G ( S p ) = 0 is that the wavenumbers S p are independent ofthe angle of incidence θ inc . We focus on showing the results for S p , rather than s p , because thenwe do not need to specify θ inc .As a specific example, let us consider circular particles with Dirichlet boundary conditions (i.e.particles with zero density or soundspeed), and the parameters T n = − J n ( ka o )H (1) n ( ka o ) , kb = 1 . , ka o = 0 . , φ = 30% , (5.2)where a o is the radius of the particle.With the above parameters, we found that truncating the matrix Ψ ( s ) , with | n | ≤ and | m | ≤ in (3.15-3.17), led to accurate results when calculating the effective wavenumbers S p , i.e. the rootsof (3.30). Numerically calculating the wavenumbers S p then leads to Figure 3.The effective wavenumbers with the lowest attenuation (smallest imaginary part) contributethe most to the transmitted wave. In Figure 3 we see two wavenumbers have lower attenuationthen the rest, both within the dashed grey circle. The blue point represents the wavenumber thatmost of the literature focuses on calculating: it has a positive real part and therefore propagatesforwards along the x − axis (into the material) as is expected for a transmitted wave. However, theother wavenumber, with negative real part, is equally as important because it actually has lowerattenuation. Figure 1 illustrates several effective wavenumbers, some travelling forward into thematerial, while others have negative phase direction (travel backwards).In Figure 3 we see what appears to be an infinite sequence of effective wavenumbers S p , where | S p | → ∞ as p → ∞ . To confirm their existence, and to find their locations as | p | → ∞ , we developasymptotic formulas in appendix B. The results of the asymptotics are summarised below.For monopole scatterers , where n = m = 0 in (3.15), equations (A 7) give the effectivewavenumbers S op at leading order: bS o ± p = σ ± p + i log (cid:32) | σ ± p | / r c (cid:33) , σ + p = θ c + 2 πp for p > − (cid:108) θ c π (cid:109) ,σ − p = θ c − π − πp for p > (cid:108) θ c π − (cid:109) , (5.3) r c e i θ c = √ π n b T H (1)0 ( kb ) e − i π , r c > , − π ≤ θ c ≤ π, (5.4) r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... and for any integer p . We use the superscript “o” to distinguish these wavenumbers for monopolescatterers from others. Even though (5.3) was deduced for large integer p , it gives remarkablyagreement with numerically calculated wavenumbers, except for the two lowest attenuatingwavenumbers, as shown in Figure 4. In the figure we denoted S o ±∗ as the effective wavenumberthat can be calculated by low volume fraction expansions [11,14]. ■■■■■■■■■ ■ ■ ■ ■ ■ ■ ■ ■ ■▲▲▲▲▲▲▲▲▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ S * o S - S - S + S + S + ■ Numeric ▲ Asymptotic - -
20 20 40 60 Re S S Figure 4.
Comparison of the asymptotic formula (5.3) , which predicts an infinite number of effective wavenumbers, withnumerical solutions for the effective wavenumbers (5.1) . The parameters used are given by (5.2) , with their definitionsexplained in (2.4–2.9). Here we chose b = 1 . , so the non-dimensional wavenumbers bS are the same as shown. Theasymptotic formula is surprisingly accurate except for the two lowest attenuating wavenumbers. The wavenumber S o ∗ canbe calculated by using low volume fraction expansions [11]. For multipole scatterers , where both n and m could potentially range from −∞ to ∞ in (3.15),we can also calculate an infinite sequence of effective wavenumbers. To show this explicitly, weconsider the limit of large bk , with | k | ∼ | S | . In the opposite limit bk (cid:28) , the Rayleigh limit, onlyone effective wavenumber is required [44,45].At leading order, the asymptotic solution of (A 11) leads to the effective wavenumbers: bS k ± p = σ ± p + i log | σ ± p − a | (cid:113) a | σ ± p | r c , (5.5) σ + p = θ c + a + 2 πp for p > (cid:108) − θ c + a π (cid:109) ,σ − p = θ c + a − π − πp for p > (cid:108) θ c + a π − (cid:109) , (5.6) r c e i θ c = − i n b ∞ (cid:88) n = −∞ T n , r c > and − π ≤ θ c ≤ π, (5.7)for integer p . This confirms that there are an infinite number of effective wavenumbers forlarge scatterers, i.e. bk (cid:29) . The distribution of these wavenumbers is similar to the monopolewavenumbers shown in Figure 4.These asymptotic formulas (5.3) and (5.5) demonstrate the existence of multiple effectivewaves in the limit of small (monopole and Dirichlet) scatterers (5.3) and large scatterers (5.5).However, neither of these formulas, nor the low volume fraction expansions of thewavenumber [11], are able to accurately estimate the low attenuating backward travellingeffective wavenumber such as S shown in Figure 3 (in this case not related to the S o ± and S k ± r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... given above). There is currently no way to analytically estimate these types of wavenumbers, eventhough they are necessary to accurately calculate transmission due to their small attenuation. Theonly approach it seems is to numerically solve (3.30).
6. Numerical results
Here we present numerical results for monopole scatterers, as these have explicit expressionsfor reflection (4.5) and the transmitted wave (4.4) (or more accurately the average scatteringcoefficients). We compare our analytic solution with a classical method that assumes only oneeffective wavenumber [11,13], and the Matching Method [18], recently proposed by the authors.It should be noted that all of these approaches aim to solve the same equation (2.14).Note that for monopole scatterers, using only one effective wavenumber s can, in somecases, lead to accurate results. However, for multipole scatterers (a more common scenariopractically) this is rarely the case because, as shown by Figure 3, there can be at least two effectivewavenumbers with low attenuation, and therefore both are needed to obtain accurate results.For the numerical examples we use the parameters T = − J ( ka o )H (1)0 ( ka o ) , b = 1 . , a o = 0 . , θ inc = π , φ = 30% , (6.1)which implies that the number fraction n ≈ . per unit area. When we choose to fix thewavenumber, as we do for Figure 5 and Figure 6, we use bk = 1 . . This leads to a wavelength( π/k ) which is roughly six times larger than the particle diameter. If the particle was, say, morethan a hundred times smaller than the wavelength, then only one effective wavenumber in thesum (4.4) would be necessary to accurately calculate A ( X ) . . . . . . . . x/b | A e i s x | (cid:12)(cid:12) P p =1 A p e i s p x (cid:12)(cid:12) Matching methodLow vol. frac. .
00 0 .
02 0 . . . . . x/b Figure 5.
Compares the absolute value of the average field A ( x ) calculated by different methods. The field A ( x ) isclosely related to the average transmitted wave [13]. The non-dimensional wavenumber kb = 1 . , the other parametersare given by (6.1) , with their definitions explained by (2.4–2.9). Using the Wiener-Hopf solution (4.4) , we approximate A ( x ) by using either effective wavenumbers s , s , . . . , s , or just 1 effective wavenumber s . The MatchingMethod also accounts for multiple effective wavenumbers, and is described in [18]. The low volume fraction methodassumes a low volume fraction expansion for just one effective wavenumber [11]. The small graph on the right is amagnification of the region around x = 0 . Close to the boundary x = 0 , both A e i s x and the low volume fractionmethod are inaccurate, which would potentially lead to inaccurate predictions for transmission and reflection. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... To start we compare the average scattering coefficient A ( x ) calculated by the Wiener-Hopfsolution (4.4) with other methods in Figure 5. The most accurate of these other methods isthe Matching Method [18,46], and it closely agrees with the Wiener-Hopf solution when using effective wavenumbers. The exception is the region close to the boundary x = 0 , where theWiener-Hopf solution experiences a rapid transition. The low volume fraction method is themost commonly used in the literature: it assumes a small particle volume fraction and just oneeffective wavenumber [11,13]. One significant conclusion we can draw from Figure 5 is that boththe low volume fraction method and A e i s x are inaccurate near the boundary x = 0 . This meansthat both of these methods lead to inaccurate reflection coefficients.In general, the Wiener-Hopf method does not lead to an explicit formula for the reflectioncoefficient (3.35), because we do not have an exact factorisation (3.23) for any truncatedsquare matrices. However, there are methods [13,18,20,29,30,34] to calculate A n ( x ) , from whichwe can obtain the reflection coefficient (2.12). The method [18] also accounts for multipleeffective wavenumbers. So one important question is: when using (2.12), how many effectivewavenumbers do we need to obtain an accurate reflection coefficient? . . . . . . − − − P | R P − R || R | Figure 6.
Demonstrates, with a log-log graph, how increasing the number of effective waves P leads to a more accuratereflection coefficient R P , when using (4.6) . The non-dimensional wavenumber is kb = 1 . , and the other parametersused are given by (6.1) , with their definitions explained by (2.4–2.9). Here R is the reflection coefficient given by (4.5) .The error | R P − R | continuously drops as P increases because of the rapid transition that occurs to A ( x ) near theboundary x = 0 , see Figure 5. However, methods such as the Matching Method [18] are able to accurately calculate thereflection coefficient without taking into account this rapid transition. In Figure 6 we show how increasing the number of effective waves P reduces the errorbetween R P (4.6) and R (4.5). To calculate a highly accurate reflection coefficient R , we coulduse either (4.5) or the Matching Method [18,46], as both give approximately the same R .Now we ask: how does the reflection coefficient (4.6), deduced via the Wiener-Hopf technique,compare with other methods across a broader range of wavenumbers. The result is shown inFigure 7, where R O is a low volume fraction expansion of just one effective wavenumber [13].The reflection coefficient R M is calculated from the Matching Method [18,46]. The general trendis clear: R O becomes more inaccurate as we increase the background wavenumber kb . On theother hand both R M and R agree closely over all k .One result to note is the “instability“ exhibited by the Wiener-Hopf solution near the boundary x = 0 , see Figure 5. This instability occurs because we represented A ( x ) as a superpositionof truncated waves, which is only accurate as long as the discarded terms are small. So, for For the low volume fraction method we used a small volume fraction expansion for the wavenumber, but we numericallyevaluated the wave amplitude. This is because the alternative, a small volume fraction expansion of the wave amplitude, ledto poor results. We use the reflection coefficient [13, equation (39)], rather than the explicit low volume fraction expansion [13, equation (40– 41) ]. This is because using equations (40 – 41) led to roughly double the error we show. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... . . . . kb | R − R O || R | . . . . × − × − × − × − × − × − kb | R − R M || R | Figure 7.
Compares different methods for calculating the reflection coefficient when varying the non-dimensionalwavenumber kb . The other parameters used are given by (6.1) , with their definitions explained in (2.4–2.9). Here R isgiven by the Wiener-Hopf solution (4.5) , R O uses a low volume fraction expansion of just one effective wavenumber [13],and R M is calculated from the Matching Method [18]. a truncation number P , we can expect the instability to occur when e i s P x is not small, i.e. x ≈ / Im s p . However, this instability does not affect the accuracy of the reflection coefficient (4.5)deduced by the Wiener-Hopf technique, as demonstrated by close agreement with the MatchingMethod in Figure 7.
7. Conclusion and Next Steps
The major result of this paper is to prove that the ensemble-averaged field in random particulatematerials consists of a superposition of waves, with complex effective wavenumbers, for onefixed incident wavenumber. These effective wavenumbers are governed by the dispersionequation (5.1) and are independent of the angle of incidence θ inc . We showed asymptoticallyin section 5 that this has an infinite number of solutions, and hence there are an infinitenumber of effective wavenumbers. The Wiener-Hopf technique also provides a simple andelegant expression for the reflection coefficient (3.35), whose form can be used to guide and assessmethods to characterise microstructure [47,48].To numerically implement the Wiener-Hopf technique, we considered particles that scatteronly in their monopole mode in section 6. There we saw that when close to the interface of thehalf-space, a large number of effective wavenumbers were necessary to reach accurate agreementwith an alternative method from the literature, the Matching Method as introduced by the authorsin [18]. To obtain a constructive method via the Wiener-Hopf technique for general scatterers,and not just monopole scatterers, will require the factorisation of a matrix-function [31], which ischallenging. For these reasons the Matching Method [18] is presently more effective than using theWiener-Hopf technique. However, there is ongoing work to use approximate methods [33,34,49]which exploit the symmetry and properties of the matrix (3.15).Moving forwards, this paper together with [18], establish accurate and robust solutions to thegoverning equation (2.14). These same methods can now be translated to three spatial dimensionsand vectorial waves (e.g. elasticity and electromagnetics), with much of the groundwork alreadyavailable [12,16,40]. Some clear challenges, that can now be addressed, are to verify the accuracyof the statistical assumptions used to deduce (2.14). These include the hole-correction and thequasicrystalline approximations. As these are now the only assumptions used, we could comparethe solution of (2.14) with multipole methods [50,51] in order to investigate their accuracy andlimits of validity. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... A. The Fourier transformed kernel ˆ ψ n ( s ) Here we calculate the Fourier transform (3.3) of ψ n ( X ) (2.13). To do so, it is simpler to use F n ( X, Y ) = ( − n H (1) n ( kR ) e i nΘ . (A 1)Note that both F n ( X, Y ) and e i ( sX + Y k sin θ inc ) satisfy wave equations, with ∇ F n ( X, Y ) = − k F n ( X, Y ) and ∇ e i ( sX + Y k sin θ inc ) = − S e i ( sX + Y k sin θ inc ) , where we used (2.13) for the first equation and (3.19) for the second equation. This means thatwe can use Green’s second identity to obtain ( k − S ) ˆ B e i ( sX + Y k sin θ inc ) F n ( X, Y )d X d Y = ˆ ∂ B (cid:34) ∂ e i ( sX + Y K sin θ inc ) ∂ n F n ( X, Y ) − e i ( sX + KY sin θ inc ) ∂ F n ( X, Y ) ∂ n (cid:35) d z, (A 2)for any area B in which the integrand is analytic, where n is the outwards pointing unit normaland d z is a differential length along the boundary ∂ B . To calculate ˆ ψ n ( s ) , we take the region B tobe defined by R ≥ b , with ( R, Θ ) being the polar coordinates of ( X, Y ) , in which case the integralover B converges because as R → ∞ we have that | e i ( sX + Y k sin θ inc ) F n ( X, Y ) | ∼ | e i sR cos Θ e i kR (1+sin Θ sin θ inc ) | (cid:112) π | k | R/ ≤ | e − R (Im k (1 −| sin θ inc | ) −| Im ( s ) | ) | (cid:112) π | k | R/ → , (A 3)exponentially fast when | Im ( s ) | < Im k (1 − | sin θ inc | ) . Under this restriction, and by assuming S (cid:54) = ± k , (A 2) then leads to ˆ ψ n ( s ) = ˆ R ≥ b e i sX + i kY sin θ inc F n ( X, Y )d X d Y = I n ( b ) k − S , (A 4)by using sX + Y k sin θ inc = RS cos( θ − θ S ) from (3.19) and I n ( R ) = ˆ π − ∂ e i SR cos( Θ − θ S ) ∂R F n ( k X ) + e i SR cos( Θ − θ S ) ∂ F n ( k X ) ∂R R d Θ = ( − n ˆ π e i SR cos( Θ − θ S ) e i nΘ (cid:2) k H (1) (cid:48) n ( kR ) − i S cos( Θ − θ S )H (1) n ( kR ) (cid:3) R d θ = ( − n ˆ π ∞ (cid:88) m = −∞ i m J m ( SR ) e i m ( Θ − θ S ) (cid:104) k e i nΘ H (1) (cid:48) n ( kR ) − i S e i ( n +1) Θ − i θ S + e i ( n − Θ + i θ S )H (1) n ( kR ) (cid:105) R d Θ = 2 π ( − i ) n R e i nθ S (cid:104) k J n ( SR )H (1) (cid:48) n ( kR ) − S J (cid:48) n ( SR )H (1) n ( kR ) (cid:105) , (A 5)where J n is the Bessel function of the first kind, and we used the Jacobi-Anger expansion one i SR cos( Θ − θ S ) , integrated over Θ and used the identity J n − ( SR ) − J n +1 ( SR ) = 2J (cid:48) n ( SR ) . Insummary ˆ ψ n ( s ) = 2 π ( − i ) n e i nθ S α − s N n ( bS ) , (A 6)when the condition (3.10) is satisfied, with N n given by (3.17).Below we establish some useful properties for ˆ ψ n ( s ) . In particular, we show that ˆ ψ n ( s ) has nobranch-points. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... The function N n ( bS ) , for integer values of n , can be expanded around S = 0 as N n ( bS ) = S | n | ∞ (cid:88) m =0 c m | n | S m , (A 7)where the c m | n | are some constants that depend on m and | n | , and the radius of convergence ofthe series above is infinite. Using (3.19) we can writee i nθ S = e i sgn( n ) | n | θ S = (cos θ S + sgn( n ) i sin θ S ) | n | = ( s + sgn( n ) i k sin θ inc ) | n | S −| n | . (A 8)Substituting (A 7) and (A 8) in (A 6) results in ˆ ψ n ( s ) = 2 π ( − i ) n α − s ( s + sgn( n ) i k sin θ inc ) | n | ∞ (cid:88) m =0 c m | n | S m . (A 9)Because S = s + k sin θ inc , we can immediately see from the above that ˆ ψ n ( s ) has no branch-points. Additionally we can establish the properties: ˆ ψ n ( s ) = ˆ ψ − n ( − s ) = ˆ ψ n ( − s ) e i nθ S ( − n . (A 10) B. Asymptotic location of the wavenumbers
Here we explicitly calculate a sequence of effective wavenumbers S p , assuming | S p | large andincreasing with p , and Im S p > , that asymptotically satisfy (5.1). A key step is to approximatethe terms appearing in (3.16), such as J n ( bS ) ∼ e i π + i nπ − i bS √ πbS and J (cid:48) n ( bS ) ∼ e − i π + i nπ − i bS √ πbS , (A 1)for large | bS | , where the terms e i bS are discarded as Im bS → ∞ . Monopole scatterers
The simplest case is for monopole scatterers, where n = m = 0 in (3.16),and the effective wavenumber S satisfies b det G = ( bS ) − ( bk ) + 2 π n b T N ( bS ) ∼ ( bS ) − c √ bS e − i bS = 0 , (A 2)where c = √ π n b T H (1)0 ( kb ) e − i π . Here we used (A 1), and ignored terms which are algebraicallysmaller than bS . To find the root of the above we substitute bS = x + i log y, (A 3)where x and y are real, and | x | and y are large with y > . This leads to ( x + i log y ) / − c e − i x y = 0 . (A 4)For the logarithm and square root we use the typical branch cut ( −∞ , and take positive valuesof the functions for positive arguments. For the above to be satisfied to leading order then x / ∼ y , which reduces the above equation to x / ∼ r c e i ( θ c − x ) y, (A 5)where we substituted c = r c e i θ c , for real scalars r c and θ c . Equating the real and imaginary partsof the above leads to x ∼ θ c + 2 πp and y ∼ r c ( θ c + 2 πp ) / for p > − θ c π (A 6) x ∼ θ c − π πp and y ∼ r c ( − θ c + 3 π − πp ) / for p < − θ c π , (A 7)for integers p . From this we can identify that, at leading order, the effective wavenumbers aregiven by (5.3). r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Multipole scatterers:
With the same method used above, we can also demonstrate the existenceof multiple effective wavenumbers for n, m = − M, − M + 1 , · · · , M in (3.16). To show thisexplicitly, we consider bk to be the same order as bS , that is | k | ∼ | S | .By considering bk large, we can approximate H (1) n ( bk ) ∼ e i ( bk − π − nπ ) (cid:114) πbk and H (1) (cid:48) n ( bk ) ∼ e i ( bk + π − nπ ) (cid:114) πbk , (A 8)combining this with (A 1) and considering | k | ∼ | S | , the term (3.16) at leading order becomes b G mn = d δ mn + c T m , (A 9)where d = ( bS ) − ( bk ) , and c = 2 n b i ( k + S ) √ kS e i b ( k − S ) . By simple rearrangement of the determinant we find that det( b G ) = d M d + c M (cid:88) m = − M T m . (A 10)Note that d (cid:54) = 0 , i.e. S (cid:54) = ± k , was necessary to reach the condition (3.10), which was used tocalculate the Fourier transforms (3.13). Taking this into consideration, and taking the limit M →∞ , the effective wavenumbers S must satisfy d + c M (cid:88) m = − M T m = 0 = ⇒ bS − bk = − n i b ∞ (cid:88) m = −∞ T m e i b ( k − S ) b √ kS . (A 11)Using an asymptotic expansion analogous to (A 3), the above leads to the effectivewavenumbers (5.5). C. Equivalent determinants
For any square matrices A and B , and scalar c , if A nm = B nm c n − m (not employing thesummation convention), then det A = det B . (A 1)This follows simply by defining the diagonal matrix C nm = δ nm c n , which leads to A = CBC − ,and det( CBC − ) = det C det B det C − = det B .Data Accessibility. We provide the code to generate all the graphs in [46].
Authors’ Contributions.
I.D.A and A.L.G. conceived of the study. A.L.G. drafted the manuscript. A.L.G.,I.D.A, and W.J.P edited the manuscript and gave final approval for publication.
Competing Interests.
We have no competing interests.
Funding.
This work was funded by EPSRC (EP/M026205/1,EP/L018039/1,EP/S019804/1) including viathe Isaac Newton Institute (EP/K032208/1,EP/R014604/1) and partial support from the UKAN grant EPSRC(EP/R005001/1).
References
1. Pinfield VJ, Challis RE. 2013 Emergence of the coherent reflected field for a single realisationof spherical scatterer locations in a solid matrix.
Journal of Physics: Conference Series , 012009.2. Mishchenko MI, Travis LD, Lacis AA. 2006
Multiple Scattering of Light by Particles: RadiativeTransfer and Coherent Backscattering . Cambridge University Press. The determinant of b G equals the product of its eigenvalues. The eigenvector ( T − M , · · · , T M ) T gives the eigenvalue d + c (cid:80) m T m , while all other eigenvalues equal d . r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
3. Mishchenko MI, Dlugach JM, Yurkin MA, Bi L, Cairns B, Liu L, Panetta RL, Travis LD, YangP, Zakharova NT. 2016 First-principles modeling of electromagnetic scattering by discrete anddiscretely heterogeneous random media.
Physics Reports , 1–75. arXiv: 1605.06452.4. Tsang L, Ishimaru A. 1987 Radiative Wave Equations for Vector Electromagnetic Propagationin Dense Nontenuous Media.
Journal of Electromagnetic Waves and Applications , 59–72.5. Tsang L, Chen CT, Chang ATC, Guo J, Ding KH. 2000 Dense media radiative transfer theorybased on quasicrystalline approximation with applications to passive microwave remotesensing of snow. Radio Science , 731–749.6. Foldy LL. 1945 The multiple scattering of waves. I. General theory of isotropic scattering byrandomly distributed scatterers. Physical Review , 107.7. Lax M. 1951 Multiple Scattering of Waves. Reviews of Modern Physics , 287–310.8. Fikioris JG, Waterman PC. 1964 Multiple Scattering of Waves. II. “Hole Corrections” in theScalar Case. Journal of Mathematical Physics , 1413–1420.9. Tsang L, Kong JA. 1983 Scattering of electromagnetic waves from a half space of denselydistributed dielectric scatterers. Radio Science , 1260–1272.10. Tishkovets VP, Petrova EV, Mishchenko MI. 2011 Scattering of electromagnetic waves byensembles of particles and discrete random media. Journal of Quantitative Spectroscopy andRadiative Transfer , 2095–2127.11. Linton CM, Martin PA. 2005 Multiple scattering by random configurations of circularcylinders: Second-order corrections for the effective wavenumber.
The Journal of the AcousticalSociety of America , 3413.12. Linton CM, Martin PA. 2006 Multiple Scattering by Multiple Spheres: A New Proof of theLloyd–Berry Formula for the Effective Wavenumber.
SIAM Journal on Applied Mathematics ,1649–1668.13. Martin PA. 2011 Multiple scattering by random configurations of circular cylinders:Reflection, transmission, and effective interface conditions. The Journal of the Acoustical Societyof America , 1685–1695.14. Norris AN, Conoir JM. 2011 Multiple scattering by cylinders immersed in fluid: High orderapproximations for the effective wavenumbers.
The Journal of the Acoustical Society of America , 104–113.15. Gower AL, Smith MJA, Parnell WJ, Abrahams ID. 2018 Reflection from a multi-speciesmaterial and its transmitted effective wavenumber.
Proc. R. Soc. A , 20170864.16. Conoir JM, Norris AN. 2010 Effective wavenumbers and reflection coefficients for an elasticmedium containing random configurations of cylindrical scatterers.
Wave Motion , 183–197.17. Norris AN, Luppé F, Conoir JM. 2012 Effective wave numbers for thermo-viscoelastic mediacontaining random configurations of spherical scatterers. The Journal of the Acoustical Society ofAmerica , 1113–1120.18. Gower AL, Parnell WJ, Abrahams ID. 2018 Multiple Waves Propagate in Random ParticulateMaterials. arXiv:1810.10816 [physics] . arXiv: 1810.10816.19. Crighton DG, Dowling A, Ffowcs Williams J, Heckl M, Leppington F. 1992
Modern Methods inAnalytical Acoustics . Springer-Verlag.20. Lawrie JB, Abrahams ID. 2007 A brief historical perspective of the Wiener–Hopf technique.
Journal of Engineering Mathematics , 351–358.21. Noble B. 1988 Methods Based on the Wiener-Hopf Technique. vol. 67. American MathematicalSociety 2nd unexpurgated edition edition.22. Martin PA, Abrahams ID, Parnell WJ. 2015 One-dimensional reflection by a semi-infiniteperiodic row of scatterers.
Wave Motion , 1–12.23. Norris A, Wickham GR. 1995 Acoustic diffraction from the junction of two flat plates. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences ,631–655.24. Haslinger SG, Movchan NV, Movchan AB, Jones IS, Craster RV. 2016 Controlling flexuralwaves in semi-infinite platonic crystals. arXiv:1609.02787 [physics] . arXiv: 1609.02787.25. Tymis N, Thompson I. 2014 Scattering by a semi-infinite lattice and the excitation of Blochwaves.
The Quarterly Journal of Mechanics and Applied Mathematics , 469–503.26. Haslinger SG, Jones IS, Movchan NV, Movchan AB. 2017 Semi-infinite herringbonewaveguides in elastic plates. arXiv:1712.01827 [physics] . arXiv: 1712.01827.27. Albani M, Capolino F. 2011 Wave dynamics by a plane wave on a half-space metamaterialmade of plasmonic nanospheres: a discrete Wiener–Hopf formulation. JOSA B , 2174–2185. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
28. Abrahams ID. 1997 On the Solution of Wiener–Hopf Problems Involving NoncommutativeMatrix Kernel Decompositions.
SIAM Journal on Applied Mathematics , 541–567.29. Abrahams ID. 1996 Radiation and scattering of waves on an elastic half-space; A non-commutative matrix Wiener-Hopf problem. Journal of the Mechanics and Physics of Solids ,2125–2154.30. Abrahams ID, Wickham GR. 1990 General Wiener–Hopf Factorization of Matrix Kernels withExponential Phase Factors. SIAM Journal on Applied Mathematics , 819–838.31. Rogosin S, Mishuris G. 2016 Constructive methods for factorization of matrix-functions. IMAJournal of Applied Mathematics , 365–391.32. Kisil A. 2018 An Iterative Wiener–Hopf Method for Triangular Matrix Functions withExponential Factors. SIAM Journal on Applied Mathematics , 45–62.33. Abrahams ID. 2000 The application of Padé approximants to Wiener-Hopf factorization. IMAJournal of Applied Mathematics , 257–281.34. Abrahams ID. 1987 Scattering of sound by two parallel semi-infinite screens. Wave Motion ,289–300.35. Ganesh M, Hawkins SC. 2010 A far-field based T-matrix method for two dimensional obstaclescattering. ANZIAM Journal , 215–230.36. Ganesh M, Hawkins SC. 2017 Algorithm 975: TMATROM—A T-Matrix Reduced Order ModelSoftware. ACM Trans. Math. Softw. , 9:1–9:18.37. Mishchenko MI. 1991 Light scattering by randomly oriented axially symmetric particles. JOSAA , 871–882.38. Mishchenko MI. 1993 Light scattering by size–shape distributions of randomly orientedaxially symmetric particles of a size comparable to a wavelength. Applied Optics , 4652–4666.39. Waterman PC. 1971 Symmetry, Unitarity, and Geometry in Electromagnetic Scattering. Physical Review D , 825–839.40. Kristensson G. 2015 Coherent scattering by a collection of randomly located obstacles – Analternative integral equation formulation. Journal of Quantitative Spectroscopy and RadiativeTransfer , 97–108.41. Bleistein N, Handelsman RA. 1986
Asymptotic expansions of integrals . Courier Corporation.42. Gokhberg IC, Krein MG. 1960 Systems of integral equations on the half-line with kernelsdepending on the difference of the arguments. , 217 – 287.43. Martin PA, Maurel A. 2008 Multiple scattering by random configurations of circular cylinders:Weak scattering without closure assumptions. Wave Motion , 865–880.44. Parnell WJ, Abrahams ID. 2010 Multiple point scattering to determine the effectivewavenumber and effective material properties of an inhomogeneous slab. Waves in Randomand Complex Media , 678–701.45. Parnell WJ, Abrahams ID, Brazier-Smith PR. 2010 Effective Properties of a CompositeHalf-Space: Exploring the Relationship Between Homogenization and Multiple-ScatteringTheories. The Quarterly Journal of Mechanics and Applied Mathematics , 145–175.46. Gower AL EffectiveWaves.jl: A package to calculate ensemble averaged waves inheterogeneous materials. https://github.com/arturgower/EffectiveWaves.jl/tree/v0.2.1 .47. Roncen R, Fellah ZEA, Simon F, Piot E, Fellah M, Ogam E, Depollier C. 2018 Bayesian inferencefor the ultrasonic characterization of rigid porous materials using reflected waves by the firstinterface. The Journal of the Acoustical Society of America , 210–221.48. Gower AL, Gower RM, Deakin J, Parnell WJ, Abrahams ID. 2018 Characterising particulaterandom media from near-surface backscattering: A machine learning approach to predictparticle size and concentration.
EPL (Europhysics Letters) , 54001.49. Veitch BH, Abrahams ID. 2007 On the commutative factorization of nxn matrix Wiener–Hopfkernels with distinct eigenvalues.
Proceedings of the Royal Society A: Mathematical, Physical andEngineering Sciences , 613–639.50. Martin PA. 2006
Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles .Cambridge: Cambridge University Press.51. Gower AL, Deakin J. 2019 MultipleScattering.jl: A Julia library for simulating, processing,and plotting multiple scattering of acoustic waves. https://github.com/jondea/MultipleScattering.jlhttps://github.com/jondea/MultipleScattering.jl