Extended Foldy–Lax Approximation on Multiple Scattering

The Foldy–Lax self-consistent system has been widely used as an efficient numerical approximation of multiple scattering of time harmonic wave through a medium with many scatterers when the relative radius of each scatterer is small and the distribution of scatterers is sparse. In this paper, an “extended” Foldy–Lax selfconsistent system including both source and dipole effects as well as corrections due to the self-interacting effects will be introduced, in which the scattering amplitudes and the corrections are determined as powers of the small scaled radius. This new approach substantially improves the accuracy of the approximation of the original Foldy–Lax approach.


Introduction
The multiple scattering by finite-size inhomogeneities has attracted interests of many researchers and the finite-size effect has been treated in different aspects. For instance, the multiple scattering has been used in time-reversal imaging [4] on reflecting more signals towards the aperture thus virtually increasing the aperture size, which is called super-resolution [2], and beyond. In the real computation, since the number of scatterers is large, proper approximating method is necessary for multiple scattering modelling.
In this paper, we consider the multiple scattering by many small scatterers in the low frequency regime. Here low frequency means low relative frequency, i.e., the situation when the wave length is large relatively to the radii of the scatterers. Without lose of generality, we assume all the scatterers are identical spheres for simplicity. This is reasonable when we concern small size point-like scatterers, instead of extended scatterer or target, otherwise we may adapt the generalized Foldy-Lax formulation devised in [9] via the coupling of a boundary integral equation and the Foldy-Lax system. Take k as the wave number of the incident wave and R as the radius of each scatterer, then low relative frequency means kR small.
We restrict our analysis on the time harmonic case. The multiple scattering problem can be written by standard inhomogeneous Helmholtz equation and the solution is represented by the well-known Lippmann-Schwinger integral equation [3], which can be solved exactly by separation of variables and using spherical harmonics, or, high-accuracy numerical techniques have also been developed to address this problem, for example [1]. All of these methods can be theoretically used to solve the multiple scattering problem accurately. These techniques have large computational time and memory overheads, and thus, in practice, have been used to model only moderate sections of the complete problem. In reality, it is necessary to model large-scale distributed scatterers by appropriate semi-analytical approximations, such as the widely accepted Foldy-Lax approximation [5,6,10,12,13,14,17]. In this approach, the free space propagator, which is the Green's function of Helmholtz equation, is inevitably used to simulate the wave propagation. Specifically, write the total field as where u i is the given incident field and u j sc is the field scattered by the j-th scatterer, M is the number of scatterers. Define the effective field which is the 'radiation incident on the n-th scatterer' in the presence of all the other scatterers. Now the problem is written as where T j is an operator relating the field incident on the j-th scatterer, u j , to the field scattered by the j-th scatterer, u j sc , hence This is the so-called Foldy-Lax self-consistent system, or simply Foldy's system. If one could solve u n , the scattered field u n sc and the total field u would be easily arrived.
The Foldy-Lax method considers only the first order approximation, in which the dipole-effects and self-interaction effects are neglected. Actually, the dipole-effects (without error estimates) has been considered by Lax [10] and also the the book of Martin [12]. However, in this paper, we will show more precisely the dipole-effect, and further, we will show how to include the selfinteraction effect in an extended Foldy-Lax approximation, thus substantially improve the approximation accuracy. In this approach, the total field in free space will be written by here α j and β j are the scattering amplitudes of the j-th scatterer,Ḡ and vector g are free space 'propagator' (which will be specified later), 'effective' wave u j and vector v j will be solved by the extended Foldy-Lax system in the form Here the first term in the summation is the source while the second term enclose the dipole effect. Furthermore, η l and θ l are correction terms due to the self-interacting effect. The extended Foldy-Lax system (1.2) is an extension of Foldy-Lax system (1.1), since (1.1) can be derived by taken η l 's and β j 's be zeros in the first equation of (1.2). With both the self-interaction effect and dipole effect being included in this new approach, the accuracy of the approximation is substantially improved. We also note that, this new formulation can be used to real application as mentioned in [7,8,9]. The paper is organized as following. In the next section, we will specify the multiple scattering problem and write the exact solution. Then we will consider the approximation of the exact solution in Section 3, and derive the extended Foldy-Lax system in Section 4. Finally the accuracy of the extended Foldy-Lax system will be given in Section 5.

