Asymptotic behavior of acoustic waves scattered by very small obstacles

The direct numerical simulation of the acoustic wave scattering created by very small obstacles is very expensive, especially in three dimensions and even more so in time domain. The use of asymptotic models is very efficient and the purpose of this work is to provide a rigorous justification of a new asymptotic model for low-cost numerical simulations. This model is based on asymptotic near-field and far-field developments that are then matched by a key procedure that we describe and demonstrate. We show that it is enough to focus on the regular part of the wave field to rigorously establish the complete asymptotic expansion. For that purpose, we provide an error estimate which is set in the whole space, including the transition region separating the near-field from the far-field area. The proof of convergence is established through Kondratiev’s seminal work on the Laplace equation and involves the Mellin transform. Numerical experiments including multiple scattering illustrate the efficiency of the resulting numerical method by delivering some comparisons with solutions computed with a finite element software. 1 General introduction Mechanical wave simulations are of great interest in many applications due to their capability of transporting information in the medium they propagate into. In particular, they are capable of detecting very small heterogeneities that we can visualize using the recording of the scattered fields generated whenever the wave field encounters obstacles [28, 29, 30, 8]. This phenomenon of multiple diffraction can be reproduced numerically and for precise calculations, the finite element method is very efficient. Indeed, finite element methods are able to capture the characteristics of the obstacles including their shape by the use of tetrahedral meshes. However, if the obstacles are very small, their use can lead to very high computational costs, especially because it is necessary to mesh very finely in the neighborhood of the obstacles. This is a considerable drawback which

Most of the time, the use of explicit time schemes is preferred to integrate wave equations. This choice is dictated by the need of avoiding inversion of the mass matrix at each time iteration, in order to limit computational costs as well as memory usage. In this case, the explicit schemes are not fully suitable for direct numerical simulation. Indeed, they are stable only by respecting the Courant-Friedrichs-Levy condition which shows a linear dependency of the time step with the space step. Globally applying this condition to a direct numerical simulation based on a locally refined grid leads to increased computational costs since the time step will be adapted to the smallest cell. This is particularly disabling when the domain has only a few small obstacles. And even if the domain contains a lot of them, it should not be forgotten that explicit time schemes can generate numerical pollution mainly due to dispersion effects caused by the application of a global time step which is only adapted to the smallest cells. This has motivated the development of time-explicit schemes using local time steps that can be locally adapted to space steps (see e.g. [16,38].) Local time stepping methods are particularly interesting when the number of small cells is high. When this is not the case, it may be more efficient to use locally implicit schemes (see e.g. [34]). These schemes result from coupling an unconditionally stable implicit scheme applied in the areas of the mesh composed of small cells with an explicit scheme elsewhere. One can also take another side and decide to work on the mathematical model itself in order to keep an explicit scheme for time integration. The work to be done on the equations then consists in erasing, as it was, the small obstacles so that they do not have to be taken into account by the mesh. The latter will no longer contain small cells and we will thus be able to integrate the equations with a larger time step. To build such a mathematical model, asymptotic methods turn out to be very efficient besides having a rigorous framework for establishing convergence results which ensure that the solution calculated by solving the asymptotic model converges towards the solution of the initial problem. Asymptotic methods have been applied several times to stationary problems (see [13,21,23,26,40,42]) and to the best of our knowledge, Mattesi and Tordeux [32,33] represents the first attempt in the time domain.
The construction of a reduced model for representing scattering problems relies on a key procedure which consists in matching the asymptotic expansions of the near and far field in order to get a full representation of the scattered field. We describe that procedure in the following section. This work includes two new contributions in addition to being carried out in the time domain and in three-dimensional space. On the one hand, we extend the computational method proposed in [32,33] to the case of several small obstacles. On the other hand, we develop a mathematical analysis which justifies the asymptotic model proposed in [32,33] for the case of a single small obstacle.
The paper is organized as follows. We begin with describing the matching procedure for the construction of the asymptotic model. This requires defining the near-field and far-field approximations which are then matched to give a representation of the diffracted field throughout the space, disregarding the obstacle which should no longer be taken into account in the computational method. Then we extend the asymptotic model proposed to the case of scattering by multiple small obstacles. We illustrate the interest of our approach by carrying out several numerical experiments which consist in solving the problem of multiple scattering by using either the asymptotic model or the complete model. The latter is solved with a discontinuous finite element method. All the calculations we have performed show that the numerical solution obtained from the asymptotic model achieves a level of precision comparable to that obtained with the finite element method when we consider small obstacles. In addition, as expected, the computational costs are clearly to the advantage of the asymptotic method, which moreover does not require any mesh and thus avoiding a step that can be difficult and often time-consuming. For example, in the case of 216 obstacles, the calculations on the reduced model take 2.56 s while the direct simulation lasts 12 h for a computation carried out in parallel on 576 processors. We finish our study by developing a mathematical analysis to rigorously justify the proposed reduced model. The analysis consists in establishing a convergence result which uses the Mellin transform (see [7]) in the formalism of Kondratiev spaces [24].

