The Dirichlet-to-Neumann map in a disk with a one-step radial potential. An analytical and numerical study

We consider the Schr\"odinger operator with a potential q on a disk and the map that associates to q the corresponding Dirichlet to Neumann (DtN) map. We give some numerical and analytical results on the range of this map and its stability, for the particular class of one-step radial potentials.


Introduction
Let Ω ⊂ R 2 be a bounded domain with smooth boundary ∂Ω. For each q ∈ L ∞ (Ω), consider the so called Dirichlet to Neumann map (DtN) given by where u is the solution of the following problem ∆u + q(x)u = 0, x ∈ Ω, u = f, ∂Ω, and ∂u ∂n | ∂Ω denotes the normal derivative of u on the boundary ∂Ω. Note that the uniqueness of u as solution of (2) requires that 0 is not a Dirichlet eigenvalue of ∆ + q. A sufficient condition to guarantee that Λ q is well-defined is to assume q(x) < λ 1 , the first Dirichlet eigenvalue of the Laplace operator in Ω, since in this case, the solution in (2) is unique. We assume that this condition holds and let us define the space L ∞ <λ 1 (Ω) = {q ∈ L ∞ (Ω), s. t. q(x) < λ 1 , a. e. }.
This has an important role in inverse problems where the aim is to recover the potential q from boundary measurements. In practice, these boundary measurements correspond to the associated DtN map and therefore, the mathematical statement of the classical inverse problem consists in the inversion of Λ. It is known that Λ is one to one as long as q ∈ L p with p > 2 (see [3], [Blȧsten, Imanuvilov and Yamamoto 2015]). Therefore, the inverse map Λ −1 can be defined in the range of Λ. There are, however, two related important and difficult questions that are not well understood: a characterization of the range of Λ and its stability i.e., a quantification of the difference of two potentials, in the L ∞ topology in terms of the distance of their associated DtN maps. Obviously, this stability will affect to the efficiency of any inversion or reconstruction algorithm to recover the potential from the DtN map (see [7], [Tejero 2016] and [8], [Tejero 2018]).
The first question, i.e. the characterization of the range of Λ is widely open. To our knowledge, the further result is due to Ingerman in [5] [Ingerman 2000], where a difficult characterization is obtained for the adherence with respect to a certain topology. Concerning the stability, there are some results when we assume that the potential q has some smoothness. In particular, if q ∈ H s (Ω) with s > 0, the following log −stability condition is known (see [3], [Blȧsten, Imanuvilov and Yamamoto 2015]), where V (t) = C log(1/t) −α for some constants C, α > 0. Stronger stability conditions are known in some particular cases. For example, in [2], [Beretta, De Hoop and Qiu 2013] it is shown that when q is piecewise constant and all the components where it takes a constant value touch the boundary, the stability is Lipschitz, i.e. there exists a constant C > 0 such that In this work we try to understand better the situation by considering the simplest case of a disk with one-step radial potentials q. More precisely, we give some results on the range of Λ and its stability when we restrict to the particular case Ω = B(0, 1) = {x ∈ R 2 : r = |x| < 1} and q ∈ F ⊂ L ∞ (Ω) given by (5) where χ (0,b) (r) is the characteristic function of the interval (0, b). Note that F is a two-parametric family depending on γ and b.
It is worth mentioning that, as we show below, the solution of (2) is unique for all b ∈ (0, 1) and γ ≥ 0, and therefore the DtN map is well defined for all these one step potentials. However, we restrict ourselves to the bounded set F to simplify.
Even in this simple case, a complete analytic answer to the previous questions is unknown. Therefore, we have considered a numerical approach based on a discrete sampling of the set F . Given an integer N > 0 we define h = 1/N and Note that F h has N (N − 1) + 1 functions from F . As h → 0 we obtain a better description of F and, in particular, we should recover the stability properties for q ∈ F . Concerning the stability of Λ, we show that it fails in the sense that inequality (4) does not hold for any continuous function V (t) with V (0) = 0. The proof relies on the ideas in [1] [Alessandrini 1988] where the analogous result is obtained for the conductivity problem.
We also obtain some partial stability results when b and γ are fixed. To state them we define the subsets F b ⊂ F , for b ∈ (0, 1), by and G γ ⊂ F , for γ ∈ [0, 1], by We prove that, for fixed b 0 > 0, Lipschitz stability holds, if we restrict ourselves to F b 0 . Therefore, the lack of stability is due to the change of b 0 rather than to changes in γ. Concerning G γ , we prove that if b ≥ b 0 > 0 there is stability of the DtN map with respect to b for potentials in G γ . This suggests a possible stability with respect to the L 1 -norm of the potentials which is sensitive to the position of the discontinuities, when consider discontinuous functions. In fact, we show numerical evidences of such stability when considering potentials in F .
For the range of Λ we give a characterization in terms of the first two eigenvalues of the DtN map. We also analyze the region where the stability constant is larger and, therefore, the potentials for which any recovering algorithm for q from the DtN map will have more difficulties.
It is worth mentioning that the results in this paper cannot be easily translated into the closely related, and more classical, conductivity problem where (2) is replaced by and the Dirichlet to Neumann map, or voltage to current map, is given by We refer to the review paper [9] [Uhlmann 2014] and the references therein for theoretical results on the DtN map in this case. The rest of this paper is divided as follows: In section 2 below we characterize the DtN map in terms of its eigenvalues using polar coordinates, in section 3 and 4 we analyze the stability and range results respectively. In section 5 we briefly describe the main conclusions and finally section 5 contains the proofs of the theorems stated in the previous sections.