Multiple Scattering Problem
The propagation of time harmonic waves in inhomogeneous medium is governed by Helmholtz equation where k is the wave number, n(x) = c 0 /c(x) is the index of refraction, which is the ratio of the wave speed in the homogeneous background to that in the inhomogeneous medium with scatterers. It will be convenient to define the scattering potential q(x) := n 2 (x) − 1, then the Helmholtz equation can be rewritten as Considering that the inhomogeneities consist of M disjoint identical spherical scatterers located in a bounded domain Ω ⊂ R 3 , each scatterer centered at ξ j with radius R, |ξ l − ξ j | > γR for any l = j, l, j = 1, . . . , M , γ > 2 is a constant due to the problem considered, which characterize the sparsity of the scatterers. Thus the scattering potential can be taken as piecewise constant where τ j is a given constant for each j and χ R is the characteristic function with radius R.
It is well known that the outgoing Green's function of Helmholtz equation is given by 4) and the solution can be written as [3,12] where the incident wave u i is an entire solution of the Helmholtz equation Using the scattering potential q in the summation form (2.3), we can accordingly write the solution (2.5) as where B = {y | |y| < R}. This equation is the well-known Lippmann-Schwinger integral equation [3] and holds for all x ∈ R 3 . One can firstly solve the Fredholm type integral equation inside each B i := {y | |y − ξ i | < R} and then evaluate the field in the free space when the field inside each scatterer is known. The equation (2.6) is defined in the whole space and the solution u is analytic in each scatterer and the free space, continuous across the boundary of each scatterer, see [15,16]. The accurate solution of this approach for large problem is formidable thus we turn to an appropriate semi-analytical approximation called Foldy-Lax approximation. However, the Foldy-Lax method considered only the first order approximation, in which the dipole-effects and self-interaction effects are neglected. The newly introduced extended Foldy-Lax method will include both the dipole-effects and self-interaction effects.
For further use, we take the gradient on both sides of (2.6) to get the gradient of u away from the boundary of each scatterer as Notice that the Green's function has only weakly singularity, it is easy to see that the integrations above are well-defined. Note also that we always put a scaling factor 1 k in front of the gradient operator here and hereafter. Particularly, we evaluate u and 1 k ∇u at each ξ l thus where we have used the notation ξ jl = ξ l − ξ j . Starting from the above integral system, we will derive an extended Foldy-Lax system when the scaled radius of each scatterer kR is small and determines the scattering amplitudes in the extended Foldy-Lax self-consistent system.

Approximation of the Solution u(x) in Free Space
For the wave field outside scatterers, by the Lippmann-Schwinger equation (2.6), we have Since u is regular in each scatterer, we will derive the asymptotic expansion of each integral in terms of the power series of the scaled radius kR. For x ∈ R 3 \ M j=1 B(ξ i , R) and |ky| ≤ kR, by Taylor expansion, For the spherical scatterers, due to the symmetrical property, it is easy to check where δ ij is the Kronecker delta function, thus substitute the above into (3.1) then we get where thus we have the formula for One can use this formula to evaluate the wave field in the free space when u(ξ j ) and 1 k ∇u(ξ j ) are known for each j, which will be given by solving extended Foldy-Lax system derived in the next section.  in which the last term distinguishes the self-interacting term.

Evaluation of the integral term for each j = l.
By Taylor expansion, for |ky| ≤ kR,

Similar calculation as in last section yields
where α j , β j are defined by (3.3).

Evaluation of the integral term with ξ l .
When |ky| ≤ kR, we use After direct calculation we find Now, substitute (4.2) and (4.3) into (4.1) we have: where α j , β j and η l are defined by (3.3) and (4.4) respectively. The value of u(ξ l ) evolves 1 k ∇u(ξ j ), we still need to evaluate the gradients of the field at each scatterer.

