ASYMPTOTIC ANALYSIS OF THE SCATTERING PROBLEM FOR THE HELMHOLTZ EQUATIONS WITH HIGH WAVE NUMBERS

. We study the asymptotic behavior of the two dimensional Helmholtz scattering problem with high wave numbers in an exterior domain, the exterior of a circle. We impose the Dirichlet boundary condition on the obstacle, which corresponds to an incidental wave. For the outer boundary, we consider the Sommerfeld conditions. Using a polar coordinates expansion, the problem is reduced to a sequence of Bessel equations. Investigating the Bessel equations mode by mode, we ﬁnd that the solution of the scattering problem converges to its limit solution at a speciﬁc rate depending on k .


1.
Introduction. In this article we study the scattering problem for the Helmholtz equation with high wave numbers, outside a circle. There is an abundant literature about the Helmholtz equations with high wave numbers, some of which is recalled below. Much less is known about studying the scattering in a bounded domain in particular at the level of numerical simulations, see the bibliographical comments below. Our interest in this problem arised from its relevance to the physics of lasers, where one is lead to study the scattering of lasers, inside a box, outside a box, or in a box with aperture (inside -outside problem); see e.g. [43]. Physicists routinely construct such "boxes" and make various measures that they would like to compare to results of calculations, in view, as in many areas of sciences and engineering, of replacing lengthy and costly laboratory experiments by numerical simulations. This article is a first very preliminary step in this direction. On the mathematical side, we have been guided in this direction by the works of [2], which  display the mathematical complexity of the asymptotic expansions that can be made at high wave numbers in the presence of a wall. In these references, as here, the applications include acoustics and electromagnetics (lasers), and such uncommon phenomena as the "whispering walls". We have been also guided by earlier works on singular perturbations and correctors techniques as explained below; see e.g. [24,25,26,29].
In this article we consider scattering by a circular object, of an incident wave, traveling to the right along the x axis (see Figures 1,2). Hence we consider the Helmholtz equation where D is the unit disk centered at (0,0), k is a sufficiently large wave number, and u (k) is the total field such that u (k) : R 2 → C. On the surface of the obstacle D, we impose the homogeneous Dirichlet boundary condition on ∂D, which corresponds to a sound-soft surface. It is well-known that the solution u (k) consists then of two parts (see Figure 2) where u (k) inc is the incident wave, and u (k) s is the scattered wave (for more details, see e.g. [2,14,39,49]). The incident solution, which satisfies (1), is given by The scattered solution also satisfies (1), and where (r, η) stands for the usual polar coordinates in the plane. We easily note that inc + u (k) s = 0 on ∂D.
Hence, the scattering problem reads as in [11,15,16,28,45,46] The conditions in (6) 3 are called the Sommerfeld radiation condition; see e.g. [47,48]. The physical meaning of the radiation condition is that the scattering field propagates away from the obstacle. Mathematically the condition ensure the uniqueness of the solution of certain exterior boundary value problems; see e.g. [14,47,48]. In mathematical physics, the Helmholtz equations often arise in the study of applied optics problems involving partial differential equations. The fundamental and classical results of solutions to Helmholtz problems can be found in [47,48]. The well-posedness of the scattering problems was carried out using a variational method in [11,12,36]. In particular, utilizing Green's functions, the authors have developed a priori estimates in [16,27]. Recently, sharp estimates were derived using Rellich's type identities in [15].
In numerical simulations, the computational domains obviously have to be bounded. One can consider absorbing boundary conditions such as the Dirichletto-Neumann map to handle the exterior boundary conditions; see e.g. [7,18,23,32] among many other articles. However, in this article, we propose a physically relevant boundary condition, which directly reflects the Sommerfeld radiation condition (6) 3 ; see e.g. (11) below.
The solutions to the Helmholtz equations corresponding to high wave numbers are highly oscillatory and the numerical simulations are very challenging because the conventional computational schemes require a very large number of degrees of freedom to represent the solutions. Given the number of degrees freedom per wavelength λ = 2π/k and per dimension d the computational cost is O(k d ). There have been efforts to reduce the number of degree of freedom. For instance, one can reduce the dimension using e.g. the boundary integral equations or the basis functions which are solutions to the Helmholtz equations like in the methods of particular solutions or non-polynomial methods. The basis functions include e.g. Fourier-Bessel, Fourier-Hankel, and fundamental solutions (see e.g. [4,5]). This leads to a computational cost O(k d−1 ) which is still large. In Figure 2 some solutions to the Helmholtz equations are simulated using the methods of particular solutions (see [6]).
To substantially reduce the cost, the so-called numerical asymptotic methods have been actively developed. See e.g. [10,21,22]. The methods require only a fixed computational cost for any k, i.e. O(1). The numerical asymptotic methods incorporate the asymptotic expansions of high frequency parts of solutions by making use of a priori knowledge regarding the oscillatory solution parts when designing the numerical schemes. The basis functions are then asymptotically derived. More precisely, they are the product of a slowly varying part and a highly oscillatory part (phase factor), which are determined a priori by geometrical optics (GO), and geometrical theory of diffraction (GTD) (see e.g. [30,31]). For instance, by an Ansatz we seek a solution of the form u(x, y) = K(x, y) exp(ikS(x, y)), s : total solution. This simulation is generated using the MPSPACK (see [6]).
where K is the amplitude, S is the phase and exp(ikS(x, y)) is the phase factor. Here, K, S are slowly varying. By the asymptotic expansions regarding the k, the function S is determined or derived associated with e.g. a reflected field, diffracted, etc. Accordingly, the basis functions have the form (7). The function K can be resolved numerically on a proper coarser grid by some standard numerical methods, e.g. conventional splines. Recently, combining boundary integral equations with asymptotic expansions the numerical asymptotic boundary integral methods are proposed (see the review article [10]). This can solve the scattering problems accurately in a computational time independent of the wave number k.
Asymptotic studies with high wave numbers have been performed in many articles. In [9,31,33], the authors studied the geometrical theory of diffraction in a local region by introducing the ray coordinates. Utilizing asymptotic expansions, the authors compared the known exact solution with ray representations to match solutions near a local region. In [3], constructing a specific Ansatz, asymptotic behaviors with high wave numbers were investigated in the local region. In addition, in [2,8,50], the authors have performed local analysis of various ray fields using different Ansatz depending on the regions. Around a convex obstacle like a circle, there are regions of shadow boundary near the point of ray tangency, near-surface region in deep shadow, and the shed creeping region. The asymptotic analysis using an Ansatz is performed in e.g. the illuminated region by geometrical optics, in the shadow region by geometrical theory of diffraction, and in the transitional shadow boundary region by Melrose-Taylor expansion (see [38]) as depicted in Figure 1.
The asymptotic analysis that we present in this article is based on the singular perturbation analysis and the boundary layer correctors technique which have been successfully applied to fluid related equations like the Navier-Stokes equations (see e.g. [17,24,25,29]). They are usually used for elliptic or parabolic operators. The asymptotic analysis performed here is new in that we deal with the Helmholtz equation. Namely, we formally define the limit solution u to the asymptotic expansion u (∞) +Θ in suitable function spaces (Sobolev spaces). On the other hand, we are able to use the exponential correctors in computations combined with e.g. finite element, finite volume methods or spectral methods. The exponential correctors are relatively easier to use than Bessel or Hankel functions, which are highly (extremely) oscillatory for a large wavenumber k. In the analysis below, for high modes |n| ≥ 2k , we indeed prove that the boundary layer correctors are of exponential type. Following the spirit of the numerical asymptotic methods, the Helmholtz solutions u can be decomposed into the sum of the oscillatory part u f and a slowly varying part u s , i.e. u = u f + u s . The former u f is approximated by the boundary layer correctors Θ and the rest u − Θ can be approximated by conventional numerical methods. The latter is considered slowly varying. The comparison of computations using the correctors with the existing computational methods is not simple; it will be made elsewhere.
Here we aim to study the asymptotic behavior of the solutions of (6) with high wave number k by introducing a corrector or correction function. We consider a reduced (limited) domain to investigate the scattering problem; see e.g. (9). The formal limit problem of (6), when k → ∞, reads We aim to verify that u (k) s converges to u (∞) s as k → ∞ in a proper space, and look for its explicit rate of convergence in terms of k. For the estimates of the scattering problem, we modify the variational method used in [45,46].
From now on and for the sake of simplicity, we use u to denote u (k) s . We rewrite the scattering problem in a reduced (limited) domain Ω as follows: where Here, L is a sufficiently large number, which may depend on k, so that the radiation condition (6) 3 implies where C is independent of k. Later, we will work on a fixed domain with the radius L 0 , which does not depend on k; see e.g. Theorem 3.1. The function g is implicitly given from the Sommerfeld condition (6) 3 on r = L. That is, g(L, η) := It thus satisfies the decaying condition (11) which will be 1164 DANIEL BOUCHE, YOUNGJOON HONG AND CHANG-YEOL JUNG used in the analysis below; for more details, see Section 3. We consider the polar coordinates expansions to study (9) as in [14,45,46]: For each n, (9) becomes where J n (k) is the Bessel function of the first kind. The boundary condition (13) 2 can be obtained using (4) and the Jacobi-Anger expansion (see e.g. [14]): We find by the Plancherel Theorem that Hence, we deduce from (11) and (15) that It is natural, to understand the scattering problem, to study (13) for each n.

