INVERSE SOURCE PROBLEMS WITHOUT (PSEUDO) CONVEXITY ASSUMPTIONS

. We study the inverse source problem for the Helmholtz equation from boundary Cauchy data with multiple wave numbers. The main goal of this paper is to study the uniqueness and increasing stability when the (pseudo)convexity or non-trapping conditions for the related hyperbolic prob- lem are not satisﬁed. We consider general elliptic equations of the second order and arbitrary observation sites. To show the uniqueness we use the analytic continuation, the Fourier transform with respect to the wave numbers and uniqueness in the lateral Cauchy problem for hyperbolic equations. Numerical examples in 2 spatial dimension support the analysis and indicate the increasing stability for large intervals of the wave numbers, while analytic proofs of the increasing stability are not available.

1. Introduction and main results. This paper studies the inverse source problem for the Helmholtz equation from boundary Cauchy data with multiple frequencies (wave numbers). The source term is supported in a bounded domain Ω and shall be recovered from full/partial boundary data on Γ ⊂ ∂Ω. Such inverse source problems have wide applications in antenna synthesis [2], biomedical imaging [1], and various kinds of tomography. It has been demonstrated that by boundary data of one single frequency, we could not identify the source function uniquely [11,Ch.4] while a family of equations (like the Helmholtz equation for various wave numbers in (0, K)) regains uniqueness. The crucial issue then is a stability estimate which is of importance in applications. Generally a stability estimate in inverse problems for elliptic equations is logarithmic which yields a robust recovery of only few source parameters and provides low resolution numerically. We mention that in [7] some uniqueness and numerical evidence were provided highlighting the improvement of multiple frequency measurement data. In [3], [4] the authors have proven the uniqueness and the first increasing stability results were obtained in more particular case by a different method of spatial Fourier analysis. In [6] we observed that the inverse source problem for Helmholtz equations, in odd dimensions, is highly connected with the identification of the initial data in the hyperbolic initial value problem by lateral Cauchy data (boundary observability in control theory) by means of the Fourier transform in time. As a consequence we also obtained most complete increasing (with growing K) stability results. In a recent paper [15] we further traced the exponential dependence on the attenuation term. We shall emphasize that in all papers mentioned above increasing stability bounds assume that the corresponding lateral Cauchy problem satisfies the so called non-trapping condition which is guaranteed by certain (pseudo)convexity of the domain and the surface with Cauchy data with respect to the hyperbolic equation.
Other better investigated uniqueness and increasing stability by higher wave numbers measurement focus on the Cauchy problem of elliptic equations and the Schrödinger potential identification. Classical Carleman estimates imply some conditional Hölder type stability estimates for solutions of the elliptic Cauchy problem. Nevertheless, F. John in [16] showed that, in the continuation problem for the Helmholtz equation from the unit disk onto any larger disk, the stability estimate is still of a logarithmic type, uniformly with respect to the wave numbers. Increasing stability for the Schrödinger potential and conductivity coefficient has been demonstrated in [13], [14] respectively. In particular, in [13] we also traced the dependence of increasing stability on the attenuation. It has been demonstrated in [9] that the stability always is improving for larger k under (pseudo) convexity conditions on the geometry of the domain and on the coefficients of the elliptic equation. On the other hand, [12] and [10] verify that in some cases increasing stability holds even when convexity and non-trapping conditions fail.
In this paper we rigorously demonstrate uniqueness in the inverse source problems without any convexity/non-trapping condition and give a strong numerical evidence of increasing stability for source identification. In absence of such conditions, increasing stability estimates for sources identification are not known.
Let Ω be a bounded Lipschitz domain in R n , n = 2, 3. We introduce another sub-domain Ω 0 such thatΩ 0 ⊂ Ω, whose boundary ∂Ω 0 is of C 2 class, and R n \Ω 0 is connected.
Let A be the second order elliptic operator and A := ∆ in R n \ Ω. Here we denote b = (b 1 , ..., b n ). We assume the ellipticity condition 0 |ξ| 2 ≤ n j,m=1 a jm (x)ξ j ξ m for some positive constant 0 and all x, ξ ∈ R n as well as the additional conditions the imaginary part Imb and (2) 0 ≤ Imc in Ω.
We consider the scattering problem where f 1 ∈ L 2 (R n ), f 0 ∈ H 1 (R n ), are real valued, and suppf 0 ∪ suppf 1 ⊂ Ω \ Ω 0 with the Neumann boundary condition where a is the symmetric matrix (a jm ), and the Sommerfeld radiation condition lim r=|x|→+∞ r n−1 We will assume that for some positive numbers C, δ depending on G 1 . In other words, G 1 is a function of (x, t) whose indicated norms with respect to x ∈ ∂Ω 0 decay exponentially with respect to t. Then g 1 (·, k) is (complex) analytic with respect to The inverse source problem of interest is to find f 0 , f 1 , g 1 entering (3)-(5) from the additional Cauchy data where Γ is a non(empty) C 2 hyper surface in ∂Ω, see Figure 1, and K 0 , K are some positive numbers satisfying K 0 < K. We present the main analytic result below: Theorem 1.1. Let the assumptions (1), (2) hold. Then a solution (f 0 , f 1 , g 1 ) to the inverse source problems (3)-(5) with the additional Cauchy data (7) is unique in each of the cases a): This result will be proven in the next section by using analyticity of u(·, k) with respect to k and an auxiliary hyperbolic initial boundary value problem obtained after the Fourier-Laplace transform of the scattering problem with respect to k.
The source terms f 0 , f 1 in the scattering problem will be the initial values for this hyperbolic problem and G 1 will be the boundary value data on ∂Ω 0 × (0, +∞). Using the John-Tataru's uniqueness of the continuation theorem one shows that the solution of the hyperbolic equation and hence the initial and boundary data are uniquely determined.
Observe that if Γ = ∂Ω, then in Theorem 1.1 one can replace the Cauchy data (7) by the Dirichlet data u(·, k) = u 0 on Γ, K 0 < k < K. Indeed, the radiating solution of the exterior Dirichlet problem is unique, so u 0 uniquely determines u outside Ω and hence u 1 = ∂ ν u is attainable on ∂Ω. Moreover, if ∂Ω is analytic and connected, and A = ∆, b 0 = 0 near ∂Ω, then in Theorem 1.1 instead of (7) it is sufficient to use u(·, k) = u 0 on Γ ⊂ ∂Ω. Indeed, then u(·, k) is (real) analytic near ∂Ω, so u on Γ by the analytic continuation uniquely determines u on ∂Ω, and as above, also ∂ ν u on Γ.
At any fixed K one expects a logarithmic stability in this inverse source problem which is quite discouraging for applications. At present, there is no analytic results showing increasing stability for larger K without convexity assumptions. In section 3 we present numerical examples of reconstruction of sources in the plane which either are characteristic functions of relatively complicated domains or smoothly spread point sources described by around 100 parameters. This reconstruction dramatically improves with growing K definitely suggesting improving stability. While we consider the Helmholtz equation with constant coefficients the shape and location of the observation curve Γ violates a non-trapping condition. Indeed, supports of sources are outside of the convex hull of Γ.
2. Proof of uniqueness. We first show solvability of the direct scattering problem and analyticity of its solution with respect to the wave number k by the Lax-Phillips method [11].
Let B be a ball containingΩ. We will use the following elliptic boundary value problem Observe that the boundary value problem (8) is elliptic and Fredholm. To show its unique solvability in the Sobolev space H 2 (Ω 1 ), it suffices to demonstrate that the homogeneous problem (f * = 0, g 1 = 0) has only the zero solution. Such argument follows from the energy integral. Indeed, the standard definition of a weak solution gives Integrating by parts, we obtain due to the conditions (1), (2) on the coefficients b.
because of our assumptions on the coefficients b 0 , c. Using the boundary conditions in (8) and taking the imaginary parts in (9) we conclude that The above inequalities (10), (11) show that the boundary integral ∂B VV is zero. Therefore V = 0 on ∂B and from the boundary condition, ∂ ν V = 0 on ∂B.
By uniqueness in the Cauchy problem for the elliptic equation (8) we have that V = 0 in Ω 1 . According to the theory of elliptic boundary value problems the operator mapping (8) has the continuous inverse. Since this operator is continuous and analytic in these spaces, by known results ( [18, pp.365 We fix a smaller ball B 1 whose closure is in B and containsΩ. Let χ be a cut off C ∞ -function with χ = 1 on Ω and χ = 0 outside B 1 . When showing solvability and analyticity of the scattering problem (3)-(5) we can assume that g 1 = 0. Indeed, u = u s + χV (0, g 1 ; k), where u s solves the scattering problem (3)-(5) with g 1 = 0 and f is replaced by which is a complex analytic function of k ∈ S with the values in H 0 (B). So it suffices to show existence and analyticity of u s . For brevity, from now on we denote u s by u.
We look for a solution solving the elliptic boundary value problem (8) and with the radiating fundamental solution Φ to the Helmholtz equation. Here f * is a function to be determined later. We extend f, f * onto Ω 0 as zero. We recall that . Defining E(f * ; k) as zero in Ω we conclude that u solves our scattering problem if and only if To show the unique solvability of (13) and analytic dependence of its solution on k ∈ S it suffices to prove uniqueness of its solution. The operator E(·; k) is analytic and compact from L 2 (B) into itself. Indeed, the well known elliptic estimates show that both w and V are continuous linear operators from Since the E(f * ) is a linear first order partial differential operator with respect to w − V , it is continuous from L 2 (B) into H 1 (B) and hence compact from L 2 (B) into L 2 (B). As a consequence, the equation (13) is Fredholm, so existence of its solution and the existence of the continuous inverse of I + E follow from the uniqueness. Moreover, since the right hand side in (13) is an analytic continuous operator (from L 2 (B) into itself) with respect to k ∈ S, again from the known results about analytic continuous operators in Banach spaces ([18, pp.365]) it follows that its (continuous) inversion exists and is analytic with respect to k ∈ S 1 (an open subset of S containing To show the uniqueness, we let f = 0 in (13) and we will prove that f * = 0. From (13) we have that f * = 0 in Ω and in B \ B 1 . Then u solves the homogeneous scattering problem (3), (4), (5), and, as we show next, u = 0 in B \ Ω 0 . Therefore, In what follows we denote B(R) = {x : |x| < R}. To demonstrate uniqueness in the direct homogeneous scattering problem (3), (4), (5) as in (9), (10), (11) we have (14) Im Im(c+ikb 0 +k 2 )uū = 0.
In the case 1) using (10), (11), we have and using the radiation condition (5) from the known results we derive that u = 0 outside B and hence by the uniqueness of the continuation for the elliptic equations in B \ Ω 0 .
In the case 2) we will use the integral representation where u * ∈ H 2 (B) is an extension of the function u from R n \ B 1 onto B. Indeed, both u * and the right hand side of (15) solve the same Helmholtz equation in R n and satisfy the Sommerfeld radiation condition, so they are equal and we have (15).
Since the fundamental solution (12) decays exponentially in |x| when 0 < k 2 (which is obvious when n = 3 and follows from known inequality |H and hence letting R → +∞ in (14) and using (10), (11) we conclude that Now we prove Theorem 1.1.
Proof. Due to the linearity it suffices to show that u 0 = u 1 = 0 on Γ, K 0 < k < K implies that f 0 = f 1 = 0, g 1 = 0. Let U be a solution to the hyperbolic mixed initial boundary value problem As known in [17] there is a unique solution U ∈ L ∞ ((0, T ); H 1 (R n \Ω 0 )) of this problem for any T > 0 and moreover by using the standard energy estimates, i.e. multiplying (17) by ∂ t U e −γ0t and integrating by parts over (R n \ Ω 0 ) × (0, t) we have for some positive γ 0 = γ 0 (U ). Then the following Fourier-Laplace transform is well defined Approximating f 0 , f 1 , g 1 by smooth functions and integrating by parts we obtain Hence u * solves (3)-(4) with k = k 1 + iγ, γ 0 < γ. In addition, u * exponentially decays for large |x|.
3. Numerical tests. In this section, we provide some examples numerically verifying the improving stability, when the wave number grows, without convexity assumptions. For simplicity's sake we choose f 0 = 0 in (3) and only consider the recovery of a real source function f 1 . We let Ω 0 to be an empty set. The numerical tests for the inverse source problems are separated into two different cases. In the first case, we consider an interior inverse source problem such that suppf 1 ⊂ Ω and the observation site Γ ⊂ ∂Ω. Whereas in the second case, we discuss an exterior inverse source problem where the compact support of f 1 is contained in the exterior domain of the convex hull of the observation site Γ which is not covered in the above discussion.
3.1. Case one: Interior inverse source problems. We firstly consider the interior inverse source problem with suppf 1 ⊂ Ω which fits well the setting in Theorem 1.1.
where k 0 is the wave number. We assume that f 1 ∈ L 2 (R 2 ) and is compactly supported in Ω. The Cauchy data is gathered at the observation site Γ = (0.3 + 2 cos θ, 0.3 + 2 sin θ) with θ ∈ (0, π) which is half of the boundary of the sphere B 2 (0.3, 0.3). For simplicity's sake, we define the following two forward operators is the fundamental solution of (22).
In order to approximate the unknown source in Ω, we discretize a subset of the domain Ω radically and spherically. More precisely, we discretize the radius equally in [0.5, 1.5] with 30 points and the angle [0, 2π] with 100 points. Mesh-points for the source function are generated by the products of radical discretization and the angle discretization with 3000 unknowns. The observation data is gathered at the Γ with equally distributed 50 angles in (0, π). The collocation method is considered to generate the discretized matrices of L 1,k and L 2,k . The same Kaczmarz-Landweber algorithm in [4,5,6,15] is implemented to recover the unknown source. The wave number k for the Kaczmarz iteration is chosen in the set K = {1, 2, . . . , 200} which is equally distributed in [1,200] and the inner Landweber iteration is fixed by 500 steps.
As shown in Figure 2 we consider a piecewise constant source function as follows We present the reconstructed results in Figure 3. The left panel presents the approximate source function by the exact measurement data for the largest wave number k = 200. The middle panel shows the error between the exact source and the approximate one. In the right panel, we collect the relative error versus the increasing wave number k. As one can observe, the relative error decreases substantially when the wave number increases. The final iterates in the right panel has a relative error of 0.1856 where the piecewise constant structure is well captured. On the other hand, we shall emphasize that because the wave number interval [1,200] might not be sufficient in this inverse source problem and the error between both sources is still comparably large.  The approximate source function, error between both sources and the relative error are presented in Figure 4. Similar summary can be concluded referring to the piecewise constant source above.  . Annular domain with different observation sites for exact measurement data. Γ = (0.3 + 2 cos θ, 0.3 + 2 sin θ), θ ∈ (0, π); Γ 1 = (0.3 + 5 cos θ, 0.3 + 5 sin θ), θ ∈ (0, 2 5 π) and Γ 2 = (0.3 + 10 cos θ, 0.3+10 sin θ), θ ∈ (0, 1 5 π). Left: the relative error slope versus the increasing wave number for the piecewise constant source; Right: the relative error slope versus the increasing wave number for the mixed-point source. The values at the tail of each line are the minimal relative errors with k = 200 or k = 400.
The proposed Kaczmarz-Landweber iteration scheme is quite robust on noisy data. If the noise is small, i.e. less than 5%, the approximate source functions mimic both exact sources in good manners as those of [4,6,15]. For simplicity's sake we skip all the discussion on the noise data but focus on different numerical tests if the observation site varies. More precisely, instead of the site Γ = (0.3 + 2 cos θ, 0.3+2 sin θ), θ ∈ (0, π) in above examples, we also consider another two sites Γ 1 = (0.3 + 5 cos θ, 0.3 + 5 sin θ), θ ∈ (0, 2 5 π) and Γ 2 = (0.3 + 10 cos θ, 0.3 + 10 sin θ), θ ∈ (0, 1 5 π) where all the traces Γ, Γ 1 , Γ 2 have the same length and the source domain is fixed as B 1.5 (0, 0) \B 0.5 (0, 0). For both the piecewise constant source (23) and the mixed-point source (24), the relative error slopes with increasing wave number are collected in Figure 5. Both panels show that the distance between the support of f 1 and the observation site is essential for the performance of the Kaczmarz-Landweber iteration scheme. When the observation site have certain distance to the compactly supported domain, the multi-frequency measurement data for wave numbers k ∈ [1,200] can not improve the resolution essentially. To obtain better accuracy one has to increase the wave number interval as shown in Figure 5 where k ∈ [1, 400] is considered for large distance cases of Γ 1 and Γ 2 . For exact data, we present the approximate source, error with respect to the exact source and the relative error slope in Figure 7. It can be observed that with a short observation site, the approximate source does not mimic the exact one accurately. One can either extend the wave number set or the length of Γ to improve the resolution.
At the same time, we focus on the influence of the distance between suppf 1 and the observation site Γ. Different from the annular domain discussed above, we    We collect the numerical results for exact measurements in Figure 10. With shorter length, the approximate source has a comparably larger relative error with respect to the exact one. Referring to the right panel of Figure 10, the relative error for the final wave number k = 200 is 0.3923.  Figure 11.
For exact data, we present the approximate source and error between both sources and the relative error slope in Figure 12. Notice that we have four times more of the observation site length compared with that of Figure 6. With longer observation site Γ, the relative error for k = 200 is much smaller compared with those in Subsection 3.1.2.

4.
Conclusion. Uniqueness can be extended onto isotropic Maxwell and linear elasticity systems by using uniqueness of the continuation results in [8] instead of [19]. The needed scattering theory is available, although not so transparent and explicit as for the Helmholtz equation.  The next really challenging issue is to obtain increasing stability bounds, similar to those in [6], [15]. At present, there is no decisive idea how to approach this question.
It is important to collect further numerical evidence of the increasing stability for more complicated geometries not satisfying non-trapping conditions.
An increasing (to optimal Lipschitz one) stability without convexity conditions has a fundamental significance, since it suggests that an arbitrary location of a measurement site Γ can still provide with a high resolution in many important applications.
One expects increasing stability in the inverse source problem when scattering is replaced by stationary waves in a bounded domain, i.e. when one replaces R n by a bounded reference domain Ω 1 and the radiation condition by one of standard boundary conditions on ∂Ω 1 (for example, by the Neumann condition which is most important in applications). In these cases there are additional difficulties due to eigenvalues of elliptic boundary value problems in bounded domains.
One can look at the different inverse source problem which is the linearized inverse problem for the Schrödinger potential: find f (supported in Ω ⊂ R 3 ) from Ω f (y) e ki|x−y| |x − y| e ki|z−y| |z − y| dy given for x, z ∈ Γ ⊂ ∂Ω, where one expects quite explicit increasing stability bounds for f for larger k.