On the Convergence of the Multipole Expansion Method

The multipole expansion method (MEM) is a spatial discretization technique that is widely used in applications that feature scattering of waves from circular cylinders. Moreover, it also serves as a key component in several other numerical methods in which scattering computations involving arbitrarily shaped objects are accelerated by enclosing the objects in artificial cylinders. A fundamental question is that of how fast the approximation error of the MEM converges to zero as the truncation number goes to infinity. Despite the fact that the MEM was introduced in 1913, and has been in widespread usage as a numerical technique since as far back as 1955, to the best of the authors' knowledge, a precise characterization of the asymptotic rate of convergence of the MEM has not been obtained. In this work, we provide a resolution to this issue. While our focus in this paper is on the Dirichlet scattering problem, this is merely for convenience and our results actually establish convergence rates that hold for all MEM formulations irrespective of the specific boundary conditions or boundary integral equation solution representation chosen.

1. Introduction. The multiple scattering of waves is an important topic that arises in a variety of scientific fields including acoustics, electromagnetics, elasticity, water waves, and quantum mechanics. In the frequency domain, it is well known that the scattering of waves from multiple disjoint circular cylinders (spheres in three dimensions) can be computed with exceptional efficiency using a meshless technique called the multipole expansion method (MEM), a spatial discretization technique that involves truncating infinite series of multipoles [28,18,36,31,18].
The idea of applying multipole expansion techniques to multiple scattering problems can be traced back to 1913 with Zȧviška introducing it in [37] to compute the scattering of waves from an array of parallel cylinders. In 1955, Row used the MEM to obtain a numerical solution to a multiple scattering problem [34,31]. Since then, the MEM has appeared in countless fundamental and applied works which feature multiple scattering from cylinders or spheres.
The MEM also serves as a key component in several other numerical methods in which scattering computations involving arbitrarily shaped objects are accelerated by enclosing the objects in artificial cylinders or spheres. For instance, MEM-based formulations appear in scattering matrix methods [28,27], T-matrix methods [22,19], and Dirichlet-to-Neumann methods [23,3].
Recently, the MEM has seen extensive use in the field of metamaterials where it has been combined with lattice summation techniques to allow for efficient computational simulations in problems featuring infinitely periodic lattices such as photonic and phononic crystals, and metasurfaces; see, for instance, [10] and the monograph [9]. In the last few years, it has also been used to simulate scattering in the context of topological insulators [8,7], and subwavelength resonance models of the cochlea [6,5]. In addition, the MEM has been used in investigations of speckle statistics and non-invasive optical focusing in random scattering media [15,30]. It is worth noting that these latter works utilized the open-source MEM scattering library µ-diff which was released in 2015 [36]. Another interesting application of µ-diff can be found in [24] where it was used to generate the training data for a neural network aimed at localizing multiple scattering objects. We use µ-diff to validate our theoretical results in this paper; in fact our results establish the convergence theory for this library. The key parameter in any multipole-based scattering technique is the truncation number. The truncation number stipulates the number of terms that should be retained when one truncates the infinite series used to represent the problem in order to obtain a finite-dimensional discretization. Convergence theories for other multipole-based scattering methods such as the fast multipole method, the T-matrix method, and the method of fundamental solutions can be found in [21,22,16].
In the case of the MEM, multiple scattering of waves between cylinders results in a coupling of the coefficients of the infinite series associated with each cylinder, and this phenomenon has a pronounced effect on the decay of the approximation error of the MEM as the truncation number N → ∞. However, despite the fact that the MEM originated over a century ago, and has been in widespread use for well over fifty years, to the best of the authors' knowledge, the asymptotic rate of convergence of the approximation error of the MEM has not been properly quantified in the literature. Henceforth, we shall use the phrase 'convergence of the MEM' to refer to the rate of convergence of the approximation error of the MEM to zero.
Numerical investigations on the convergence of the MEM have recently been undertaken in [18,17]. There have been several other works wherein numerical investigations have been performed to ascertain the performance of iterative methods applied to the MEM system of equations [4,11,28]; it should be noted, however, that these iterative methods are converging to the approximate solution given by the MEM, hence, the underlying question of the asymptotic rate of convergence of the MEM solution to the true solution remains unanswered. In this paper we present a resolution to the long-standing problem of quantifying the decay of the MEM approximation error. The system of equations we consider first arose in Row's 1955 paper [34,Eq. (3) and Eq. (5)].
Let {Ω p } M p=1 be a set of disjoint circular cylinders in R 2 . Let the incident field be given by a point source located at x 0 . Denote by O p the center of cylinder Ω p , and by a p its radius. Denote by d pq = |O p − O q | the distance between the centers of cylinders Ω p and Ω q . Denote by d px0 = |O p − x 0 | the distance between the center of cylinder Ω p and the point source located at x 0 . Denote by Denote by E(N ) the approximation error of the MEM for the M cylinder system. The following theorem provides an asymptotic bound on the convergence of E(N ).
for a point source, for a plane wave. (1.1) The notation ' ', which will be made precise later, can be interpreted as signifying the left hand side is less than the right hand side up to an asymptotically irrelevant sub-exponentially increasing factor. In the point source case, the first term on the right hand side of the inequality represents the approximation error directly associated with the incident field, while the second term represents the approximation error associated with the geometry; thus, the bound indicates that if the point source is sufficiently close to one of the cylinders, then it is this placement of the point source that ultimately dictates the rate of convergence, otherwise the convergence is dictated by the pair of cylinders that maximizes the expression a p /(d pq − a q ). Numerical simulations verify that γ 1 (N ) provides a tight characterization of the convergence of the MEM as long as the cylinders are not too close together. When some cylinders are, in fact, in close proximity to one another, near-trapping of energy occurs as waves repeatedly reflect among theses cylinders, and thus it takes longer for energy to leak away to infinity. This behavior manifests itself in a decrease in the rate of convergence of the MEM. This decrease in the rate of convergence is captured by γ 1 (N ), but having said that, it transpires that the bound becomes overly pessimistic. The issue is that, while the a p /(d pq − a q ) term in γ 1 (N ) characterizes interactions among the cylinders, it doesn't fully account for the phenomenon of multiple scattering. However, the derivation of Theorem 1.1 relies on the asymptotic analysis of explicit expressions, and in the case of fully accounted for multiple scattering, analogous closed-form expressions are not available. To deal with this issue, and obtain a closed-form expression that yields a more accurate estimate of the convergence of the MEM in the case of closely spaced cylinders, we develop an approximation that accounts for first-order scattering effects while neglecting higher order multiple scattering effects. We then derive the rate of convergence of this approximation.
Denote by E (1) (N ) the approximation error associated with the first-order scattering approximation of the MEM just discussed. Since first-order scattering effects strongly dominate over higher order scattering effects, we expect that E (1) (N ) accurately characterizes the converge of the MEM, unless some of the cylinders are very close together. We have the following result for the convergence of E (1) (N ).
for a point source, for a plane wave. (1.2) Numerical simulations confirm that γ 2 (N ) is indeed a far more accurate estimator of the convergence of the MEM than γ 1 (N ) for the closely spaced case. When some cylinders are placed very close together, possibly almost touching, this approximation degrades somewhat as the higher multiple scattering effects become too significant to safely neglect. We discuss the possibility of accounting for higher multiple scattering effects by connecting our first-order scattering approximation with an iterative technique called the method of reflections [20,14]. Both Theorem 1.1 and Theorem 1.2 were derived for the case of a scattering problem featuring Dirichlet boundary conditions for which an indirect boundary integral equation solution representation was chosen. It is straightforward to show that if one were to consider a different set of boundary conditions, or a different boundary integral equation solution representation, different sub-exponential factors would arise during the derivations, but the expressions for γ 1 (N ) and γ 2 (N ) would ultimately remain the same. Thus, our theory holds not just for the setting we consider specifically, instead it holds for all boundary conditions and all boundary integral equation solution representations.
This paper is structured as follows. In Section 2, we introduce some notational conventions and function spaces. In Section 3, we describe the setting of the multiple scattering problem, and provide a representation of it in terms of an indirect boundary integral equation. In Section 4, we apply a MEM discretization to the boundary integral equation. In Section 5, we present our convergence analysis of the MEM along with numerical simulations that validate the theoretical findings.

