Explicit tight bounds on the stably recoverable information for the inverse source problem

For the inverse source problem with the two-dimensional Helmholtz equation, the singular values of the 'source-to-near field' forward operator reveal a sharp frequency cut-off in the stably recoverable information on the source. We prove and numerically validate an explicit, tight lower bound for the spectral location of this cut-off. We also conjecture and support numerically a tight upper bound for the cut-off. The bounds are expressed in terms of zeros of Bessel functions of the first and second kind.


Introduction
We treat the single-frequency inverse source problem for the Helmholtz equation in the plane, illustrated in figure 1. Fix a positive constant wavenumber k = 2π/λ, where λ is the operating frequency, and let D 0 and D be open disks in R 2 centered at the origin and with radii R 0 and R ≥ R 0 , respectively. Write ∆ = ∂ 2 x 1 + ∂ 2 x 2 for the Laplacian, and consider the Helmholtz problem (∆ + k 2 )u = s in R 2 , lim |x|→∞ |x|(∂ |x| − ik)u(x) = 0, uniformly for x/|x| ∈ S 1 , for some source s ∈ L 2 (D 0 ) extended by zero to the whole plane. The second condition in (1) is the outgoing Sommerfeld radiation condition in the plane. The inverse source problem, ISP, is now given a single measurement U ∈ L 2 (∂D), find a source s ∈ L 2 (D 0 ) such that there is a function ('radiated field') u satisfying u| ∂D = U and satisfying the system (1).
The ISP arises naturally in inverse acoustic and electromagnetic scattering, and has been devoted a substantial body of literature. The ISP is treated, e.g., in the multifrequency regime by Bao, Lin and Triki (2010), and with far-field measurement data by Griesmaier, Hanke and Raasch (2012); see also El Badia and Nara (2011). It occurs in antenna synthesis and diagnostics (Persson and Gustafsson, 2005;Jørgensen et al., 2010), the analytic continuation of solutions of exterior scattering problems (Sternin and Shatalov, 1994;Zaridze, 1998;Bliznyuk, Pogorzelski and Cable, 2005;Karamehmedović, 2015), and in linearized inverse obstacle scattering problems.
In terms of the forward operator F : s → U , described in detail in section 2, solving the ISP amounts to solving F s = U for s ∈ L 2 (D 0 ).
This problem is ill-posed, since ker Figure 2. A schematic of the singular value spectrum σ m and σ −1 m of the forward operator F and its pseudoinverse F † , respectively, as function of angular frequency m of right singular vectors of F . The lower bandwidth bound B − is given in Theorem 1. The upper bandwidth bound B + is predicted in Conjecture 1. Both B − and B + are validated numerically in section 3.
noisy and sampled over a finite set of points. A common regularizing measure is to look for the minimum-L 2 -norm, or minimum-energy, solution of (2), which is given by Another regularization scheme uses a truncated singular value decomposition (TSVD) of the forward operator F . Here, s is approximated by a finite sum of the form N 1 σ −1 m (U, φ) m ψ m , with (σ m , ψ m , φ m ) a singular system of F . Our aim is to estimate the maximal amount of information about any source s ∈ L 2 (D 0 ) that can be stably recovered in principle, that is, regardless of the sampling frequency in the measurement and of the choice of the regularisation scheme. By 'stably recoverable information' we mean 'information recoverable robustly to noise,' and we refer to figure 2 for a more precise definition. A non-asymptotic analysis of the singular values σ m of the forward operator F , performed in section 2, reveals a low-pass filter behavior with well-defined passband and stopband. This turns out to be true also when the singular values are ordered according to increasing angular frequency of the right singular vectors of F , that is, of the singular vectors defined at the measurement boundary ∂D. In this case, the singular values within the passband generally do not increase or decrease monotonically, and the singular values in the stopband, still ordered according to angular frequency m, are monotonic functions of m. We call the bandwidth B of the forward operator F the singular value index (angular frequency m of a right singular vector of F ) at which the singular value spectrum of F becomes strictly decreasing as function of nonnegative m: With this in mind, we define the stably recoverable information on a source s to be the projection of s onto the singular subspace of F defined by |m| ≤ B. Then, finding the maximal amount of stably recoverable information about any source s, regardless of measurement sampling quality and of regularization scheme, amounts to estimating the bandwidth B of the forward operator F .
To simplify the notation, write κ 0 = kR 0 and κ = kR for the size parameters of the source support and of the measurement boundary, respectively. Also, for integer m, write j m,1 and y m,1 for the first positive zero of the Bessel function J m of the first kind, respectively Bessel function Y m of the second kind, and order m. It is wellknown (Magnus, Oberhettinger and Soni, 1966, p. 146) that j m,1 > 0 for all m ∈ N 0 . Our main result, proved in Section 2.2, is Theorem 1. The bandwidth B of the forward operator F : s → U associated with the Helmholtz problem (1) and measurement at ∂D is bounded from below by For convenience, in Section 2.2 we also show that the bandwidth bound of Theorem 1 can be expressed explicitly in the source size parameter κ 0 : Corollary 1. For sufficiently large κ 0 , we have Finally, the general form of the result in Theorem 1, as well as extensive numerical experimentation, lead us to conjecture a tight upper bound on the bandwitdth B: In section 2 we analyze the singular value spectrum of the forward operator F . In particular, we prove Theorem 1 and Corollary 1 in section 2.2. We validate the bounds B − and B + on the bandwidth B numerically in section 3, and discuss some implications of Theorem 1 in section 4. A conclusion and suggestions for further work are given in section 5.