Evaluation of the gradient of the solution
Rewrite the gradient of the solution at ξ l in the second equation of (2.7) as we then derive the asymptotic expansion of each term above. For the integral term of each j = l, similar calculation like the one in last subsection yields where the coefficients α j and β j are defined by (3.3). The integral term with ξ l in (4.6) is a little bit different from the corresponding term in (4.1), since the integration here involves the gradient of the Green's function but the original integration deals only G. The gradient of Green's function is by direct calculation we find Next, substitute (4.7) and (4.8) into (4.6): where α j , β j and θ l are defined by (3.3) and (4.9) respectively.

The extended Foldy-Lax system
Combining (4.5) and (4.10) together we have It is obvious that the number of scatterers M is no more than the order of ( 1 kR ) 3 , then the higher order term in above system is of O((kR) 6 ). Omitting the higher order term, we have then we have this is a self-consistent system, which is called extended Foldy-Lax system. Here α j and β j , defined in (3.3), are the scattering amplitudes of the j-th scatterer, the first term in the summation is the source while the second term enclose the dipole effect. Furthermore, η l and θ l , defined in (4.4) and (4.9), are correction terms due to the self-interacting effect. Recall the field outside the scatterers is given by (3.4), which can be approximated as , one can evaluate u in the free space once the extended Foldy-Lax system (4.11) has been solved. Note that the linear system (4.11) can be rewritten in matrix form, with coefficient matrix being diagonal dominant when kR is sufficiently small thus uniquely solvable.
Foldy [5] and Lax [10] assume that β j = 0, η l = 0 and θ l = 0 in (4.11)-(4.12). The main feature of this approach is that the extended system includes the selfinteraction effect and dipole effect. We further note that, although the above computation is derived for spherical scatterers, theoretically it can be directly generalized to more general shape of scatterers, by adapting the computation in Section 3 with B replaced accordingly. However the computation will depend on the configuration of different shapes, and may only be very much practical for simple geometry.

Accuracy Check
In this section, we will check the accuracy of the the approximation given by extended Foldy-Lax self-consistent system. To achieve this purpose, we compare the difference from exact wave field at each scatter given by Lippmann-Schwinger integral equation to the approximate field derived from the extended Foldy-Lax self-consistent system [11]. Our main result in this section is the following theorem. Theorem 1. Give a bounded domain Ω ⊂ R 3 enclosing M disjoint identical spherical scatterers, each one centered at ξ j with radius R, |ξ l − ξ j | > γR for l = j, l, j = 1, . . . , M and γ > 2. Assume kR be small and kR · M ∼ 1. Let u(ξ l ) and 1 k ∇u(ξ l ) be the exact wave field and its gradient at ξ l given by (2.7), u l and v l be the approximate field defined recursively by (4.11). Then we have Proof. The theorem can be proved by the method used in [11]. First, we introduce the intermediate quantities α jḠ (ξ jl )u(ξ j ) + β j g(ξ jl ) · 1 k ∇u(ξ j ) + η l · u(ξ l ), (5.2) regime is much better than the original Foldy-Lax approach. Recall that in the same setting of Theorem 1, the result in [11] is M l=1 u(ξ l ) − u l 2 ≤ O(kR) 1+ 4 5 /γ 2 for the original Foldy-Lax approach, where γ > 2 characterize sparsity of distribution of scatterers. In the present analysis, we don't ask the sparsity but included the dipole effects and the self-interaction corrections, which give an upper bound of the error as order O(kR) 11 . This suggest that the dipoleeffects and self-interaction effects are essential on improving the accuracy of multiple scattering approximation, despite they have been neglected in the original Foldy-Lax approach which has been widely used. This new formulation is expected to be used in [7,8,9] or other applications, to improve the numerical accuracy of small scatterer multiple scattering in wave simulation.