Preliminaries. Denote by
Denote by Γ p := ∂Ω p the boundary of cylinder Ω p , and by Γ := M p=1 Γ p the boundary of the full system of M cylinders. Denote by (r p (x), θ p (x)) the polar coordinates of the point x ∈ R 2 with respect to a polar coordinate system whose origin is located at the center of cylinder Ω p ; see [36, Fig. 1].
The setting we use for the MEM formulation of a scattering problem is the space of periodic functions on M cylinders, which is a natural generalization of the standard space of periodic functions on [−π, π] This manuscript is for review purposes only.
Let f be a periodic function that has an expansion in terms of these basis functions: The coefficient f p m can be viewed as the m-th generalized Fourier coefficient associated with cylinder Ω p . Now, letting s ∈ R, we denote by H s (Γ) the fractional Sobolev space containing those f that satisfy The inner product on H s (Γ) is and the duality pairing on We define l s to be the space consisting of doubly indexed sequences (α p m ) p∈N M m∈Z that satisfy Thus, we have that f ∈ H s (Γ) if and only if its associated Fourier coefficients (f p m ) p∈N M m∈Z ∈ l s . As we primarly work in terms of Fourier coefficients in this paper, for convenience, we abuse notation and use the same notation for the norm of a function and the norm of its associated doubly indexed sequence of Fourier coefficients. Specifically, when we write ||(f p m ) p∈N M m∈Z || s , we mean ||(f p m ) p∈N M m∈Z || ls . Likewise, for the operator norm of an operator A : l s → l t , when we write ||A|| s,t , we mean ||A|| ls→lt .
For our MEM approximation of functions f ∈ H s (Γ), we introduce the finite-dimensional spaces T N (Γ) ⊂ H s (Γ) defined, for N ≥ 0, as We note the following identities where J m is the Bessel function of order m, and H m is the Hankel function of the first kind of order m: 3. Problem setting and boundary integral equation formulation. Recall that we are concerned with the scattering of waves by M disjoint circular cylinders {Ω p } M p=1 . Denote by Ω := M p=1 Ω p the set of all cylinders, and by Ω + := R 2 \ Ω the region exterior to the cylinders. Let k > 0 be the wavenumber in Ω + . We assume that k 2 is not a Dirichlet eigenvalue of −∆ inside Ω p , for p ∈ N M ; this condition ensures our problem is well-posed [18,Prop. 2]. To derive the asymptotic rate of convergence of the MEM, it suffices to consider the Helmholtz equation with Dirichlet boundary conditions for the total field u: in Ω + , u = 0, on Γ,