Description of the matching of asymptotic expansions
To simplify the introduction of the matching method, we limit ourselves to the case of a single obstacle defined as the small sphere with a near-zero radius . We will generalize the method to the case of several obstacles by the end of the paper. The corresponding scattering problem posed in Ω = R 3 ∖ reads as: is the spatial variable defined in R 3 , ≥ 0 denotes the time variable and ∆ = 2 + 2 + 2 is the Laplace operator in the cartesian coordinate system.

Asymptotic parameterization
The construction of the approximation of the diffracted field involves a certain number of physical parameters which will intervene in the characterization of the wave field. First of all, we introduce the parameter which defines the radius of the obstacle. It is small in the sense that it is very small in front of the mean wavelength of the wave field. The latter is defined as the distance between two maximum values of the Fourier transform of the signal. The Fourier transform is most of the time centered around a characteristic angular frequency which defines a characteristic wavelength = where denotes the propagation velocity of the wave. By small obstacle, we then mean that is very small in front of . Now that the asymptotic parameters are defined, we move on the parameterization of the exterior domain surrounding the scatterer. For any x in the exterior of , we set: Then we have: x is located in the so-called far field region where the propagation phenomenon is preponderant while the obstacle has little impact on the scattered field. The scatterer acts like a sourcepoint so that the scattered field can be represented by an asymptotic expansion set in 3 .
-if = 1, x = X; x is located in the so-called near-field region where the obstacle has strong impact on the scattered field while the propagation phenomenon is negligible. A quasi-static behavior of the scattered field is actually expected from its asymptotic expansion. -if 0 < < 1, x is in the so-called matching area where both the propagation and the obstacle effects have to be taken into account. Hence this is here where the far field expansion has to match the near-field one to ensure that the final asymptotic representation provides a reliable representation of the scattered field. For this, it is crucial to see that if x is of order ℓ 1− then x/ is of order −1 ℓ 1− . The first one tends to zero as tends to zero whereas the second tends towards infinity.
Remark 2.2. It is worth noting that it is which acts as the small parameter. Hence in the following, the scattered field will be given as a series involving powers of .

The far-field approximation
In the far-field region, we seek a solution to problem (2.1) in the following form ( Fig. 1): 3) Figure 1. Far-field and near-field domains Ω f and Ω n and their limit Ω ⋆ andΩ.
where each term : R 3,* × R + −→ R is a space-time function defined in the domain R 3 ∖ {0}. Following the methodology described in [9], we have: -The first term 0 is ∞ and is defined as a solution to the wave equation set in the whole space (i.e. without obstacle) (2.4) The function 0 actually represents the regular part of the wave field and is given thanks to the theory of separation of variables as The functions , are the spherical angular harmonics and are given by where is the Legendre polynomial of degree given by More details on Legendre polynomials are available for instance in ( [37], see p. 35), in ( [21], see p. 47) or in ( [19], see p. 353). The radial function , has the following expression where Λ 0, , : R −→ R is the amplitude of the regular mode , which takes the form (2.11) We adopt the notation R 3,* to refer to the space R 3 free of origin. Each term is represented as