Notations.
We define here, for later use, some weighted inner products and weighted norms: where s is an integer. For (r, η) ∈ Ω, we define where s is an integer andv is the complex conjugate of v. The weighted Sobolev where l ≥ 0 is an integer.
2. Asymptotic analysis. We study (13) in two different cases; the high mode case (|n| ≥ 2k ) and the low mode case (|n| < 2k ) where x which is called a floor function is given by Indeed, the sign of the term (n 2 /r 2 − k 2 ) ≈ (n 2 − k 2 ), near r = 1, in (13) is different where n varies from 0 to ∞. Hence, this requires to split the modes into two separate cases. We use different approaches for each of the cases since we can use exponentially decaying asymptotic expansion only for the high modes. For the low mode, due to the lack of decay (or the poor decay), we consider an explicit approach utilizing the Hankel functions and make use of their properties.
2.1. Low modes. We first consider the low modes, for which |n| < 2k . By settingr = kr, (13) 1 becomes We note that (21) is the Bessel equation. A fundamental solution of (21) is where H (1) n and H (2) n are the Hankel functions of the first and second kind, respectively; for more details on the Hankel functions, see e.g. [1,14,40,41]. In addition, C 1 and C 2 are constants. For each n, we introduce a so-called corrector function ϕ n which corrects the boundary conditions so that u n − ϕ n = 0 on ∂D. For alternate uses of the correctors, see e.g. [26] in fluid dynamics where the equations of the corrector play the role of the Prandtl equation (see e.g. [42,44]). We then construct, for each n such that |n| < 2k , the equations for ϕ n using (13) and (21): We first choose a solution of (23) 1 and (23) 2 without considering the boundary condition at r = L: Utilizing the solution (24), the Sommerfeld boundary value g * n (L) is determined: Lemma 2.1. For r > 1, we have where ϕ is as in (24) and C is independent of k and r. Here c 0 > 4 is a constant.
Proof. From Appendix A.4, we have the following asymptotic properties of the Hankel functions where |arg(z)| ≤ π − δ, δ is an arbitrarily small number, and z is complex. In addition, from Appendix A.5, we have for complex ρ, with |ρ| < r. Hence, if ρ = ic 0 and 4 < c 0 is a real positive constant, we derive