The Dirichlet to Neumann map
In this section we characterize the Dirichlet to Neumann map in the case of a disk. System (2) in polar coordinates reads where v(r, θ) = u(r cos θ, r sin θ) and g(θ) = f (cos θ, sin θ) is a periodic function.
An orthonormal basis in L 2 (0, 2π) is given by {e inθ } n∈Z . Here we use this complex basis to simplify the notation but in the analysis below we only consider the subspace of real valued functions. Therefore, any function g ∈ L 2 (0, 2π) can be written as and g 2 L 2 (0,2π) = n∈Z |g n | 2 . Associated to this basis we define the usual Hilbert spaces H α # , for α > 0, as The Dirichlet to Neumann map in this case is defined as where v is the unique solution of (11).
In the above basis the Dirichlet to Neumann map turns out to be diagonal. In fact, we have the following result: Theorem 1 Let Ω be the unit disk and q ∈ F . Then, where and J n (r) are the Bessel functions of first kind.
For r ∈ (0, b) we solve the ODE with the boundary data at r = 0, while for r ∈ (b, 1) we use the boundary data at r = 1. In the first case the ODE is the Bessel ODE or orden 0 and therefore we have where J 0 is the Bessel function of the first kind and A 0 , C 0 are constants. These are computed by imposing continuity of a 0 and a 0 at r = b. In this way, we obtain b . Solving the system for A 0 and C 0 and taking into account that Λ q (1) = ∂v ∂r (1, θ) = a 0 (1) = C 0 we easily obtain (14).
Similarly, to compute c n in (14) we have to consider g(θ) = e inθ in (11) and therefore we assume that the solution v(r, θ) can be written in separate variables, i.e. v(r, θ) = a n (r)e inθ . Then, a n must satisfy r 2 a n + ra n + (r 2 q(r) − n 2 ) a n = 0, 0 < r < 1, a n (1) = 1, lim r−→0,r>0 a n (r) < ∞, n ≥ 1. (20) As in the case of c 0 , for r ∈ (0, b) we solve the ODE with the boundary data at r = 0, while for r ∈ (b, 1) we use the boundary data at r = 1. We have a n (r) = A n J n ( √ γr), r ∈ (0, b), C n (r n − r −n ) + r n , r ∈ (b, 1), where A n , C n are constants. These are computed by imposing continuity of a n and a n at r = b. In this way, we obtain Solving the system for A n and C n we obtain in particular We simplify this expression using the well known identity and we obtain, Now, taking into account that Λ q (e inθ ) = ∂v ∂r (1, θ) = a n (1)e inθ = (2nC n + n)e inθ we easily obtain (14).