Spectral analysis of the forward operator
The function (i/4)H (1) 0 (k|x|), x ∈ R 2 , is the radial outgoing fundamental solution of the Helmholtz operator in the plane, with singularity at the origin. Recall that H (1) 0 = J 0 + iY 0 is the Hankel function of zero order and of the first kind. As in Bao, Lin and Triki (2010), introduce the forward operator that maps sources s to the traces at ∂D of the corresponding radiated fields. It is wellknown (Bao, Lin and Triki, 2010) where H (2) 0 = J 0 − iY 0 is the Hankel function of zero order and of the second kind.

A singular system of F
Bao, Lin and Triki (2010) derived a singular system of the forward operator F . We here slightly improve a part of their Proposition 2.1: Lemma 1. The forward operator F admits the singular value decomposition and Our slight improvement of Proposition 2.1 of Bao, Lin and Triki (2010) consists in explicitly evaluating the integral R 0 =0 J 2 m (k ), occurring in σ m and ψ m , in terms of A m (κ 0 ). This explicit evaluation is crucial to our proof of Theorem 1. We also note that our expressions for the singular vectors φ m , as well as the singular values σ m , differ from Bao, Lin and Triki (2010) in that they are only proportional to those given in that reference.
Proof of Lemma 1. For s ∈ L 2 (D 0 ) and y ∈ D 0 we have A special case of the Graf addition theorem (Abramowitz and Stegun, 1972, Eq. 9.1.79, p. 363) reads Similar to Bao, Lin and Triki (2010), inserting this in (4) we get This gives an eigendecomposition of the operator F * F ; to normalize the eigenvectors, we note that Gradsteyn and Ryzhik (2007, Eq. 5.54.2, p. 629) gives and the recursion formula for cylinder functions (Gradsteyn and Ryzhik, 2007, Eq. 8.471.1, p. 926) implies Thus, and F * F admits the spectral decomposition Evidently, σ 2 0 has multiplicity one and all the other eigenvalues σ 2 m , m ∈ N, have multiplicity two. The lemma now follows from Theorem 4.7 on p. 100 of Colton and Kress (2013); it here just remains to compute Figure 3 shows the first 71 nonnegative-index singular values of the forward operator F with size parameters κ = κ 0 = 10π. Clearly, the forward operator is a low-pass filter with respect to the singular values σ m , with bandwidth B = 27. We quantify the frequency response of this filter in section 2.2.

Proof of Theorem 1 and of Corollary 1
In this section we prove the lower bound B − on the bandwidth B given in Theorem 1, and the approximate value of B − given in Corollary 1. For completeness, we first prove that the distance between the zeros of the function 0 ≤ µ → J µ (κ 0 ) is greater than 1.
We can now link the variation of the function m → A m (κ 0 ) with that of the Bessel function of the first kind. Fix m ∈ N 0 . Lemma 3. If J ξ (κ 0 ) = 0 for some ξ ∈ [m, m + 1] then A m (κ 0 ) ≤ A m+1 (κ 0 ).

Proof. The recursion formula (5) implies
For any fixed positive argument x, the function R µ → J µ (x) is differentiable and not identically zero. Thus, by assumption, and by lemma 2, this function changes sign precisely once in the interval [m, m + 1], so Remark 1. Clearly, the function N 0 m → |H (1) m (κ)| 2 is positive-valued. It is also strictly increasing, as can be seen from Nicholson's integral for |H (1) m (κ)| 2 (Watson, 1945, pp. 441-444), where K 0 is the modified Bessel function of the second kind. This namely implies since both K 0 and the hyperbolic sine are positive over positive reals.
The above discussion suffices for a proof of the lower bound B − .