1)
This manuscript is for review purposes only.
where u s is the scattered field, and u inc is the incident wavefield given by for a point source. (3.2) Here, H 0 is the Hankel function of the first kind of order zero, and β = [cos(β), sin(β)] T , withβ ∈ [0, 2π], is the direction of propagation of the plane wave. Finally, we require that the scattered field satisfies the Sommerfeld radiation condition: where ∂/∂r is the radial derivative. Denote by G k the outgoing fundamental solution of the associated Helmholtz equation: For φ ∈ H −1/2 (Γ), we introduce the single layer potential S k defined as Upon taking the trace of the single layer potential, we obtain the single layer boundary integral operator We consider an indirect boundary integral equation representation for the solution of (3.1): in Ω + .
Upon taking the trace of this equation, we obtain the following boundary integral equation: where we denote by f := −u inc for convenience. Our aim is to ascertain a precise characterization of the asymptotic rate of convergence of the MEM applied to this equation. In the sequel we suppress the wavenumber dependence of V k for clarity and simply write V .
4. Spatial discretization with the MEM. To obtain a MEM discretization, first we have to represent the scattering problem in terms of infinite series of multipoles. Once this representation has been obtained, we truncate the infinite series to obtain a finite-dimensional discretized problem that can be solved numerically. We obtain a weak formulation of the multiple scattering problem by multiplying (3.6) by a test function and integrating over the boundary of the cylinders.
the problem becomes a matter of finding the doubly indexed sequence of coefficients (φ p m ) p∈N M m∈Z . Specifically, substituting the above expansions into (4.1) and using Galerkin orthogonality the problem becomes: Define the infinite dimensional per-cylinder coefficient vectors, for p ∈ N M , by Then the problem can be expressed in terms of infinite block matrices and vectors: Here, The elements of the matrices V pq are given by The elements V pq mn have explicit representations [36,18]: Here, θ pq is the angle between cylinder Ω p and cylinder Ω q ; see [36, Figure 1] or [31, Figure 2.1]. The incident field coefficients also have explicit representations [36,Prop. 4 and Prop. 5]: for a plane wave, On a historical note, representations (4.4) and (4.5) for V pq mn and f p m , respectively, appear as Equation (5) in Row's 1955 paper [34]. Equation (3) in Row's paper corresponds to (4.3). It is common to apply a diagonal preconditioner to the system of equations (4.3), as it vastly improves the performance of iterative solvers such as GMRES [28,36]. We apply the same preconditioner in this paper, however, our motivation for applying the preconditioner is rather to facilitate the convergence analysis in Section 5. In any case, we multiply both sides of (4.3) by and obtain the following equivalent problem: This manuscript is for review purposes only.
Here, G = BF , and W = BV = I + A, with The elements of the matrices I pp are given by The matrices A pq and vectors g p are given by Therefore, by (4.9), (4.6), and (4.5) we have: for a point source. (4.10) Likewise, by (4.8), (4.6), and (4.4) we have: So far we have just expressed the continuous problem (4.1) in a different form. The next stage in the MEM discretization procedure consists of truncating the infinite-dimensional block matrices and block vectors to obtain a finite-dimensional discretized problem. However, in order to perform a convergence analysis, rather than directly working with the finite-dimensional truncated objects, we use infinite-dimensional versions of them in which the elements that fall outside the truncation range are set to 0. Throughout this paper, we use tildes to denote the effectively finite-dimensional truncated MEM matrices and vectors associated with the infinite-dimensional matrices and vectors of the original problem.
We denote by N ≥ 0 the MEM truncation number. Definẽ Note that the various mathematical objects in the discrete problem have a dependence on N but we regularly suppress this in the sequel. Truncation of the matrices and vectors that comprise the block matrices and block vectors in (4.7) leads to the following discrete problem: Here,W =Ĩ +Ã, with The elements of the matricesĨ pp are given bỹ The elements of the matricesÃ pq are given bỹ (4.13) The elements of the vectorsg p are given bỹ (4.14) In practise, one numerically solves the linear system of equations in (4.12) to obtain the MEM approximate solutionΦ(N ). The approximation error E(N ) is the difference between the solution of the original problem Φ and the MEM approximate solutionΦ(N ): The question of precisely how fast E(N ) decays to zero as N → ∞ is a fundamental aspect of the MEM that has not been properly addressed to date in the literature and is the focus of the next section.
5. MEM convergence theory. Before we begin, we introduce some notation for the purposes of clarity. The functions involved in our MEM convergence theory frequently feature rates of growth or decay that are at least exponential with respect to some variable of interest. In light of this, algebraic factors in these functions are asymptotically irrelevant; they ultimately lead to arbitrarily small corrections to the rate of convergence. This motivates the introduction of the following notational convention, since it allows us to absorb algebraic factors. We use the notation a b, (5.1) to signify that a is bounded by b up to some function that increases sub-exponentially with respect to m. To be specific, we define a sub-exponentially increasing function as any function that increases with respect to m slower than e Cm , for C > 0, as m → ∞.
As an example of this notation let us consider the large order asymptotics of the Bessel function J m and Hankel function H m of order m as we will require these later. For m ∈ Z \ {0} and r > 0, the following super-exponential uniform bounds hold for some constants c J , c H , C J , C H ; see [1,Section 9.3] or [16]. The notation of (5.1) allows us to 'disregard' the algebraic factors and say Stirling's approximation will also be required later [33]. Stirling's approximation states that for m > 0, it holds that √ 2πm m+1/2 e −m ≤ m! ≤ em m+1/2 e −m .
Using the notation introduced above, we can write this as Upon taking norms and applying the triangle inequality, we obtain We can also immediately say Our plan in what follows is to derive the asymptotic bound for E(N ) given in Theorem 1.1 by explicitly estimating the right hand side of (5.5) as N → ∞. It transpires that this bound provides a tight characterization of the approximation error when the cylinders are not too close together, however, it becomes somewhat pessimistic when some cylinders are, in fact, in close proximity to one another. Thus, once we have obtained an asymptotic bound on E(N ) using (5.5), we will return to the closely spaced case and consider a first-order scattering approximation based on (5.4); this approximation allows us to derive the bound given in Theorem 1.2 which provides a more accurate representation of the convergence of the MEM in this particular regime.
Proof of Theorem 1.1. It is well-known that when k 2 is not an interior Dirichlet eigenvalue for −∆ inside Ω p , for p ∈ N M , the operator A is compact due to the positive distance between any two cylinders in the system, and hence the fact that I + A is invertible with bounded inverse follows by the Fredholm theory; see [18,Prop. 2], [28,Sec. 5.2] or [2,Thm. 2]. Recall from (4.13) thatÃ(N ) is simply A with the elements of its sub-matrices set to zero outside a finite range. Hence, for sufficiently large N , I +Ã(N ) is also invertible and we have for some positive constant C. In Lemma 5.1 and Lemma 5.7, we will establish that for a plane wave, In the case of a plane wave incident field, the bound on ||G −G(N )|| −1/2 decays super-exponentially as N → ∞, and therefore this bound will always be dominated by the bound on ||A −Ã(N )|| −1/2,−1/2 which decays merely exponentially. For an incident field due to a point source on the other hand, if the point source is sufficiently close to one of the cylinders, the bound on ||G −G(N )|| −1/2 can be larger than ||A −Ã(N )|| −1/2 as N → ∞, and thus it can't be neglected. Bearing this in mind, and substituting (5.6), (5.7), and (5.8) into (5.5) gives (1.1).
for a plane wave, for a point source. Recalling (4.14), we have Now, by the definition of |g p m | in (4.10), and using the relations in (2.3), In this case of a plane wave, by (5.10), this means for some sufficiently large constant C, as N → ∞. Then, upon absorbing the algebraic factor using (5.1), we obtain the result in (5.9) for the case of a plane wave. The result for the point source case can be obtained in a fashion, albeit using (5.11) instead of (5.10).
Before proving Lemma 5.7, we need the following two lemmas. This manuscript is for review purposes only.
Proof. First, define C pq := π 2 a p a q /4. Then, using the relations in (2.3), and the fact that H m−n (x) ≤ H m+n (x), m, n ≥ 0, (5.14) it is straightforward to show that Then, upon applying the uniform bounds for the Bessel and Hankel functions (5.2), and absorbing constant factors using (5.1), we have that for m, n > 0, In the next lemma, we derive an explicit representation of  Then, denoting by z pq := (a q /d pq ) 2 , we get To find a simplified expression for the series in the above expression, consider that where (a) n is the Pochhammer symbol, otherwise known as the rising factorial. This is essentially the definition of the hypergeometric function 2 F 1 for a particular set of parameters: 2 F 1 (m + 1, m + 1; 1; z pq ) = 1 + ∞ n=1 (m + 1) n (m + 1) n (1) n (z pq ) n n! .

