A Proof that Multiple Waves Propagate in Ensemble-Averaged Particulate Materials

TL;DR

This paper proves that multiple effective wavenumbers exist in ensemble-averaged particulate materials, showing that wave propagation involves several waves and that accurate transmission calculations require considering many of these waves.

Contribution

It provides a rigorous proof of the non-uniqueness of effective wavenumbers and analyzes their significance in wave propagation through particulate media.

Findings

Multiple effective wavenumbers exist in such materials.

Only a few dominate the wave field in most regimes.

Accurate reflection and transmission calculations need many effective waves.

Abstract

Effective medium theory aims to describe a complex inhomogeneous material in terms of a few important macroscopic parameters. To characterise wave propagation through an inhomogeneous material, the most crucial parameter is the effective wavenumber. For this reason, there are many published studies on how to calculate a single effective wavenumber. Here we present a proof that there does not exist a unique effective wavenumber; instead, there are an infinite number of such (complex) wavenumbers. We show that in most parameter regimes only a small number of these effective wavenumbers make a significant contribution to the wave field. However, to accurately calculate the reflection and transmission coefficients, a large number of the (highly attenuating) effective waves is required. For clarity, we present results for scalar (acoustic) waves for a two-dimensional material filled (over a half space) with randomly distributed circular cylindrical inclusions. We calculate the effective medium by ensemble averaging over all possible inhomogeneities. The proof is based on the application of the Wiener-Hopf technique and makes no assumption on the wavelength, particle boundary conditions/size, or volume fraction. This technique provides a simple formula for the reflection coefficient, which can be explicitly evaluated for monopole scatterers. We compare results with an alternative numerical matching method.

Equations

$${\nabla }^{2}u+k^{2}u=0,\text{(in the background material)},$$

$${\nabla }^{2}u+k_o^{2}u=0,\text{(inside a particle)},$$

$$u(x,y)=u_{\mathrm{inc}}(x,y)+\sum _j{u}_j(x,y).$$

$$u_{\mathrm{inc}}(x,y)={\text{e}}^{\text{i}(\alpha x+\beta y)},\text{with}(\alpha ,\beta )=(k\cos {\theta }_{\mathrm{inc}},k\sin {\theta }_{\mathrm{inc}}),$$

$$\mathrm{Re}k>0\text{and}\mathrm{Im}k>0.$$

$$b=\text{the minimum distance between particle centres},$$

$$\mathfrak{n}=\text{number of particles per unit area},$$

$$T_n=\text{the coefficients of the particle’s T-matrix},$$

$$\varphi =\frac{\pi \mathfrak{n}{b}^2}{4}=\text{particle area fraction}.$$

$$⟨u_1(x_1+X,y_1+Y){⟩}_{(x_1,y_1)}=\sum _{n=-\infty }^{\infty }{A}_n(x_1){\text{e}}^{\text{i}\beta y_1}{\mathrm{H}}_n^{(1)}(kR){\text{e}}^{\text{i}n\Theta },$$

$$⟨u(x,y)⟩=u_{\mathrm{inc}}(x,y)+\mathfrak{R}{\text{e}}^{-\text{i}\alpha x+\text{i}\beta y}\text{for}x<-b,$$

$$\mathfrak{R}={\text{e}}^{\text{i}\alpha x}\mathfrak{n}\sum _{n=-\infty }^{\infty }{\int }_0^{\infty }{A}_n(x_1){\psi }_n(x_1-x)d{x}_1,$$

$${\psi }_n(X)={\int }_{Y^2>b^2-X^2}{\text{e}}^{\text{i}\beta Y}(-1{)}^n{\mathrm{H}}_n^{(1)}(kR){\text{e}}^{\text{i}n\Theta }dy.$$

$$\begin{gathered} \mathfrak{n}{T}_m\sum _{n=-\infty }^{\infty }{\int }_0^{\infty }{A}_n(x_2){\psi }_{n-m}(x_2-x_1)\mathrm{d}{x}_2 \\ =A_m(x_1)-{\text{e}}^{\text{i}\alpha x_1}{T}_m{\text{e}}^{\text{i}m(\pi /2-{\theta }_{\mathrm{inc}})},\text{for}{x}_1\ge 0, \end{gathered}$$