Remark 2
In this proof of Theorem 1 we do not use the restriction γ ≤ 1 that satisfy the potentials in F . In fact, the statement of the theorem still holds for any step potential, as in F , but with any arbitrary large γ ≥ 0.

Stability
In this section we focus on the stability results for the map Λ. Some results are analytical and they are stated as theorems. The proofs are given in the appendix below. We divide this section in three subsections where we consider the negative stability result for q ∈ F norm and some partial results when we consider the subsets F b and G γ defined in (7) and (8).

Stability for q ∈ F
The first result in this section is the lack of any stability property when q ∈ F . In particular, we prove that inequality (4) fails, for any continuous function V (t) with V (0) = 0.
This result contradicts any possible stability result of the DtN map at q 0 ∈ F . Roughly speaking the idea is that the eigenvalues of Λ, given in Theorem 1 above, depend continuously on b, unlike the L ∞ norm of the potentials. A detailed proof of Theorem 3 is given in the appendix below.

Partial stability
We give now two partial stability results when we fix b and γ, respectively. Theorem 4 Given b ∈ (0, 1) and q 1 , q 2 ∈ F b , we have On the other hand, given γ ∈ (0, 1] and q 1 , q 2 ∈ G γ , we have where b = min{b 1 , b 2 }. The proof of this theorem is in the appendix below. Inequality (22) provides a Lipschitz stability result for Λ when b is fixed. This shows that the lack of Lipschitz stability is related to variations in the position of the discontinuity, which is the main idea in the negative result given in Theorem 3.
A numerical quantification of this Lipschitz stability for b fixed is easily obtained. We fix b = b 0 and consider and, for q 0 ∈ F h,b 0 , then C 2 (h, q 0 , b 0 ) remains bounded as h → 0 for all q 0 ∈ F h . In Figure 1 we show the behavior of C 2 (h, q 0 , b 0 ) when h = 10 −4 for different values of b 0 .
To illustrate the behavior with respect to b → 0 we plot in the left hand side of Figure 1 the graphs of the functions We see that both constants become larger for small values of b. We also see that both graphs are close in this logarithmic scale. However, the range of the interval [C 2,max (b), C 2,min (b)] is not small, as showed in the right hand side of Figure 1.
Concerning inequality (23) in Theorem 4, it provides a stability result of Λ with respect to the position of the discontinuity. In particular, this means that we can expect some Lipschitz stability if we consider a norm for the potentials that is sensitive to the position of the discontinuity. This is not the case for the L ∞ norm but it is true for the L p -norm for some 1 ≤ p < ∞.
In particular, is bounded as h → 0 and b ≥ b 0 > 0. In Figure 2 we show the values when h = 10 −4 . We observe that the constant blows up as b → 0.  As F is a bi-parametric family of potentials, it is natural to check if we can characterize the family {c n } n≥0 with only the first two coefficients c 0 and c 1 . In this section we give numerical evidences of the following facts:

Range of the DtN map
1. The first two coefficients c 0 and c 1 in (15)-(16) are the most sensitive with respect to (b, γ) and, therefore, the more relevant ones to identify b and γ from the DtN map.

The function
is injective. This means, in particular, that the DtN map can be characterized by the coefficients c 0 and c 1 , when restricted to functions in F h . We also illustrate the set of possible coefficients c 0 , c 1 .
3. The lack of stability for Λ is associated to higher density of points in the range of Λ h . This occurs when either b or γ are close to zero.