Hence we have
To progress further we need a large argument asymptotic approximation of the hypergeometric function 2 F 1 for the above set of parameters. This can be found in Lemma A.1 in the appendix, where we prove that as m → ∞, it holds that Therefore, as m → ∞, we find that Remark 5.4. It is worth nothing that if one were to choose a different set of boundary conditions than the Dirichlet conditions specified in (3.1), or choose a different boundary integral equation solution representation for the problem, the result in Lemma 5.2 would not change. This is due to the fact that when we use the large order asymptotics of the Bessel and Hankel functions, the expressions obtained with those other choices differ from the expressions in Lemma 5.2 only by asymptotically irrelevant subexponentially increasing factors which are ultimately absorbed using the notation (5.1). It is for this reason that our convergence theory holds for all MEM formulations irrespective of the specific boundary conditions or boundary integral equation solution representation chosen. On a related note, recall that the specific boundary integral formulation we work with in this paper is an indirect single layer potential solution representation with a diagonal preconditioner applied. Interestingly, it has been shown [35] that every boundary integral equation shares the same spectral properties after applying the relevant diagonal preconditioner, which results in iterative Krylov subspace solvers such as GMRES having the same rate of convergence, irrespective of the particular boundary integral equation formulation chosen.
Remark 5.5. Results somewhat similar to the ones we derived in Lemma 5.2 and Lemma 5.3 were obtained in [12] in the context of spectral and condition number estimates of the single layer operator in dense media in the low-frequency regime, using an approach that doesn't involve hypergeometric functions. The fact that the low-frequency results in [12] are similar to our results based on large order asymptotics is not so surprising when one considers that (iπa p /2)J m (ka p )H m (ka p ) ∼ a p /2m, when either k → 0 while m > 0, or k > 0 while m → ∞, that is, the small argument asymptotics coincide with the large order asymptotics [16,Sec. 2.1]. See also [13] for related spectral and condition number estimates of the single layer operator, this time in dilute media.
The results established in the previous two lemmas are succinctly summarized in the following corollary which will be used in Lemma 5.7. Corollary 5.6. As N → ∞, it holds that This manuscript is for review purposes only.
Proof. By Lemma 5.2 and Lemma 5.3, and noting that 0 < a q < d pq for p, q ∈ N M , we have The result for the other series can be obtained in a similar fashion.
Proof. Since the operator norm is dominated by the Hilbert-Schmidt norm, we have We decompose the two inner series as follows: since we have from (4.13) that A pq mn =Ã pq mn for m, n ∈ Z N , andÃ pq mn = 0 otherwise. Consider, for a moment, the second series on the right hand side of the above expression. It holds that Substituting this into (5.16), and applying Corollary 5.6, we obtain m∈Z n∈Z Finally, substituting this result into (5.15), we find that In Figure 1, we provide convergence plots that demonstrate the accuracy of the bound γ 1 (N ) from Theorem 1.1, for an M = 3 three cylinder system, with radii (a 1 , a 2 , a 3 ) = (2, 1, 0.5), in the case of a point source incident wavefield, with the source located far away from the cylinders. The approximation error E(N ) was computed using the MEM scattering library µ-diff [36]. For each subplot, E(N ) is given by whereW is defined after (4.12), and N ref is taken five higher than the largest truncation number N used in the subplot. Due to the fast convergence of the MEM, this choice of N ref is sufficient forΦ(N ref ) to accurately approximate Φ. The block matrices and vectors used inΦ(N ) were zero-padded to make them align correctly with the corresponding block matrices and vectors used in the reference solutionΦ(N ref ). As sub-exponential factors are asymptotically irrelevant, we can disregard the (1 + |m| 2 ) −1/2 term in the fractional Sobolev norm l s (2.2) and perform l 0 based computations.
The plots in the first, second, and third columns correspond to closely spaced cylinders, moderately far apart cylinders, and far apart cylinders, respectively. With regards to wavenumbers, in each column: the first row features a wavenumber of k = 0.6; the second row features a wavenumber of k = 3; the third row features a wavenumber of k = 15. These wavenumbers have been so chosen because they imply regimes in which the wavelength is: large with respect to the diameter of the mid-sized cylinder (first row); around the same size as the diameter of the mid-sized cylinder (second row); smaller than the diameter of the mid-sized cylinder (third row). Hence we are analyzing the performance of the bound in a range of representative settings. Note that in the second and third columns, the last row is missing. This is because a very large truncation number N is required to reach the asymptotic regime in the case of the large wavenumber k = 15, which results in numerical precision issues.
It is clear that γ 1 (N ) accurately characterizes the convergence of the approximation error E(N ) in the moderately far apart and far apart regimes. However, it leaves something to be desired in the closely spaced regime; in this case, it is overly pessimistic. To obtain a more precise characterization of the approximation error in this regime we need to derive an estimate based on (5.4) rather than (5.5), since the former expression takes account of the multiple scattering of the incident wavefield among the cylinders. But this expression is problematic since it features the full solution of the untruncated problem Φ on the right-hand side, a term for which a closed-form expression does not exist. To overcome this difficulty, in the next section we develop a first-order scattering approximation which allows for a more accurate characterization of the convergence of the MEM in all regimes, but particularly in the closely spaced regime in which γ 1 (N ) becomes overly pessimistic.