$$\begin{gathered} \mathfrak{n}{T}_m\sum _{n=-\infty }^{\infty }{\int }_0^{\infty }{A}_n(x_2){\psi }_{n-m}(x_2-x_1)\mathrm{d}{x}_2 \\ =\{\begin{array}{ll}{A}_m(x_1)-{\text{e}}^{\text{i}\alpha x_1}{T}_m{\text{e}}^{\text{i}m(\pi /2-{\theta }_{\mathrm{inc}})},& x_1\ge 0,\\ D_m(x_1),& x_1<0,\end{array} \end{gathered}$$

$$D_0(x)=T_0\mathfrak{R}{\text{e}}^{-\text{i}\alpha x}.$$

$$\hat{f}(s)={\int }_{-\infty }^{\infty }f(x){\text{e}}^{\text{i}sx}dx\text{with}f(x)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\hat{f}(s){\text{e}}^{-\text{i}sx}dx,$$

$${\hat{A}}_n^{+}(s)={\int }_0^{\infty }{A}_n(x){\text{e}}^{\text{i}sx}dx,{\hat{D}}_n^{-}(s)={\int }_{-\infty }^0{D}_n(x){\text{e}}^{\text{i}sx}dx.$$

$$|A_n(x)|<{\text{e}}^{-xc}\text{for}x\to \infty ,$$

$$|D_n(x)|<{\text{e}}^{xc}\text{for}x\to -\infty ,$$

$$|\mathrm{Im}s|<c.$$

$${\hat{A}}_n^{+}(s)=\mathcal{O}(|s{|}^{-1})\text{and}{\hat{D}}_n^{-}(s)=\mathcal{O}(|s{|}^{-1})\text{for}|s|\to \infty ,$$

$$\begin{gathered} {\int }_0^{\infty }{A}_n(x_2){\int }_{-\infty }^{\infty }{\psi }_{n-m}(x_1-x_2){\text{e}}^{\text{i}s{x}_1}\mathrm{d}{x}_1\mathrm{d}{x}_2 \\ ={\hat{A}}_n^{+}(s){\hat{\psi }}_{n-m}(s), \end{gathered}$$

$$|\mathrm{Im}s|<(1-|\sin {\theta }_{\mathrm{inc}}|)\mathrm{Im}k,$$

$$\begin{gathered} {\int }_{-\infty }^0{D}_m(x_1){\text{e}}^{\text{i}s{x}_1}d{x}_1+{\int }_0^{\infty }{A}_m(x_1){\text{e}}^{\text{i}s{x}_1}\mathrm{d}{x}_1 \\ -{\text{e}}^{\text{i}m(\pi /2-{\theta }_{\mathrm{inc}})}{T}_m{\int }_0^{\infty }{\text{e}}^{\text{i}{x}_1(s+\alpha )}\mathrm{d}{x}_1={\hat{D}}_m^{-}(s)+{\hat{A}}_m^{+}(s)-T_m\frac{\text{i}{\text{e}}^{\text{i}m(\pi /2-{\theta }_{\mathrm{inc}})}}{(s+\alpha {)}^{+}}, \end{gathered}$$

$$|\mathrm{Im}s|<ϵ,\text{where}ϵ=\min \{c,(1-|\sin {\theta }_{\mathrm{inc}}|)\mathrm{Im}k,\mathrm{Im}\alpha \}.$$

$$\frac{\mathbf{\Psi }(s){\hat{\mathbf{A}}}^{+}(s)}{s^2-{\alpha }^2}=-{\hat{\mathbf{D}}}^{-}(s)+\frac{\mathbf{B}}{s+\alpha },$$

$$B_m=\text{i}{T}_m{\text{e}}^{\text{i}m(\pi /2-{\theta }_{\mathrm{inc}})},$$

$${\Psi }_{mn}(s)=G_{mn}(S)(-\text{i}{)}^{n-m}{\text{e}}^{\text{i}(n-m){\theta }_S},$$

$$G_{mn}(S)=(s^2-{\alpha }^2){\delta }_{mn}+2\pi \mathfrak{n}{T}_m{\mathrm{N}}_{n-m}(bS),$$

$${\mathrm{N}}_m(bS)=bk{\mathrm{J}}_m(bS){\mathrm{H}}_m^{(1)\prime }(bk)-bS{\mathrm{J}}_m^{\prime }(bS){\mathrm{H}}_m^{(1)}(bk).$$

$${\Psi }_{mn}(s)=(s^2-{\alpha }^2)\left[{\delta }_{mn}-\mathfrak{n}{T}_m{\hat{\psi }}_{n-m}(s)\right],$$

$$s=S\cos {\theta }_S\text{with}S\sin {\theta }_S=k\sin {\theta }_{\mathrm{inc}}.$$

$${\theta }_S=\arctan \left(\frac{k\sin {\theta }_{\mathrm{inc}}}{s}\right),S=\sqrt{s^2+(k\sin {\theta }_{\mathrm{inc}}{)}^2},$$

Full text

1 Introduction

Materials comprising particles or inclusions that are randomly distributed inside a uniform host medium occur frequently in the world around us. They occur as synthetically fabricated media and also in nature. Common examples include composites, emulsions, suspensions, complex gases, and polymers. Understanding how electromagnetic, elastic, or acoustic waves propagate through these materials is necessary in order to characterise the properties of these materials, and also to design new materials that can control wave propagation.

The wave scattered from a particulate material will be influenced by the positions and properties of all particles, which are usually unknown. However, this scattered field, averaged over space or over time, depends only on the average particle properties. Many measurement systems perform averaging over space, if the receivers or incident wavelength are large enough [1], or over time [2]. In most cases, this averaging process is the same as averaging over all possible particle configurations. Such systems are sometimes called ergodic [2, 3]. In this paper, we focus on ensemble averaged waves, satisfying the scalar wave equation in two-dimensions, reflecting from, and propagating in, a half-space particulate material. In certain scenarios, such as light scattering, it is easier to measure the average intensity of the wave. However, even in these cases, the ensemble-averaged field is often needed as a first step [4, 5].

One driving principle, often used in the literature, is that the ensemble-averaged wave itself satisfies a wave equation with a single effective wavenumber [6, 7, 8]. Reducing a inhomogeneous material, with many unknowns, down to one effective wavenumber is attractive as it greatly reduces the complexity of the problem. For this reason many papers have attempted to deduce this unique effective wavenumber from first principles in electromagnetism [9, 10, 3], acoustics [11, 12, 13, 14, 15] and elasticity [16, 17]. See [18] for a short overview of the history of this topic, including typical statistical assumptions employed within the methods, such as hole-correction and the quasi-crystalline approximation, which we also adopt here.

The assumption that the ensemble averaged wave field satisfies a wave equation, with an effective wavenumber, has never been fully justified. Here we prove that there does not exist a unique effective wavenumber but instead there are an infinite number of them. Gower et al. [18] first showed that there exist many effective wavenumbers, and provided a technique, the Matching Method, to efficiently calculate the effective wave field. In the present paper and [18], we show that for some parameter regimes, at least two effective wavenumbers are needed to obtain accurate results, when compared with numerical simulations. We also provide examples of how a single effective wave approximation leads to inaccurate results for both transmission and reflection for a halfspace filled with particles, see Figure 1.

Although the Matching Method developed in [18] gave accurate results, when compared to numerical methods and known asymptotic limits, the limitations of the method were not immediately clear. Here however we illustrate that the Matching Method is robust, because combining many effective wavenumbers is not just a good approximation, it is an analytical solution to the integral equation governing the ensemble averaged wave field. We prove this by employing the Wiener-Hopf technique and then, for clarity, illustrate the solution for particles that scatterer only in their monopole mode. The Wiener-Hopf technique also gives a simple and elegant expression for the reflection coefficient.

The Wiener-Hopf technique is a powerful tool to solve a diverse range of wave scattering problems, see [19, Chapter 5. Wiener-Hopf Technique] and [20, 21] for an introduction. It is especially useful for semi-infinite domains [22, 23, 24, 25, 26, 27] and boundary value problems of mixed type. In this work, the Wiener-Hopf technique clearly reveals the form of the analytic solution, but to compute the solution would require an analytic factorisation of a matrix-function. To explicitly perform this factorisation is difficult [28, 29, 30, 31]. Indeed this is often the hardest aspect of employing the Wiener-Hopf technique, although there exist approximate methods for this purpose [28, 32, 33, 34]. We do not focus in this article on these analytic factorisations, as there already exists a method to compute the required solution [18]. Instead, the present work acts as proof that the Matching Method [18] faithfully reproduces the form of the analytic solution.

Figure 1 shows the main setup and result of this paper: an incident plane wave excites the half-space $x>0$ filled with ensemble-averaged particles (the blue region), which generates a reflected wave and many effective transmitted waves. The $\bm{s}_{p}$ are the transmitted wavevectors, and the smaller the length of the vector, the faster that effective wave attenuates as it propagates further into the material.

The paper begins by summarising the equations that govern ensemble averaged waves in two-dimensions in section 2. Following this, in section 3 we apply the Wiener-Hopf technique to the governing integral equation and deduce that the solution is a superposition of plane waves, each with a different effective wavenumber. A simple expression for the reflection coefficient is also derived. In section 4 we specialise the results for particles that scatter only in the monopole mode, which leads to a closed form analytic solution.

The dispersion relation (44), derived in section 3, admits an infinite number of solutions, the effective wavenumbers. In section 5, we deduce asymptotic forms for the effective wavenumbers in both a low and high frequency limit. In section 6 we compare numerical results for monopole scatterers, using the Wiener-Hopf technique, with classical methods that assume only one effective wavenumber [11, 13], and the Matching Method introduced in [18]. In general, when comparing predicted reflection coefficients, the Wiener-Hopf and Matching Method agree well, whereas the classical single-effective-wavenumber method can disagree by anywhere up to $20\%$. These results are discussed in section 7 together with anticipated future steps.

2 Waves in ensemble averaged particles

Consider a region filled with particles or inclusions that are uniformly distributed. The field $u$ is governed by the scalar wave equations:

[TABLE]

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. For a distribution of particles, or multi-species, see [15].

Our goal is to calculate the ensemble average field $\langle u(x,y)\rangle$, that is, the field averaged over all possible particle positions and orientations. For clarity, and ease of exposition, we consider that the particles are equally likely to be located anywhere except that they cannot overlap (this is often called the hole correction assumption). We also assume the quasi-crystalline approximation; for details on this, and for further details on deducing the results in this section, see [11, 15, 18].

By splitting the total steady wave field $u(x,y)$ into a sum of the incident wave $u_{\mathrm{inc}}(x,y)$ and waves scattered by each particle, the $j$th scattered wave being $u_{j}(x,y)$, we can write:

[TABLE]

A simple and useful scenario to consider is when all particles are placed only within the half-space111The case where particles can be placed anywhere in the plane can lead to ill-defined integrals [11]. $x>0$, which are then excited by a plane wave, with implicit time dependence $\textrm{e}^{\textrm{i}\omega t}$, incident from a homogeneous region:

[TABLE]

where we restrict the incident angle $-\frac{\pi}{2}<\theta_{\mathrm{inc}}<\frac{\pi}{2}$, as shown in Figure 1, and consider a slightly dissipative medium with

[TABLE]

This dissipation will facilitate the use of the Wiener-Hopf technique, and after reaching the solution we can take $k$ to be real222Assuming Im $k=\epsilon>0$, rather than $\epsilon\geq 0$, will facilitate calculating certain integrals that appear below. However, after reaching a solution, we can take the limit $\epsilon\to 0$ to recover the physically viable solution for Im $k=0$..

To describe the particulate medium we employ the following notation:

[TABLE]

Although the area fraction $\phi$, normally called the volume fraction, is a combination of other parameters, it is useful because it is non-dimensional. If we let $a_{o}$ be the maximum distance from the particle’s centre to its boundary, then we can set $b=\gamma a_{o}$, where $\gamma\geq 2$ so as to avoid two particles overlapping. The volume fraction that does not include the exclusion zone $\phi^{\prime}$, as used in [18, equation (4.7)], is then $\phi^{\prime}=4\phi/\gamma^{2}$.

The $T_{n}$ are the coefficients of a diagonal T-matrix [35, 36, 37, 38, 39]. The T-matrix determines how the particle scatters waves, and so depends on the particle’s shape and boundary conditions. A diagonal T-matrix can be used to represent either a radially symmetric particle, or particles averaged over their orientation, assuming the orientations have a random uniform distribution.

The results of ensemble averaging (3) from first principals are deduced in a number of references[18, 15] and so deails of this procedure are omitted here for brevity. To represent the ensemble averaged scattered wave from a particle, whose centre is fixed at $(x_{1},y_{1})$, we use

[TABLE]

for $R:=\sqrt{X^{2}+Y^{2}}>b/2$, so that $(X,Y)$ is on the outside of this particle, with $(R,\Theta)$ being the polar coordinates of $(X,Y)$, $\mathrm{H}^{(1)}_{n}$ are Hankel functions of the first kind, and $A_{n}$ is some field we want to determine333The factor $\textrm{e}^{\textrm{i}\beta y_{1}}$ appears due to the translational invariance of $\langle u_{1}(x_{1},y_{1})\rangle$ in $y_{1}$, which is a result of the material being statistically homogeneous, see [15] for details..

By choosing444We define the reflection coefficient only for $x<-b$, instead of $x<-b/2$, so that we can use $\psi_{n}$ in the formula for $\mathfrak{R}$, which will in turn facilitate calculating $\mathfrak{R}$. $x<-b$, which is outside of the region filled with particles, then taking the ensemble average on both sides of (3) results in equation (6.7) of [18], given by:

[TABLE]

which is the incident wave plus an effective reflected wave with reflection coefficient:

[TABLE]

where we assumed particles are distributed according to a uniform distribution, and the kernel $\psi_{n}$ is given by

[TABLE]

Later we show that, as expected, $\mathfrak{R}$ is independent of $x$.

The system governing $A_{m}(x)$ is given by equation (4.7) of [18]:

[TABLE]

for all integers $m$. Kristensson [40, equation (15)] presents an equivalent integral equation for electromagnetism and particles in a slab.

Our main aim is to reach an exact solution for $A_{n}(x)$ by employing the Wiener-Hopf technique to (14). We show how this also leads to simple solutions for the reflection coefficient by using (11). We acknowledge the authors of [11], as they noticed that (14) is a Wiener-Hopf integral equation; but apparently did not follow the steps indicated in the following sections555However, they were unable to solve it because, it seems, of a mistake in the integrand of equation (37) of [11]; they used $\textrm{e}^{\textrm{i}\beta Y}$ where they should have used $\cos(\beta Y)$..

3 Applying the Wiener-Hopf technique

Equation (14) is convolution integral equation with a difference kernel. This means applying a Fourier transforms can lead to elegant and simple solutions. To facilitate, we must analytically extend (14) for all $x_{1}\in\mathbb{R}$ by defining

[TABLE]

for integers $m$, where if the $A_{n}(x)$ were known for $x>0$, then the $D_{n}(x)$ would be given from the left hand-side. Note that the kernel $\psi_{n}$ defined in (13) is already analytic in the domain $\mathbb{R}$.

The field $D_{0}(x)$ is not just an abstract construct, it is closely related to the reflected wave: by directly comparing (15) with the reflection coefficient (12), for $x<-b$, we find that

[TABLE]

To solve (15) we employ the Fourier transform and its inverse, which we define as

[TABLE]

for any smooth function $f$. We then define

[TABLE]

We can determine where $\hat{A}_{n}^{+}$ and $\hat{D}_{n}^{-}$ are analytic by assuming666The solutions for $A_{n}(x)$ and $D_{n}(x)$, in the next section, show that these assumptions do hold. that

[TABLE]

for some (possibly small) positive constant $c$. This leads to $\hat{A}_{n}^{+}(s)$ being analytic for $\mathrm{Im}\,s>-c$, while $\hat{D}_{n}^{-}(s)$ is analytic for $\mathrm{Im}\,s<c$. In other words, both $\hat{A}_{n}^{+}(s)$ and $\hat{D}_{n}^{-}(s)$ are analytic in the overlapping strip

[TABLE]

To apply the Wiener-Hopf technique we also need to specify the large $s$ behaviour for both $\hat{A}_{n}^{+}(s)$ and $\hat{D}_{n}^{+}(s)$. To achieve this, we assume, on physical grounds, that $A_{n}(x)$ is bounded when $x\to 0^{+}$, and $D_{n}(x)$ is bounded when $x\to 0^{-}$. Then, it can be shown [21, 41] that

[TABLE]

in their respective half-planes of analyticity.

Applying a Fourier transform to both sides of equation (15), the left-hand side becomes

[TABLE]

in which $\hat{\psi}_{n}(s)$ is well defined (i.e. analytic) for $s$ in the strip:

[TABLE]

see appendix A for details. The right-hand side of (15) becomes

[TABLE]

where for the last step we assumed Im $(s+\alpha)>0$, which is why we use the superscript $+$ on $(s+\alpha)^{+}$. This assumption, together with (21) and (24), is satisfied if

[TABLE]

If (26) is satisfied then we can combine (23),(25) and (71), to obtain the Fourier transform of (15) in matrix form:

[TABLE]

where $\hat{\bm{A}}^{+}(s)$ and $\hat{\bm{D}}^{-}(s)$ are vectors with components $\hat{A}_{n}^{+}(s)$ and $\hat{D}_{n}^{-}(s)$, respectively and

[TABLE]

where, for reference,

[TABLE]

and $\hat{\psi}_{n-m}(s)$ is given by (71). In the above $\theta_{S}$ and $S$ are chosen to satisfy

[TABLE]

Later we identify $S$ and $\theta_{S}$ as the effective wavenumber and transmission angle. The above does not determine the sign of $S$ for any given complex $s$. To fully determine $S$ and $\theta_{S}$, we take $\operatorname{sgn}(\mathrm{Re}\,s)=\operatorname{sgn}(\mathrm{Re}\,S)$ which together with (33) leads to

[TABLE]

where both $S$ and $\theta_{S}$, when considered as functions777With our choice of branch-cut, the simplest way to compute $S$, with most software packages, is to use $S=\operatorname{sgn}(\mathrm{Re}\,s)\sqrt{s^{2}+(k\sin\theta_{\mathrm{inc}})^{2}}$ with the default cut location for the square root function, i.e. along the real negative line. of $s$, contain branch-points at $s=\pm\textrm{i}k\sin\theta_{\mathrm{inc}}$ with finite branch-cut running between $-\textrm{i}k\sin\theta_{\mathrm{inc}}$ and $\textrm{i}k\sin\theta_{\mathrm{inc}}$. However, $\bm{\Psi}(s)$ is an entire matrix function having only zeros in $s$ and no branch-points; see the end of appendix A for details.

Determining the roots of $\det\bm{\Psi}({s})=0$ will be a key step in solving (27), and so the following identities will be useful

[TABLE]

where (35) results from (32) and (75). Equation (36) then follows from using (35)1, (29), and appendix C.

3.1 Multiple waves solution

To solve (27), we use a matrix product factorisation [42] of the form:

[TABLE]

where $\bm{\Psi}^{-}({s})$, and its inverse, are analytic in Im $s<\epsilon$, and $\bm{\Psi}^{+}({s})$, and its inverse, are analytic for Im $s>-\epsilon$. See (26) for the definition of $\epsilon$.

For our purposes, it is enough to know that such a factorisation exists [42], as this will lead to a proof that $\bm{A}(x)$ is a sum of attenuating plane waves.

Multiplying both sides of (27) by $[\bm{\Psi}^{-}({s})]^{-1}$ and by $({s}-{\alpha})_{-}$ leads to

[TABLE]

where $(s+\alpha)_{+}$ is analytic for $\mathrm{Im}\,s>-\mathrm{Im}\,\alpha$, while $(s-\alpha)_{-}$ is analytic for $\mathrm{Im}\,s<\mathrm{Im}\,\alpha$. We need to rewrite the last term above as a sum of a function which is analytic in the upper half-plane (Im $s>-\epsilon$) and another analytic in the lower half-plane. This is achieved below:

[TABLE]

where we define

[TABLE]

so that $\bm{g}^{-}(s)$ does not have a pole at $s=-\alpha$ and is therefore analytic for $\mathrm{Im}\,s<\epsilon$.

Substituting (39) into (38) leads to

[TABLE]

Because both sides are analytic in the strip $|\mathrm{Im}\,s|<{\epsilon}$, we can equate each side to $\bm{E}(s)$, some analytic function in the strip. Further, as the left-hand side (right-hand side) of (40) is analytic for $\mathrm{Im}\,s>\epsilon$ ( $\mathrm{Im}\,s<-\epsilon$), we can analytically continue $\bm{E}(s)$ for all $s$, i.e. $\bm{E}(s)$ is entire.

To determine $\bm{E}(s)$ we need to estimate its behaviour as $|s|\to\infty$. From (22) we have that $A^{+}(s)=\mathcal{(}|s|^{-1})$ as $|s|\to\infty$ in the upper half-plane, and from (29 - 31):

[TABLE]

for $s$ in the strip (26). From this we know that the factors $\bm{\Psi}^{+}(s)$ and $\bm{\Psi}^{-}(s)$ must be $\mathcal{O}(|s|)$ as $|s|\to\infty$, in their respective half-planes of analyticity [28]. So, the left hand-side of (40) behaves as $\mathcal{O}(|s|^{-1})$ as $|s|\to\infty$ in $\mathrm{Im}\,s>-\epsilon$. We can therefore use Liouville’s theorem to conclude that $\bm{E}(s)\equiv 0$, which means the Wiener-Hopf equation (40) is formally equivalent to

[TABLE]

Let $\mathbf{C}^{+}(s)$ be the cofactor matrix of $\bm{\Psi}^{+}(s)$, so that

[TABLE]

From the property (36)1 we can write $\det\bm{\Psi}(s)=f(s^{2})$ for some function $f$. Then, for every root $s=s_{p}$ of $\det\bm{\Psi}(s)$, with Im $s_{p}>0$, we have that $-s_{p}$ is also a root, and vice-versa. From here onwards we assume:

[TABLE]

For any truncated matrix $\bm{\Psi}(s)$, i.e. evaluating $m,n=-M,\ldots,M$ in (29), the roots $s_{p}$ are discrete. In section 5 we demonstrate asymptotically that they are indeed discrete for the limits of low and high wavenumber $k$. For the numerical results presented in this paper, we numerically solve the above dispersion relation for the truncating the matrix $\bm{\Psi}(s)$, and then increase $M$ until the roots converge (typically no more than $M=4$ was required).

Given $\det\bm{\Psi}(s)=\det\bm{\Psi}^{-}(s)\det\bm{\Psi}^{+}(s)$, every root of $\det\bm{\Psi}(s)$ must either be a root of $\det\bm{\Psi}^{-}(s)$ or a root of $\det\bm{\Psi}^{+}(s)$. For $[\bm{\Psi}^{+}(s)]^{-1}$ to be analytic in the upper half-plane, $\det\bm{\Psi}^{+}(s)$ must only have roots $s=-s_{p}$. As a consequence, $\det\bm{\Psi}^{-}(s)$ only has roots $s=s_{p}$.

To use the residue theorem below, we need to calculate $\det\bm{\Psi}^{+}(s)$ for $s$ close to the root $-s_{p}$, in the form

[TABLE]

where we use $\frac{d\det\bm{\Psi}}{ds}(-s_{p})$ instead of $\frac{d\det\bm{\Psi}^{+}}{ds}(-s_{p})\det\bm{\Psi}^{-}(-s_{p})$, because it is more difficult to numerically evaluate $\frac{d\det\bm{\Psi}^{+}}{ds}(-s_{p})$.

Using the above, and that $\mathbf{C}^{+}(S)$ is analytic for Im $s>-{\epsilon}$, we can apply an inverse Fourier transform (17)2 to both sides of (42) and using residue calculus we find

[TABLE]

For $x>0$, the integral over $s\in[-\infty,\infty]$ in (46) is, by Jordan’s lemma, the same as a clockwise integral over the closed contour $\mathcal{C}_{\bm{A}}$ which surrounds the poles $-s_{1},\,-s_{2},\ldots$, i.e. roots of (44), as shown by Figure 2. Note that the cofactor matrix $\mathbf{C}^{+}(s)$ contains no poles and so does not contribute additional residual terms. The yellow striped region in Figure 2 is the domain where $\bm{\Psi}$ is analytic. On the other hand, for $x<0$, the integral (46) is the same as an integral over the counter-clockwise closed contour within the region Im $s>0$ (not shown in Figure 2). The integrand has no poles in this domain and hence evaluates to zero.

Likewise, by applying an inverse Fourier transform to  (43), we obtain:

[TABLE]

For $x<0$ the above integral is the same as a counter-clockwise closed integral over $\mathcal{C}_{\bm{D}}$ which surrounds the pole $s=\alpha$ (recalling that Im $\alpha>0$), as shown in Figure 2. The result is just the residue at this pole. That is, the function $\bm{\Psi}^{-}(s)\bm{g}^{-}(s)$ contains no other singularities within $\mathrm{Im}\,s>0$. On the other hand, for $x>0$ the integral is the same as a closed clockwise integral around the region Im $s<0$ which evaluates to zero, as there are no singularities in this region (not shown in Figure 2).

Clearly (46) shows that $\bm{A}(x)$ is a sum of plane waves with different effective wavenumbers $s_{p}$, each satisfying (44). In section 5 we discuss these roots in more detail, and in section 6, we see that usually only a few effective wavenumbers are required to obtain accurate results.

3.2 Reflection coefficient

By substituting (48) in (16) leads to

[TABLE]

Alternatively, the reflection coefficient can be calculated from (12) by employing the form of $\bm{A}(x)$ from (46), which is the more common approach. To simplify, we use

[TABLE]

which then implies that $\psi_{n}(x_{1}-x)=\frac{2}{\alpha}\textrm{i}^{n}\textrm{e}^{-\textrm{i}n\theta_{\mathrm{inc}}}\textrm{e}^{\textrm{i}\alpha(x_{1}-x)}$ for $x_{1}\geq x$. The above is shown in [43, equation (37)] and [11, equation (65)]. This result together with (46) substituted into (12) leads to the form

[TABLE]

where we used that Im $s_{p}>0$. The above agrees with [13, equation (39)] and888When taking a zero thickness boundary layer, i.e. $J=0$, and appropriate substitutions. [18, equation (6.9)].

4 Monopole scatterers

For particles that scatter only in their monopole mode, i.e. the scattered waves are angularly symmetric about each particle, we can easily calculate the factorisation (37). This type of scattered wave tends to dominate in the long wavelength limit for scatterers with Dirichlet boundary conditions. In acoustics, these correspond to particles with low density or low sound speed.

Once we know the factorisation (37), we can then calculate the average scattering coefficient (46) and average reflection coefficient (49). We will compare both of these against predictions from other methods in section 6.

4.1 Wiener-Hopf factorisation

For scalar problems, there are well known techniques to factorise $\Psi_{00}(s)=\Psi_{00}^{-}(s)\Psi_{00}^{+}(s)$, such as Cauchy’s integral formulation, for details see [19, Section 5. Wiener-Hopf Technique] and [21].

For monopole scatterers we use $S^{2}-k^{2}=s^{2}-\alpha^{2}$ and rewrite

[TABLE]

with $\mathrm{N}_{0}(bS)$ given by (31). Then, because $q(s)\to 1$ as $|s|\to\infty$, we can factorise $q(s)=q^{-}(s)q^{+}(s)$ using

[TABLE]

where the integral path for $q^{+}(s)$ ($q^{-}(s)$) has to be in the strip where $q(s)$ is analytic, with the path for $q^{+}(s)$ ($q^{-}(s)$) passing below (above) $z$. We then have999Note that the factors $q^{+}(s)$ and $q^{-}(s)$ are singularity and pole free in their respective regions of analyticity, and so their inverses $[q^{+}(s)]^{-1}$ and $[q^{-}(s)]^{-1}$ have the same property. that

[TABLE]

where (54)3 holds if $-s$ is below the integration path of (53) and $s$ is above the integration path of (52). From (46) we see that we need only evaluate $\Psi^{+}_{00}(s)$, and therefore $q^{+}(s)$, for $s=s_{1},s_{2},\ldots,s_{p}$ where as $p$ increases, the $s_{p}$ become more distant from the real line. Then for large $z$, by inspection of (31), we have that

[TABLE]

and therefore we can accurately approximate the integral (52) by truncating the integration domain for large $z$.

4.2 Explicit solution for monopole scatterers.

For monopole scatterers $A_{n}(x)=D_{n}(x)=0$ for $|n|>0$. Using this in (28 - 31) leads to all vectors and matrices having only one component, given by setting $n=m=0$. In this case $\bm{A}$ (46) reduces to

[TABLE]

for $x>0$, where we used (54), $\mathbf{C}^{+}(s)=1$, $\bm{B}=\textrm{i}T_{0}$, and $\frac{d\Psi_{00}}{ds}(-s)=-\frac{d\Psi_{00}}{ds}(s)$ for every $s$. Likewise for (49) we arrive at

[TABLE]

Alternatively, using (51), we can calculate the contribution of $P$ effective waves to the reflection coefficient

[TABLE]

where the error $|\mathfrak{R}^{P}-\mathfrak{R}|$ then indicates how many effective waves are needed to accurately describe the field near the boundary $x=0$.

5 Multiple effective wavenumbers

Equation (46) clearly shows that $\bm{A}(x)$ is a sum of attenuating plane waves, each with a different effective wavenumber $s_{p}$. These $s_{p}$ satisfy the dispersion equation (44):

[TABLE]

with $\bm{\Psi}$ given by (30) and the first identity follows from (36).

An important conclusion from $\det\mathbf{G}(S_{p})=0$ is that the wavenumbers $S_{p}$ are independent of the angle of incidence $\theta_{\mathrm{inc}}$. We focus on showing the results for $S_{p}$, rather than $s_{p}$, because then we do not need to specify $\theta_{\mathrm{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

[TABLE]

where $a_{o}$ is the radius of the particle.

With the above parameters, we found that truncating the matrix $\bm{\Psi}(s)$, with $|n|\leq 3$ and $|m|\leq 3$ in (29-31), led to accurate results when calculating the effective wavenumbers $S_{p}$, i.e. the roots of (44). Numerically calculating the wavenumbers $S_{p}$ then leads to Figure 3.

The effective wavenumbers with the lowest attenuation (smallest imaginary part) contribute the most to the transmitted wave. In Figure 3 we see two wavenumbers have lower attenuation then the rest, both within the dashed grey circle. The blue point represents the wavenumber that most of the literature focuses on calculating: it has a positive real part and therefore propagates forwards along the $x-$axis (into the material) as is expected for a transmitted wave. However, the other wavenumber, with negative real part, is equally as important because it actually has lower attenuation. Figure 1 illustrates several effective wavenumbers, some travelling forward into the material, while others have negative phase direction (travel backwards).

In Figure 3 we see what appears to be an infinite sequence of effective wavenumbers $S_{p}$, where $|S_{p}|\to\infty$ as $p\to\infty$. To confirm their existence, and to find their locations as $|p|\to\infty$, we develop asymptotic formulas in appendix B. The results of the asymptotics are summarised below.

For monopole scatterers, where $n=m=0$ in (29), equations (82) give the effective wavenumbers $S_{p}^{o}$ at leading order:

[TABLE]

and for any integer $p$. We use the superscript “o” to distinguish these wavenumbers for monopole scatterers from others. Even though (60) was deduced for large integer $p$, it gives remarkably agreement with numerically calculated wavenumbers, except for the two lowest attenuating wavenumbers, as shown in Figure 4. In the figure we denoted $S_{*}^{o\pm}$ as the effective wavenumber that can be calculated by low volume fraction expansions [11, 14].

For multipole scatterers, where both $n$ and $m$ could potentially range from $-\infty$ to $\infty$ in (29), we can also calculate an infinite sequence of effective wavenumbers. To show this explicitly, we consider the limit of large $bk$, with $|k|\sim|S|$. In the opposite limit $bk\ll 1$, the Rayleigh limit, only one effective wavenumber is required [44, 45].

At leading order, the asymptotic solution of (86) leads to the effective wavenumbers:

[TABLE]

for integer $p$. This confirms that there are an infinite number of effective wavenumbers for large scatterers, i.e. $bk\gg 1$. The distribution of these wavenumbers is similar to the monopole wavenumbers shown in Figure 4.

These asymptotic formulas (60) and (62) demonstrate the existence of multiple effective waves in the limit of small (monopole and Dirichlet) scatterers (60) and large scatterers (62). However, neither of these formulas, nor the low volume fraction expansions of the wavenumber [11], are able to accurately estimate the low attenuating backward travelling effective wavenumber such as $S_{1}$ shown in Figure 3 (in this case not related to the $S_{1}^{o\pm}$ and $S_{1}^{k\pm}$ given above). There is currently no way to analytically estimate these types of wavenumbers, even though they are necessary to accurately calculate transmission due to their small attenuation. The only approach it seems is to numerically solve (44).

6 Numerical results

Here we present numerical results for monopole scatterers, as these have explicit expressions for reflection (56) and the transmitted wave (55) (or more accurately the average scattering coefficients). We compare our analytic solution with a classical method that assumes only one effective wavenumber [11, 13], and the Matching Method [18], recently proposed by the authors. It should be noted that all of these approaches aim to solve the same equation (14).

Note that for monopole scatterers, using only one effective wavenumber $s_{1}$ can, in some cases, lead to accurate results. However, for multipole scatterers (a more common scenario practically) this is rarely the case because, as shown by Figure 3, there can be at least two effective wavenumbers with low attenuation, and therefore both are needed to obtain accurate results.

For the numerical examples we use the parameters

[TABLE]

which implies that the number fraction $\mathfrak{n}\approx 0.38$ per unit area. When we choose to fix the wavenumber, as we do for Figure 5 and Figure 6, we use $bk=1.001$. This leads to a wavelength ($2\pi/k$) which is roughly six times larger than the particle diameter. If the particle was, say, more than a hundred times smaller than the wavelength, then only one effective wavenumber in the sum (55) would be necessary to accurately calculate $A_{0}(X)$.

To start we compare the average scattering coefficient $A_{0}(x)$ calculated by the Wiener-Hopf solution (55) with other methods in Figure 5. The most accurate of these other methods is the Matching Method [18, 46], and it closely agrees with the Wiener-Hopf solution when using $352$ effective wavenumbers. The exception is the region close to the boundary $x=0$, where the Wiener-Hopf solution experiences a rapid transition. The low volume fraction method is the most commonly used in the literature: it assumes a small particle volume fraction101010For the low volume fraction method we used a small volume fraction expansion for the wavenumber, but we numerically evaluated the wave amplitude. This is because the alternative, a small volume fraction expansion of the wave amplitude, led to poor results. and just one effective wavenumber [11, 13]. One significant conclusion we can draw from Figure 5 is that both the low volume fraction method and $A_{0}^{1}\textrm{e}^{\textrm{i}s_{1}x}$ are inaccurate near the boundary $x=0$. This means that both of these methods lead to inaccurate reflection coefficients.

In general, the Wiener-Hopf method does not lead to an explicit formula for the reflection coefficient (49), because we do not have an exact factorisation (37) for any truncated square matrices. However, there are methods [18, 13, 20, 29, 30, 34] to calculate $A_{n}(x)$, from which we can obtain the reflection coefficient (12). The method [18] also accounts for multiple effective wavenumbers. So one important question is: when using (12), how many effective wavenumbers do we need to obtain an accurate reflection coefficient?

In Figure 6 we show how increasing the number of effective waves $P$ reduces the error between $\mathfrak{R}^{P}$ (57) and $\mathfrak{R}$ (56). To calculate a highly accurate reflection coefficient $\mathfrak{R}$, we could use either (56) or the Matching Method [18, 46], as both give approximately the same $\mathfrak{R}$.

Now we ask: how does the reflection coefficient (57), deduced via the Wiener-Hopf technique, compare with other methods across a broader range of wavenumbers. The result is shown in Figure 7, where $\mathfrak{R}^{O}$ is a low volume fraction expansion111111We 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. of just one effective wavenumber [13]. The reflection coefficient $\mathfrak{R}^{M}$ is calculated from the Matching Method [18, 46]. The general trend is clear: $\mathfrak{R}^{O}$ becomes more inaccurate as we increase the background wavenumber $kb$. On the other hand both $\mathfrak{R}^{M}$ and $\mathfrak{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_{0}(x)$ as a superposition of truncated waves, which is only accurate as long as the discarded terms are small. So, for a truncation number $P$, we can expect the instability to occur when $\textrm{e}^{\textrm{i}\mathbf{s}_{P}x}$ is not small, i.e. $x\approx 1/\mathrm{Im}\,s_{p}$. However, this instability does not affect the accuracy of the reflection coefficient (56) deduced by the Wiener-Hopf technique, as demonstrated by close agreement with the Matching Method in Figure 7.

7 Conclusion and Next Steps

The major result of this paper is to prove that the ensemble-averaged field in random particulate materials consists of a superposition of waves, with complex effective wavenumbers, for one fixed incident wavenumber. These effective wavenumbers are governed by the dispersion equation (58) and are independent of the angle of incidence $\theta_{\mathrm{inc}}$. We showed asymptotically in  section 5 that this has an infinite number of solutions, and hence there are an infinite number of effective wavenumbers. The Wiener-Hopf technique also provides a simple and elegant expression for the reflection coefficient (49), whose form can be used to guide and assess methods to characterise microstructure [47, 48].

To numerically implement the Wiener-Hopf technique, we considered particles that scatter only in their monopole mode in section 6. There we saw that when close to the interface of the half-space, a large number of effective wavenumbers were necessary to reach accurate agreement with an alternative method from the literature, the Matching Method as introduced by the authors in [18]. To obtain a constructive method via the Wiener-Hopf technique for general scatterers, and not just monopole scatterers, will require the factorisation of a matrix-function [31], which is challenging. For these reasons the Matching Method [18] is presently more effective than using the Wiener-Hopf technique. However, there is ongoing work to use approximate methods [49, 33, 34] which exploit the symmetry and properties of the matrix (29).

Moving forwards, this paper together with [18], establish accurate and robust solutions to the governing equation (14). These same methods can now be translated to three spatial dimensions and vectorial waves (e.g. elasticity and electromagnetics), with much of the groundwork already available [12, 40, 16]. Some clear challenges, that can now be addressed, are to verify the accuracy of the statistical assumptions used to deduce (14). These include the hole-correction and the quasicrystalline approximations. As these are now the only assumptions used, we could compare the solution of (14) with multipole methods [50, 51] in order to investigate their accuracy and limits of validity.

Appendix A The Fourier transformed kernel $\hat{\psi}_{n}(s)$

Here we calculate the Fourier transform (17) of $\psi_{n}(X)$ (13). To do so, it is simpler to use

[TABLE]

Note that both $\mathrm{F}_{n}(X,Y)$ and $\textrm{e}^{\textrm{i}({s}X+Yk\sin\theta_{\mathrm{inc}})}$ satisfy wave equations, with

[TABLE]

where we used (13)2 for the first equation and (33) for the second equation. This means that we can use Green’s second identity to obtain

[TABLE]

for any area $\mathcal{B}$ in which the integrand is analytic, where $\bm{n}$ is the outwards pointing unit normal and $\mathrm{d}z$ is a differential length along the boundary $\partial\mathcal{B}$. To calculate $\hat{\psi}_{n}({s})$, we take the region $\mathcal{B}$ to be defined by $R\geq b$, with $(R,\Theta)$ being the polar coordinates of $(X,Y)$, in which case the integral over $\mathcal{B}$ converges because as $R\to\infty$ we have that

[TABLE]

exponentially fast when $|\mathrm{Im}\,({s})|<\mathrm{Im}\,k(1-|\sin\theta_{\mathrm{inc}}|)$. Under this restriction, and by assuming $S\not=\pm k$, (67) then leads to

[TABLE]

by using ${s}X+Yk\sin\theta_{\mathrm{inc}}=RS\cos(\theta-\theta_{S})$ from (33) and

[TABLE]

where $\mathrm{J}_{n}$ is the Bessel function of the first kind, and we used the Jacobi-Anger expansion on $\textrm{e}^{\textrm{i}SR\cos(\Theta-\theta_{S})}$, integrated over $\Theta$ and used the identity $\mathrm{J}_{n-1}(SR)-\mathrm{J}_{n+1}(SR)=2\mathrm{J}_{n}^{\prime}(SR)$. In summary

[TABLE]

when the condition (24) is satisfied, with $\mathrm{N}_{n}$ given by (31).

Below we establish some useful properties for $\hat{\psi}_{n}(s)$. In particular, we show that $\hat{\psi}_{n}(s)$ has no branch-points.

The function $\mathrm{N}_{n}(bS)$, for integer values of $n$, can be expanded around $S=0$ as

[TABLE]

where the $c_{m|n|}$ are some constants that depend on $m$ and $|n|$, and the radius of convergence of the series above is infinite. Using (33) we can write

[TABLE]

Substituting (72) and (73) in (71) results in

[TABLE]

Because $S^{2}=s^{2}+k^{2}\sin^{2}\theta_{\mathrm{inc}}$, we can immediately see from the above that $\hat{\psi}_{n}(s)$ has no branch-points. Additionally we can establish the properties:

[TABLE]

Appendix B Asymptotic location of the wavenumbers

Here we explicitly calculate a sequence of effective wavenumbers $S_{p}$, assuming $|S_{p}|$ large and increasing with $p$, and $\mathrm{Im}\,S_{p}>0$, that asymptotically satisfy (58). A key step is to approximate the terms appearing in (30), such as

[TABLE]

for large $|bS|$, where the terms $\textrm{e}^{\textrm{i}bS}$ are discarded as Im $bS\to\infty$.

Monopole scatterers

The simplest case is for monopole scatterers, where $n=m=0$ in (30), and the effective wavenumber $S$ satisfies

[TABLE]

where $c=\sqrt{2\pi}\mathfrak{n}b^{2}T_{0}\mathrm{H}^{(1)}_{0}(kb)\textrm{e}^{-\frac{\textrm{i}\pi}{4}}$. Here we used (76), and ignored terms which are algebraically smaller than $bS$. To find the root of the above we substitute

[TABLE]

where $x$ and $y$ are real, and $|x|$ and $y$ are large with $y>1$. This leads to

[TABLE]

For the logarithm and square root we use the typical branch cut $(-\infty,0)$ and take positive values of the functions for positive arguments. For the above to be satisfied to leading order then $x^{3/2}\sim y$, which reduces the above equation to

[TABLE]

where we substituted $c=r_{c}\textrm{e}^{\textrm{i}\theta_{c}}$, for real scalars $r_{c}$ and $\theta_{c}$. Equating the real and imaginary parts of the above leads to

[TABLE]

for integers $p$. From this we can identify that, at leading order, the effective wavenumbers are given by (60).

Multipole scatterers:

With the same method used above, we can also demonstrate the existence of multiple effective wavenumbers for $n,m=-M,-M+1,\cdots,M$ in (30). To show this explicitly, we consider $bk$ to be the same order as $bS$, that is $|k|\sim|S|$.

By considering $bk$ large, we can approximate

[TABLE]

combining this with (76) and considering $|k|\sim|S|$, the term (30) at leading order becomes

[TABLE]

where

[TABLE]

By simple rearrangement of the determinant we find that121212The determinant of $b^{2}\bm{G}$ equals the product of its eigenvalues. The eigenvector $(T_{-M},\cdots,T_{M})^{\mathrm{T}}$ gives the eigenvalue $d_{0}+c_{0}\sum_{m}T_{m}$, while all other eigenvalues equal $d_{0}$.

[TABLE]

Note that $d_{0}\not=0$, i.e. $S\not=\pm k$, was necessary to reach the condition (24), which was used to calculate the Fourier transforms (27). Taking this into consideration, and taking the limit $M\to\infty$, the effective wavenumbers $S$ must satisfy

[TABLE]

Using an asymptotic expansion analogous to (78), the above leads to the effective wavenumbers (62).

Appendix C Equivalent determinants

For any square matrices $\mathbf{A}$ and $\bm{B}$, and scalar $c$, if $A_{nm}=B_{nm}c^{n-m}$ (not employing the summation convention), then

[TABLE]

This follows simply by defining the diagonal matrix $C_{nm}=\delta_{nm}c^{n}$, which leads to $\bm{A}=\bm{C}\bm{B}\bm{C}^{-1}$, and $\det(\bm{C}\bm{B}\bm{C}^{-1})=\det\bm{C}\det\bm{B}\det\bm{C}^{-1}=\det\bm{B}$.

\dataccess

We provide the code to generate all the graphs in [46].

\aucontribute

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

We have no competing interests.

\funding

This work was funded by EPSRC (EP/M026205/1,EP/L018039/1,EP/S019804/1) including via the Isaac Newton Institute (EP/K032208/1,EP/R014604/1) and partial support from the UKAN grant EPSRC (EP/R005001/1).