DANIEL BOUCHE, YOUNGJOON HONG AND CHANG-YEOL JUNG
where r > 1. Using (27) with n = 0, since kr is large, we find that For n ≥ 0, using the property in page 263 of [13], we deduce from (29), (30) and (31) that By the definition of the Hankel function, we easily see that Recalling the expression (24) of the correctors, we find By Stirling approximation (see Appendix A), we find 2 n n! (kc 0 ) n ≤ 2 n en n+1/2 e −n (kc 0 ) n = en We then deduce that ≤ C k Since c 0 > 4, summing the geometric series, we find We note that the case n ≤ −1 can be treated in a similar way using the properties of the Bessel functions listed in Appendix A.3, so that (26) is now proved.

High modes.
We now consider the high modes case, that is |n| ≥ 2k . It turns out that the correctors for the high modes possess an exponentially decaying property, and we can construct the error equations without handling special functions such as the Bessel or Hankel functions. Defining ε −2 := n 2 − k 2 in (13), we see that ε is small for large k, large n and that 0 < ε ≤ C/k. Indeed, we note that We introduce the variable ξ such that Here, we assume that L is sufficiently large, and as we said in Lemma 2.2, L = C 1 k α and α > 0. Utilizing the stretched variable, we rewrite (13) 1 as This implies We introduce the zeroth order corrector θ 0 n which corrects the boundary value at ξ = 0. Setting u n = θ 0 n in (54), we find the dominating differential operators at order O(ε −2 ) and establish the equations of θ 0 n : The explicit solution of (55) is easily found To achieve the main result of this article, we introduce additional correctors up to order O(ε 2 ). Setting u n = 4 s=0 ε s θ s n , (54) becomes For −1 ≤ l ≤ 2, we collect the corresponding terms in (57) at each order of ε l and establish the following corrector equations: where θ j n = 0 for j < 0. We impose the proper boundary conditions in (58) by using (13). For 1 ≤ j ≤ 4, we propose the corrector equations of θ j n : n + n 2ξ2 θ j−4 n , θ j n = 0, atξ = 0, θ j n → 0 asξ → ∞.