5.2.
A first-order scattering approximation for closely spaced cylinders. Recall that the MEM formulation of the scattering problem involves finding Φ such that Note that if we disregard the operator A, we simply have Φ = G. This is essentially a MEM problem in which the incident field scatters off the cylinders but the subsequent multiple scattering interactions among the cylinders are neglected. As we are already using Φ to represent the full solution, let us denote byΦ the solution to this single-scattering problem: This means the original solution can be decomposed as Φ = G + Φ diff . Substituting this expression into (5.4) and applying the triangle inequality gives We can write this as  Here, E (1) (N ) is the component of the bound (5.19) associated with the first-order scattering event, while E (diff) (N ) is associated with subsequent multiple scattering interactions. In some sense, the term ||(A −Ã(N ))G|| −1/2 in E (1) (N ) can be viewed as a measure of the approximation error associated with the incident wavefield striking cylinder Ω p and being transferred once to cylinders Ω q , for p, q ∈ N M with q = p. This is in contrast to the term ||(A −Ã(N ))Φ|| −1/2 in (5.4) which can be viewed as a measure of the approximation error associated with the incident wavefield striking the cylinders and undergoing an infinite number of multiple reflections among them. As long as the cylinders are not close enough together such that the effect of multiple scattering becomes comparable to the initial single scattering effect, E (1) (N ) provides the dominant contribution to the bound (5.19). Hence, the bound on E(N ) in (5.4) is well approximated by E (1) (N ) which accurately characterizes the convergence. Thus, our plan now is to explicitly derive the rate of convergence of E (1) (N ). Numerical simulations, which we present later, verify that E (1) (N ) does indeed describe the convergence of E(N ) in the closely spaced regime in which the bound from Theorem 1.1 becomes overly pessimistic. Of course, if some cylinders are brought very close together, this first-order scattering approximation breaks down, since in that event E (diff) (N ) can become quite significant. In any case, we will shortly derive the rate of convergence of E (1) (N ) for the cases of point source and plane wave incident wavefields.
First, however, we make a brief remark about how the approximation we have just described can be connected to the method of reflections [20], which is also known as the boundary decomposition method [14]. This approach amounts to treating the scattering of waves between objects in an iterative fashion. In fact, there are several different types of methods of reflections, the one we consider below is known as the parallel method of reflections. Denote by Ω + p := R 2 \Ω p the region exterior to cylinder Ω p , for p ∈ N M . The scattered field in a Helmholtz multiple scattering problem can be decomposed as u s = p∈N M u s p where u s p satisfies the Sommerfeld radiation condition, and Now, in terms of the method of reflections, the single-scattering problem is associated with (5.20). In this problem only the initial single-scattering event is considered; the subsequent multiple scattering events are represented by (5.21). So the method of reflections problem (5.20) corresponds directly to the MEM single-scattering problem (5.18).
This connection between our first-order scattering approximation and the method of reflections suggests that it may be possible to characterize the higher order multiple scattering effects in the MEM in an iterative fashion, similar to how (5.21) iteratively provides the higher order effects for the method of reflections. However, this is not entirely straightforward, as there are certain conditions that need to be met so that the method of reflections series solution converges [14,25]. A means of overcoming the convergence issue could be to use the averaged parallel method of reflections (APMR) introduced in [29], which can be viewed as a relaxtion applied to the standard parallel method of reflections that leads to a convergent series. We now derive the convergence of the first-order scattering approximation E (1) (N ).
Proof of Theorem 1.2. The proof is practically the same as that of Theorem 1.1, the difference being that instead of estimating ||A −Ã(N )|| −1/2,−1/2 , we have to estimate ||(A −Ã(N ))G|| −1/2 . This can be accomplished in a similar manner to how ||A −Ã(N )|| −1/2,−1/2 was estimated in Lemma 5.7, so we only highlight the key differences. It is straightforward to show that The inner summations can decomposed as m∈Z n∈Z (5.23) We can bound m∈Z c N n∈Z (1 + |m| 2 ) −1/2 |A pq mn | 2 |g q n | 2 as outlined in Lemma 5.2, in the process also absorbing the algebraic factor (1 + |m| 2 ) −1/2 , to get σ pq (m, n), (5.24) where, this time, σ pq (m, n) depends on the incident wavefield. First, we deal with the point source case, for which we have .
By the same approach used in Lemma 5.3, we find that as m → ∞, Denote by z pq := (eka 2 q /(2d pq )) 2 . Applying Stirling's approximation (5.3) to some of the terms in σ pq , just as we did in Lemma 5.3, we find that (m + n)! m!n! 2 1 (n!) 2 e 2n (z pq ) n , so we have once again reduced the problem of bounding ∞ n=1 σ pq (m, n) to that of finding a large argument asymptotic expansion of a hypergeometric function. In Lemma A.2, we show that, as x → ∞, it holds that 2 F 3 (m + 1, m + 1; 1, 1, 1; z) Exp(4(m 2 z) 1/4 ). Therefore, as m → ∞, we have where we absorbed the root-exponential term using (5.1). Once again, using (5.22), (5.23), and (5.24), we proceed along the sames lines as Corollary 5.6 to obtain that, as N → ∞, In Figure 2, we plot the convergence of the approximation error E(N ) of the MEM, along with the bound γ 2 (N ) on the first-order scattering approximation derived in Theorem 1.2, for the case of a point source located far from the cylinders; the plot for the plane wave case is very similar so we omit it. We also show the overly pessimistic bound γ 1 (N ) that was derived in Theorem 1.1. The setting in these plots is the same as in Figure 1, which we recall shows the convergence for the case of low (first row), medium (second row), and high (third row) wavenumbers when the cylinders are close together, a moderate distance apart, and far apart. It is clear that γ 2 (N ) characterizes the convergence of the MEM better than γ 1 (N ) in all regimes, with this improvement being particularly noticeable in the closely spaced regime.
Note also that γ 2 (N ) may in fact slightly overestimate the convergence in the closely spaced regime, which is to be expected since in this case, while the first-order initial scattering event has a dominating effect on the convergence, the higher order multiple scattering effects, which correspond to the neglected term E (diff) in (5.19), are also starting to become noticeable.
While the asymptotic rate of convergence of the MEM is wavenumber independent, it can be seen from Figure 2 that the range of values of N for which γ 2 (N ) accurately characterizes the convergence does have a dependence on k; as k increases, a larger N is required before the asymptotic regime is reached. This is because the functions |J m (ka p )| and |H m (ka p )| only begin to reach their asymptotic rates of convergence when m ≈ ka p . A similar phenomenon has been observed with the method of fundamental solutions [16], which is another meshless method featuring solution representations comprised of Bessel and Hankel functions.
An interesting question is that of when the first-order scattering approximation breaks down. Numerical investigations show that the approximation is accurate as long as the distance between the closest points on the largest cylinder and the mid-sized cylinder is greater than approximately 5/4 times the radius of the mid-sized cylinder, irrespective of the wavenumber. Equivalently, for the wavenumbers k = 0.6, 3, and 15 used in Figure 2, the first-order scattering approximation is accurate as long as the distance between the closest points on the largest cylinder and the mid-sized cylinder is greater than approximately 1/8, 5/8, and 3 wavelengths, respectively. If the distance between the cylinders is any less than this, an approximation that takes into account higher order scattering effects would be necessary to accurately characterize the convergence of the MEM. 6. Conclusion. In this work we have provided a resolution to the long-standing problem of characterizing the asymptotic rate of convergence of the approximation error of the MEM by performing a detailed convergence analysis. The system of equations we considered first arose in Row's 1955 paper [34,Eq. (3) and Eq. (5)]. We began by deriving a bound that is tight as long as the cylinders are not too close together. To handle the case when some cylinders are, in fact, in close proximity to one another, we formulated a first-order scattering approximation for the MEM approximation error. This approximation accounts for the initial scattering event, namely, an incident wavefield impinging on each of the cylinders which then gets reflected onto each of the other cylinders. Meanwhile, higher order repeated multiple reflections of waves among the cylinders are neglected. We derived explicit bounds on the rate of convergence of this approximation, for the cases of both point-source and plane wave incident wavefields.
While our estimates were derived based on an indirect boundary integral equation solution representation applied to a Dirichlet scattering problem, this was merely for convenience. The convergence of the MEM for other boundary conditions or boundary integral equation solution representations differs from the case we considered only by sub-exponentially increasing factors which are asymptotically irrelevant. Thus, ours is a general theory of the MEM convergence that holds for all boundary conditions and boundary integral equation solution representations.
While the primary aim of this paper was to address the long-standing question on the asymptotic convergence of the MEM, there are several avenues worthy of investigation in terms of future research. Firstly, one could explore the generalization of the approach outlined in this paper to the case of spheres in three dimensions. The three dimensional case is more complicated as the solution representation features not only Bessel/Hankel functions and complex exponentials, but also associated Legendre functions. Moreover, more complicated addition theorems are required [31]. Due to extra difficulties such as these, the hypergeometric functions that could potentially arise during the analysis of the three dimensional case may turn out to be too complicated to analyze asymptotically. If this is the case, it would be interesting to investigate whether the approach employed in [12] could be of use, as in this paper expressions somewhat similar to ours were derived without making use of hypergeometric functions.
One could also attempt to provide a more accurate characterization of the convergence of the MEM in the case of cylinders that are almost touching; we conjectured that this may be possible by connecting our approach with a technique known as the method of reflections. Finally, since the MEM often features as a building block in other numerical methods such as, for instance, the MEM-based lattice summation techniques that arise in the field of photonic and phononic metamaterials [9], the framework we outlined in this paper could also prove helpful in obtaining rates of convergence in those approaches.
7. Acknowledgements. The authors wish to thank the reviewers for their insightful comments and suggestions that helped improve and clarify this manuscript.