The matching procedure
To make the asymptotic representation of the scattered field self-content, we have to realize the matching between both developments previously introduced. For that purpose, we go back to the relation that links the variables x and X. Given that x = X, we have that x tends to be close to 0 while X tends to infinity. This leads us to consider an expansion of each in the vicinity of the origin as well as a development of in a neighborhood of infinity. For this, we use the spherical coordinate system and write down: Now we move on to the explicit writing of the first terms of (2.26). In the vicinity of the origin, the far-field term 0 is a regular solution to the wave equation (2.4) and we have By identifying the different terms, we get that 0,0 ( , , ) = 0 (0, ) and 0,1 ( , , The matching conditions (2.23) and (2.25) then imply that 0,0 ( , , ) = 0,0 ( , , ) = 0 (0, ) and 0, ( , , ) = 0 for > 0.
In particular, Regarding 0 , it is a solution to the Laplace equation (2.17) which admits the following modal expansion [35] Taking (2.28) into account and given that 0,0 = 1 √ 4 , we deduce that 0 = √ 4 0 (0, ) and = 0 for ≥ 1. It follows from (2.19), that = − . We then get Next, we compare (2.31) with the term 1 in (2.12). It is worth noting that the term , (Λ 1, , , ) in (2.12) is computed in the vicinity of x = 0 according to the representation formula: Then, remarking that Λ 1, , ≡ 0, barring ( , ) = (0, 0), we get It follows that Λ 1,0,0 ( ) = −(4 ) 3/2 0 (0, ) and consequently we obtain Similarly, we get the first order near-field term while the second order far-field term is given by We then end up with the second order far-field expansion: where 0 is the regular solution to the wave equation in the whole space (2.4) (i.e. without obstacle). Moreover, we have This last expression can be injected in (2.37) to get a second order approximation of the wave field: A second order far-field approximation of the exact solution is given by Remark 2.4. This formal result can be made rigorous by an error analysis, [40]. Given that ( = ,2 + ( 3 )) and ,2 = mod ,2 + ( 3 ), we obtain that for all x ̸ = 0, there exist x > 0 and x > 0 such that for all ∈]0, x [.
It is worth noting that (2.37) and (2.39) do not involve the obstacle, which explains why using such representation of the scattered field will not require considering the obstacle as a geometrical object. By this way, there will be no need to reproduce the surface of the scatterer by a set of points, the latter being modeled as a point source.

Application to multiple scattering
We consider the scattering problem (2.1) created by spheres of radius and center x ∈ R 3 , 1 ≤ ≤ .
The propagation domain is denoted by The total field is decomposed as the sum of an incident field 0 defined by (2.4) and a finite superposition of scattered fields : There are different ways for computing an approximation of the solution of this problem. The field can be approximated by a mono-polar source of amplitude Λ in the following way: where inc,n is the field illuminating the obstacle (x ). Then, the approximate solution depends on the definition of inc, . Here, we address the following possibilities: -Adopting the first order Born approximation, the interactions between obstacles are neglected. In that case, This leads to an explicit formula of from: -The Foldy-Lax model, written classically for harmonic waves [9,11,12,18,31], includes the interactions between obstacles. Here for a given obstacle (x ), the incident field is the sum of the incident wave 0 and the fields scattered by all the other obstacles: with the corresponding amplitude computed according to the definition (3.2) Numerically, the unknown functions ↦ → Λ ( ) will be approximated on a regular grid in time: = with ∈ N. In the following, we denote by (Λ ) the numerical approximation of Λ ( ). It is worth noting that the numerical solution of system (3.5) requires some evaluations at time − , − for all ̸ = . These are discrete times which may not be on the time grid. An interpolation technique is then used to evaluate the right hand side with respect to the values on the grid.
In the above expression,̃︀ Λ is an interpolation of Λ given by In the following numerical study, we will compare results obtained with -The Born approximation.
-The Foldy Lax model. -A Direct Numerical Simulation performed with the numerical Library Hou10ni [6]. This software package is able to simulate time-dependent acoustic wave propagation in 3D. The space discretization is based upon the Interior Penalty Discontinuous Galerkin method [2,4,5,15,20] and the time integration is carried out with a Leap Frog scheme. There is a need in truncating the computational domain and for this, we use a first order absorbing boundary condition which is set on the exterior boundary. The code integrates the -adaptivity option, which allows to use different orders of approximation ranging from = 1 to = 6. In this way, the computational costs of the direct simulation method are minimized by adopting a refined mesh in the vicinity of obstacles and larger cells elsewhere.
In each numerical experiment, we consider an incident wave which is analytically defined by The wave speed is equal to 1. The real numbers 0 = 5 is a characteristic frequency whereas 0 = 1.2 0 is a characteristic time. Most of the energy of the signal corresponds to the wave length: First numerical experiment: the case of five spheres. We consider two configurations with five spherical obstacles with radius = 0.1 and = 0.01 centered in In terms of ratio of characteristic length, the ratio / is equal to 0.13 in the first configuration and to 0.013 in the second configuration.
For the first configuration, the direct numerical simulation has involved 88 073 elements distributed as follows 108 2-elements, It is worth noting here that reducing the diameters of the sphere by a factor 10 only slightly increased the number of elements, since we only need to small elements close from the obstacle. However, this divided the time step by a factor 20.
The results are depicted for the first configuration in Figure 2. We notice that the numerical solution of the Foldy-Lax model is much closer to the direct numerical solution than the solution of the Born approximation. This is particularly true in the center of the computational domain where a diffusion phenomenon can be observed. This is related to multiscattering which is not correctly taken into account by the Born model.
For the second configuration, the results are depicted in Figure 3. In this case, the radius of the obstacles and consequently the ratio / are smaller. Both approximations give excellent qualitative results.
Second numerical experiment: the case of 216 spherical obstacles evenly distributed. We consider here the case of a network of 6 by 6 by 6 spheres with radius = .01 whose centers x , , lie on a regular grid of space step ℎ = .5. More precisely, we for 0 ≤ , , ≤ 5,  This case is only affordable with direct numerical simulation performed on a supercomputer. On a laptop it would have taken almost two years. Since the simulation is performed in time-domain the memory issue is less crucial.
A clear difference can be observed close to the boundary of the computational domain. This is due to the approximate outgoing wave condition which does not filter the reflections completely, so we can see a spurious wave which pollutes the finite element calculation.
The Foldy-Lax model took 2.56 s to solve the problem on a personal computer (1.8 GHz simple core with 4 Gb of RAM). These preliminary results are presented in Figure 4. We can observe the solution at times = 2 and = 3. Only 36 spheres are visible since we are representing the solution on the plane = 0.25.
It is worth mentioning that thanks to Foldy-Lax model, we have been able to perform simulations of acoustic wave diffraction by dense arrays of small obstacles up to 10 000 in number. These simulations could not be carried out with the finite element code.

Statement of the result
In Section 2, we have formally written a matched asymptotic expansion providing a representation of the scattered field in presence of a small obstacle. In the following, we justify this development by proving that the Taylor expansion of in (2.20) converges. For > 0, the terms are singular functions and according to (2.12) and (2.13), they are represented by a finite sum (see [32,33] for details). As far as convergence is concerned, we will thus focus on the regular term  with It has been proved that (−1) = 0 for any odd positive integer number , see [32]. It follows the following Theorem The proof of (4.4) of Lemma 4.1 is simply related to an explicit computation. Details can be found in [32]. The rest of this section is devoted to the proof of (4.5) of Lemma 4.1.

Prerequisite for the mathematical analysis
The key tool for this work is

Mellin transform: definition and properties
In this sub-section, we summarize the classic Mellin transformation results, whose proofs have been gathered in [7] (easily accessible online). The reader may also refer to [37].
In what follows, is a real and is a positive integer. The Kondratiev spaces (see for instance [24], [10] and [14]) are defined by where the notation · {ℓ} stands for the differential operator These spaces are equipped with the Hilbertian norms (4.12) We denote by ∈ C complex number = + , ∈ R and ∈ R, Using the density of (]0, +∞[) in , the Mellin transform can be extended to . Moreover, the Mellin transform ↦ → ℳ ( ) is defined for any ∈ and for any ∈ C where C denotes the line in the complex plane On the other hand, the Mellin transform is an isomorphism from ontoˆ, whereˆstands for In the above definition, 2 (C ) denotes the space of functions on C equipped with the norm and satisfies for ℛ ( ) = < 0 Proposition 4.7. Let be a positive integer, 1 and 2 be two real numbers such that 1 < 2 . Let [ 1, 2 ] := 1 ∩ 2 . Then for any ∈ [ 1 , 2 ], we have ⊂ [ 1,2] .
To prove Theorem 4.2 and thereby estimate the rest u , we are going to apply the following theorem coming from singularity theory, that is: Theorem 4.8. Let 0 < 1 < 2 be three real numbers. Let ∈ 0 1 and C [ 1,2] be the strip of the complex plane defined by We assume that (i) the Mellin transform ℳ : C 1 → C admits an analytical continuationˆin C ] 0, 2[ ; (ii) there exists > 0 such that for any = + ∈ C ] 0, 2[ and | | > 1 then for any ∈] 0 , 2 [, ∈ 1 , and satisfies

Energy estimates
In this paragraph, we recall some elementary properties of the wave equation such as energy conservation or the principle of finite velocity propagation. Let us consider the solution 0 ∈ ∞ (R 3 × R + ) to the Cauchy problem (2.4) defined from a pair of initial data ( 0 , 1 ) satisfying Hypothesis 2.1. The wave field 0 propagates with finite speed and its energy is preserved (see [39], Chapt. 2, Sect. 6), that is From (4.26), we deduce that: For any pair of Cauchy data ( 0 , 1 ) satisfying Hypothesis 2.1, we have for all ≥ 0 Proof. Since the initial data are in ∞ , any time derivative = 0 satisfies the wave equation (2.4). For ≥ 1, we can thus apply (4.26) to Property 4.10. For any pair of Cauchy data ( 0 , 1 ) satisfying Hypothesis 2.1, the support of (·, ) is enclosed in ⋆ + and Proof. From the finite speed propagation principle, we have that 0 (x, ) = 0 except if |x| ≤ ≤ + ⋆ . The Poincaré inequality implies that Thus, the previous inequality extends to 3 : To complete the proof of Property 4.10, all that remains is to combine (4.26) with (4.31).

Spectral decomposition
In what follows, we recall some classical results about separation of variables for writing regular solutions of the wave equation (see e.g. [17,22,36]).
Any regular solution of the acoustic wave equation can be decomposed into This writing is valid for any value of and the series converges in 2 (R 3 ). The terms , are the spectral coefficients of 0 in the orthogonal basis , of 2 (S) previously introduced in (2.6); we have: where S denotes the unit sphere parameterized with the angular space variables ( , ). It is then important to list the properties of the functions , . We have: Proof. By construction, the set of functions , forms an orthonormal basis of 2 (S). Indeed according to (2.6), we have The corresponding Legendre polynomials are orthogonal and satisfy for any ∈ [0, ] (see e.g. [37] p. 37) By applying the change of variable = cos , we get the normalized Legendre polynomials which satisfy This ends the proof of (i). Still by definition (see (2.7)), (1) = 0, ∀ ̸ = 0. We thus have (ii) according to (2.6). By definition (see (2.6) and (2.7)), we have because (1) = 1. Indeed, following [37] (see p. 35), we have and by using the recurrence relation (4.38) for = 1, we obtain (1) = −1 (1). It is then sufficient to see that 0 (1) = 1 for getting (iii).
By orthogonality, we have for all > 0 which gives (iv) from (4.39). This completes the proof of Property 4.11.
In [27], the basis functions , of 2 (S) are defined as eigenfunctions of the Laplace-Beltrami operator which reads in spherical coordinates as We actually have The partial differential equation can be diagonalized by using the spectral coefficients , . According to (4.41) together with the fact that the functions , are orthogonal in 2 (S), the coefficients , satisfy for any > 0 and ≥ 0 We then have Proposition 4.13. The space 2 (R 3 ) can be characterized as By switching sum and integral terms (following Fubini theorem), it holds that Then, considering that

Convergence proof
We remind that u is given by the modal expansion (4.3). It is represented by a series to which we want to give meaning. We begin by studying the terms of this series precisely and then we deal with the convergence of this series in the classical sense and not only in 2 (R 3 ). We follow a two-step approach: -The first step will consist in estimating precisely the spectral coefficients , for every ∈ N and − ≤ ≤ .
-The second step consists in processing the infinite summation and demonstrate its convergence.

Control of the terms of the series
Here we will apply the Theorem 4.8 . After testing hypotheses (i) and (ii), we will obtain the behavior close to =0 of , which is described in the following proposition: Proposition 4.14. Let be an even integer greater or equal to 2. For any integer ≥ , ∈ N with − ≤ ≤ , we have  where denotes the integer such that − 2 ∈ [− 7 2 ; − 3 2 [. It is explicitly given by = ( 7 4 + 2 )) with = ℜ( ). The formula (4.59) makes it possible to defineˆ, (·, ) as a meromorphic function of C.

Lemma 4.15. For any
where is a function of R + × R −→ R + locally bounded in .
Proof. We remind that = + . We first remark that formula (4.59) is also valid for = 2 Considering the imaginary part of , we have .

End of the proof of Theorem 4.2
Let S be the unit sphere. We establish preliminary results dealing with the Laplace Beltrami operator. we have ∈ ∞ (S) and the following estimate holds: Proof. We setŷ ∈ S and we introduce a coordinate change that associatesẑ ∈ S tox ∈ S with Next, we note that and we conclude that which ends the proof thanks to (4.78).
Lemma 4.18. Let ≥ 2 be an even integer. We have Proof. According to Parseval equality (4.46), we get From Proposition 4.14, we deduce that (4.86) We apply Parseval theorem once again to get  This estimate holds true for any > 0 and we deduce a non-optimal estimate of u for any even integer set in the following proposition.

Conclusion and perspectives
We have proposed a new solution methodology for solving 3D multiple scattering problems when the size of the obstacles is small with respect to the characteristic wavelength. This work contains two key results. First, it validates an asymptotic representation of the field diffracted by a small obstacle illuminated by an incident acoustic wave of very large characteristic length in front of the radius of the obstacle. The proof is based on Kondratiev's theory and the extensive use of the Mellin transform. To the best of our knowledge, this result is the first one addressing the case of time-dependent wave problems. Second, this work shows the potential of the proposed asymptotic model to numerically simulate the effects of a large number of small obstacles on an incident wave. The proposed asymptotic method is confronted with an advanced direct simulation method based on discontinuous finite elements, the time integration being performed with a leapfrog scheme. In a rather simple case where the number of obstacles is limited to 5, the asymptotic method is validated by comparison with the finite element method. Then we treat the case of 216 obstacles to show that the asymptotic method continues to deliver an accurate solution for a very short computation time while the finite element method reaches its limits. It would also be interesting to consider the case of penetrable obstacles. A more complicated case would be to treat the case of elastoacoustic interactions in order to quickly detect small defects in large structures. It would also be of great interest to compare our solution methodology with those involving highorder radiation conditions [1,41] and to investigate possible extensions of [3] in the time domain. Indeed, now that we dispose of an accurate representation of the scattered field, we should be able to construct an OSRC in the time domain. Finally, the method we have proposed can be extended to other wave equations. Recently, it has been developed for Maxwell's equations in harmonic regime [25,26]. The time-dependent case is clearly more technical but it is quite possible, at least when considering low order asymptotic models.