(59)
Here, θ j n = 0 for j < 0. An explicit solution θ j n is given in the following lemma.
Lemma 2.3. The corrector θ j n in (59) is given by where the A l j (ξ)'s are polynomials with respect toξ with coefficients independent of n, k, ε such that A l j (0) = 0. Proof. The proof is by induction on m. At m = 0, (60) for j = 1, 2 can be recursively verified by using the solution θ 0 n given in (56). We assume that (60) holds at order m ≥ 0. Then, at order m + 1, from (59) we have where |n| ≥ 2k . We write the Sommerfeld condition of θ n : for |n| ≥ 2k . From (54), (55) 1 , (59) 1 , and (62), we deduce that Regarding (13) and (64), and the radiation condition in (63), we find the equations of w n for |n| ≥ 2k : r ∂ r w n + n 2 r 2 w n − k 2 w n = R n :=R n r 2 , 1 < r < L, w n = 0, at r = 1, ∂ r w n − ikw n = g n − g * n =: h n , at r = L, whereR n is as in (65).

2.4.
Estimates of the correctors. We now look for uniform estimates of g * n and R n , which are independent of n and k. These results will be needed in Section 3.
Definition 2.4. A function or a constant depending on ε such asg ε ,g ε is called an exponentially small term, and denoted e.s.t., if there exist β, β > 0 such that for any q ≥ 0, there exists a constant c β,β ,q > 0 independent of ε such that Lemma 2.5. Let g * n be as in (63). Then, for |n| ≥ 2k , we find where C and c are constants independent of n, L, and k. Here 1 ≤ L 0 < L.
Lemma 2.6. LetR n be as in (65). Then, for m ≥ 0 and |n| ≥ 2k , we have r m 2R n 2 where C is independent of n and k.

3.
A priori estimates. In this section, we find a priori estimates for the solutions of the error equations (49) and (66) motivated from [46] and [45] (see also [11], [15], [19], [20] and [37]). The main purpose of this sections is to show that as k → ∞ in a proper space. We first consider the equations of w n where h n := g n − g * n as in (13) 3 , (25) and (63), and forR n as in (65). For the sake of simplicity, we drop the index n and we will write (·, ·) r s instead of (·, ·) r s ,(1,L) where s is an integer, see (17). We multiply (76) by rv with v ∈ H 1 r (1, L) and integrate over (1, L). We then derive the weak form of (76): To find w ∈ H 1 r (1, L) such that We first set v = w in (78), and look for the real and imaginary parts, separately. For the real part, we find and for the imaginary part, we find Then, (79) yields Here and after the δ j are proper positive constants which will be determined later on. From (80), we deduce that Combining (81) and (82), we obtain It remains to bound the term k 2 w 2 r in (83). We now set v = 2(r − 1)∂ r w in (78).
Hence, we obtain from (107) and (109) that We deduce from (102) where |n| ≥ 2k and C is independent of n and k. We finally arrive at the main theorem.
Theorem 3.1. Let u = u (k) s be the solution of (6) and Ω 0 = Ω * \ D where D is the unit disk centered at the origin and Ω * is a disk centered at (0, 0) with radius L 0 , 1 L 0 , which does not depend on k. Then, the following estimates hold: where Θ is the corrector given by where the ϕ n and θ n are as in (24) and (62), respectively.
Proof. We recall that Ω and B are as in (9) and (10), where L = C 1 k ≤ Ck − 7 4 + C L k 2 + where C is independent of n and k. Then, since w r,Ω0 ≤ w r,Ω where L 0 < L, the theorem follows.