Sensitivity of c n
To analyze the relevance and sensitivity of the coefficients c n = c n (b, γ) to identify the parameters (b, γ) we have computed their range when (b, γ) ∈ [0, 1] × [0, 1], and the norm of their gradients. As we see in (4.1) the range decreases for large n. This means that, for larger values of n, the variability of c n is smaller and they are likely to be less relevant to identify q. However, even if the range of c n becomes smaller for large n they could be more sensitive to small perturbations in (b, γ) and this would make them useful to distinguish different potentials. But this is not the case. In Figure  3 we show that for the given values of γ = 0.1, 0.34, 0.67, 0.99 and b ∈ (0, 1] the gradients of the first two coefficients, with respect to (b, γ), are larger than the others. Therefore, we conclude that the two first coefficients c 0 and c 1 are the most sensitive, and therefore relevant, to identify the potential q.
We also see in Figure 3 that these gradients are very small for b << 1. This means, in particular, that identifying potentials with small b from the DtN map should be more difficult.  Table 1: Range of the coefficients, i.e. for each c n the range is defined as max q∈F h c n − min q∈F h c n . Figure 3: Norm of the gradient of the coefficients c n (γ, b) in terms of b ∈ (0, 1) for different values of γ. We see that the gradients of higher coefficients n ≥ 2 are smaller than those of the first two ones. We also observe that these gradients become small for small values of b. Coordinates lines for fixed γ and b are given in Figure 5. We observe that R(Λ h ) is a convex set between the curves r low : {(c 0 (γ, 1), c 1 (γ, 1)), with γ ∈ [0, 1]}, Note also that, in the c 0 , c 1 plane, the length of the coordinates lines associated to b constant are segments that become smaller as b → 0. Analogously, the length of those associated to constant γ become smaller as γ → 0. Thus, the region where either b or γ are small produces the higher density of points in the range of Λ h . This corresponds to the upper left part of its range (see Figure 4). On the other hand, this Figure provides a numerical evidence of the injectivity of Λ h too. In fact, any point inside R(Λ h ) is the intersection of two coordinates lines associated to some unique b 0 and γ 0 . The higher density of points in the upper left hand side of the range of Λ h should correspond to potentials q with large stability constant C 2 (h, q). In Figure 6 we show the level sets of C 2 (h, q) for h = 10 −4 and different q ∈ F h . The region with larger constant corresponds to small values of b (upper right figure) and larger values of c 1 (upper left and lower figures). On the other hand, the region with lower stability constant is for b close to b = 1, which corresponds to the lower part of the range of Λ h when c 0 is small. Figure 6: Level sets of the C 2 (b, γ) for q ∈ F h and h = 10 −4 in terms of (b, γ) (upper left) and in terms of (c 0 , c 1 ) (upper right). A zoom of the upper left region in this last figure is in the lower figure. Regions separated by the level sets are indicated: region I corresponds to the potentials with stability constant larger that 10 7 , region II corresponds to those with stability constant lower that 10 7 but larger than 10 6 , and so on.
It is interesting to analyze the set of potentials with the same coefficient c 0 or c 1 . We give in Figure 7 the coordinates lines of the inverse map (Λ h ) −1 . When increasing the value of either c 0 (light lines) or c 1 (dark lines) we obtain lines closer to the left part of the (b, γ) region. We see that the angle between coordinate lines becomes very small for b small. In this region, close points could be the intersection of coordinates lines associated to not so close parameters (b, γ). This agrees with the region where the stability constant is larger.

Conclusions
We have considered the relation between the potential in the Schrödinger equation and the associated DtN map in one of the simplest situations, i.e. for the subset of radial one step potentials in dimension 2. In particular we have focused on two difficult problems: the stability of the map Λ (defined in 3) and its range. In this case, the map Λ is easily characterized in terms of the Bessel functions and this allows us to give some analytical and numerical results on these problems. We have proved the lack of any possible stability result, by adapting the argument in [1] [Alessandrini 1988] for the conductivity problem. We have also obtained some partial Lipschitz stability when the position of the discontinuity is fixed in the potential and numerical evidences of the stability with respect to the L 1 norm. Finally, we have characterized numerically the range of Λ in terms of the first two eigenvalues of the DtN map and given some insight in the regions where stability of Λ is worse.
The following lemma will be used in the proof of Theorem 4.
Proof. The previous estimates are consequence of (33) and the inequality cos r − cos s = 2 sin s + r 2 sin s − r 2 ≤ s 2 − r 2 2 .