Numerical validation
We here compute the bandwidth B, as well as the bandwidth bounds B − and B + of Theorem 1 and Conjecture 1, respectively, for 300 values of the size parameters κ = κ 0 uniformly distributed over the interval κ ∈ [2, 100π]. Recall that κ = kR = 2πR/λ and κ 0 = kR 0 = 2πR 0 /λ, where R is the radius of the sampling circle ∂D, R 0 is the radius of the source domain, and λ is the operating wavelength. Thus, we consider 300 values of the relative wavelength λ/R = λ/R 0 distributed nonuniformly over the interval λ/R ∈ [1/50, π]. Figure 5 shows the errors ε ± = B ± − B and the relative errors ε rel,± = |B ± − B|/B in the estimated bandwidth as function of the problem size parameter κ.
For the two lowest considered values of κ, we find that B = 0 and B − = 0; there, we set ε rel,− = 0. Both B and B − are positive for higher considered values of κ. In particular, there is zero bandwidth for κ smaller than some threshold value between approx. 1.7 and approx. 2.7, and for such size parameters κ the inverse source problem is, from the viewpoint of the bandwidth of the singular values, similar to the inverse heat conduction problem. Over the considered interval for κ, the mean errors are ε − = −1.68, ε + = 3.02, and the maximum absolute errors are max |ε − | = 3, max |ε + | = 4. The relative error in B − is below 5% for κ ≥ 24.7461, i.e., for R/λ ≥ 3.94, and B + is below 5% for κ ≥ 45.7181, i.e., for λ/R ≥ 7.28. We find both B, B − and B + to be approximately linear functions of κ in the given interval, with least-squares fits summarized in Table 3. Figure 6 shows errors in the approximations B ± (for definition of B ± see proof of Corollary 1 and the paragraph immediately following it, on page 9.) The approximate expression for the lower bound shows almost the same small error as the lower bound itself, and the approximate expression B + ≈ κ 0 for the upper bound, while simple, has error below 5% only for problem size parameters of approx. 175 or higher when κ = κ 0 is maintained.
Our bounds B ± are independent of the radius of the measurement surface, and we next validate this property numerically. Figure 7 shows the first 71 nonnegative-index

Discussion
The bandwidth estimates B ± are directly applicable as optimal filter estimates in the numerical solution of the inverse source problem in terms of a truncated singular value decomposition (TSVD) of the forward operator. Next, it has been amply observed in the literature concerning the single-frequency inverse source problem that the numerical stability of the solution increases with the operating frequency. Theorem 1 confirms and explicitly quantifies this increase in numerical stability, also for non-asymptotic frequencies.
Theorem 1 of course has direct implications for the maximum achievable stable resolution of the reconstruction in the inverse source problem. Detailed analysis of this resolution requires an investigation of the pointwise behavior of the left singular vectors of the forward operator. While we here do not perform such analysis, we do note that the left singular vectors tend to be supported near the origin for low values of index m, and near the measurement boundary for high index values. This means the amplified noise produces is a 'wall of non-information' near the measurement boundary and blocks faithful reconstruction of the source inside D.
As shown in Bao, Lin and Triki (2010) and in Section 2 here, the right singular vectors (defined over the measurement boundary) of the forward operator are proportional to exp(imθ), m ∈ Z. This means that the bandwidth index B is approximately the angular frequency of the highest-frequency data component that can be stably inverted. Thus, the sampling theorem (Shannon, 1949) is directly applicable with Theorem 1 to give the following: in case the radiated field u is sampled equidistantly at the boundary ∂D, any angular sampling rate greater than approximately ∆θ ≈ π/B ≤ π/B − = π/argmin m∈N 0 {j m,1 ≥ κ 0 } is excessive due to the limited bandwidth of the forward operator. the inversion to noise generally improving as the source support approaches the measurement boundary.

Conclusion and further work
We analyzed the singular values of the forward operator associated with the singlefrequency inverse source problem for the Helmholtz equation in the plane. In particular, we considered bounds on the information content that is preserved by the forward operator, proving a tight lower bound and conjecturing a tight upper bound on the singular value index of the highest-frequency data component that is stably recoverable. The bounds were expressed in terms of the zeros of Bessel functions of the first and the second kind. We validated both bounds numerically, establishing concrete estimates on the stably recoverable information in the inverse source problem regardless of the data sampling rate and the choice of regularization. The result can be used directly, e.g., to estimate optimal TSVD filters and data sampling rates.
Proving the statement in Conjecture 1 is a natural next step. Also, it would complete the picture to supplement the results on the bandwidth with a more precise description of the general levels and decay rates of the singular values as function of the size parameters of the source support and of the measurement boundary, individually or in relation to one another. Finally, a spectral analysis of the forward operator in dimension greater than 2 will be interesting.