Phenomena

for a sound speed c(x). This equation can be considered as a model for other hyperbolic equations, and the methods presented here can in some cases be extended to study wave phenomena in other fields such as electromagnetism or elasticity.

We will mostly be interested in inverse problems for the wave equation. In these problems, one has access to certain measurements of waves (the solutions u) on the surface of a medium, and one would like to determine material parameters (the sound speed c) of the interior of the medium from these boundary measurements. A typical field where such problems arise is seismic imaging, where one wishes to determine the interior structure of Earth by making various measurements of waves at the surface. We will not describe seismic imaging applications in more detail here, since they are discussed elsewhere in this volume.

Another feature in this chapter is that we will consistently consider anisotropic materials, where the sound speed depends on the direction of propagation. This means that the scalar sound speed c(x), where $x = (x^{1},x^{2},\ldots,x^{n}) \in \Omega \subset \mathbb{R}^{n}$ , is replaced by a positive definite symmetric matrix $(g^{jk}(x))_{j,k=1}^{n}$ , and the wave equation becomes

$\displaystyle{ \partial _{t}^{2}u(x,t) -\sum _{ j,k=1}^{n}g^{jk}(x) \frac{\partial ^{2}u} {\partial x^{j}\partial x^{k}}(x,t) = 0. }$

Anisotropic materials appear frequently in applications such as in seismic imaging.

It will be convenient to interpret the anisotropic sound speed (g ^jk) as the inverse of a Riemannian metric, thus modeling the medium as a Riemannian manifold . The benefits of such an approach are twofold. First, the well-established methods of Riemannian geometry become available to study the problems, and second, this provides an efficient way of dealing with the invariance under changes of coordinates present in many anisotropic wave imaging problems. The second point means that in inverse problems in anisotropic media, one can often only expect to recover the matrix (g ^jk) up to a change of coordinates given by some diffeomorphism. In practice, this ambiguity could be removed by some a priori knowledge of the medium properties (such as the medium being in fact isotropic, see section “From Boundary Distance Functions to Riemannian Metric”).

2 Background

This chapter contains three parts which discuss different topics related to wave imaging. The first part considers the inverse problem of determining a sound speed in a wave equation from the response operator, also known as the hyperbolic Dirichlet-to-Neumann map, by using the boundary control method; see [5, 7, 42]. The second part considers other types of boundary measurements of waves, namely, the scattering relation and boundary distance function, and discusses corresponding inverse problems. The third part is somewhat different in nature and does not consider any inverse problems but rather gives an introduction to the use of curvelet decompositions in wave imaging for nonsmooth sound speeds. We briefly describe these three topics.

Wave Imaging and Boundary Control Method

Let us consider an isotropic wave equation. Let $\Omega \subset \mathbb{R}^{n}$ be an open, bounded set with smooth boundary $\partial \Omega$ , and let c(x) be a scalar-valued positive function in $C^{\infty }(\overline{\Omega })$ modeling the wave speed in $\Omega$ . First, we consider the wave equation

$\displaystyle\begin{array}{rcl} & \partial _{t}^{2}u(x,t) - c(x)^{2}\Delta u(x,t) = 0\quad \mbox{ in }\Omega \times \mathbb{R}_{+},& \\ & u\vert _{t=0} = 0,\quad u_{t}\vert _{t=0} = 0, & \\ & c(x)^{-n+1}\partial _{\mathbf{n}}u = f(x,t)\quad \mbox{ in }\partial \Omega \times \mathbb{R}_{+}, &{}\end{array}$

(1)

where $\partial _{\mathbf{n}}$ denotes the Euclidean normal derivative and n is the unit interior normal. We denote by $u^{f} = u^{f}(x,t)$ the solution of (1) corresponding to the boundary source term f.

Let us assume that the domain $\Omega \subset \mathbb{R}^{n}$ is known. The inverse problem is to reconstruct the wave speed c(x) when we are given the set

$\displaystyle\begin{array}{rcl} \left \{\left (f\vert _{\partial \Omega \times (0,2T)},\,u^{f}\vert _{ \partial \Omega \times (0,2T)}\right ):\ f \in C_{0}^{\infty }(\partial \Omega \times \mathbb{R}_{ +})\right \},& & {}\\ \end{array}$

that is, the Cauchy data of solutions corresponding to all possible boundary sources $f \in C_{0}^{\infty }(\partial \Omega \times \mathbb{R}_{+})$ , T ∈ (0, ∞]. If T = ∞, then this data is equivalent to the response operator

$\displaystyle\begin{array}{rcl} \Lambda _{\Omega }: f\mapsto u^{f}\vert _{ \partial \Omega \times \mathbb{R}_{+}},& &{}\end{array}$

(2)

which is also called the nonstationary Neumann-to-Dirichlet map . Physically, $\Lambda _{\Omega }f$ describes the measurement of the medium response to any applied boundary source f, and it is equivalent to various physical measurements. For instance, measuring how much energy is needed to force the boundary value $c(x)^{-n+1}\partial _{\mathbf{n}}u\vert _{\partial \Omega \times \mathbb{R}_{+}}$ to be equal to any given boundary value $f \in C_{0}^{\infty }(\partial \Omega \times \mathbb{R}_{+})$ is equivalent to measuring the map $\Lambda _{\Omega }$ on $\partial \Omega \times \mathbb{R}_{+}$ ; see [42, 44]. Measuring $\Lambda _{\Omega }$ is also equivalent to measuring the corresponding Neumann-to-Dirichlet map for the heat or the Schrödinger equations or measuring the eigenvalues and the boundary values of the normalized eigenfunctions of the elliptic operator $-c(x)^{2}\Delta$ ; see [44].

The inverse problems for the wave equation and the equivalent inverse problems for the heat or the Schrödinger equations go back to works of M. Krein at the end of the 1950s, who used the causality principle in dealing with the one-dimensional inverse problem for an inhomogeneous string, $u_{tt} - c^{2}(x)u_{xx} = 0$ ; see, for example, [46]. In his works, causality was transformed into analyticity of the Fourier transform of the solution. A more straightforward hyperbolic version of the method was suggested by A. Blagovestchenskii at the end of 1960s to 1970s [12, 13]. The multidimensional case was studied by M. Belishev [4] in the late 1980s who understood the role of the PDE control for these problems and developed the boundary control method for hyperbolic inverse problems in domains of Euclidean space. Of crucial importance for the boundary control method was the result of D. Tataru in 1995 [77, 79] concerning a Holmgren-type uniqueness theorem for nonanalytic coefficients. The boundary control method was extended to the anisotropic case by M. Belishev and Y. Kurylev [7]. The geometric version of the boundary control method which we consider in this chapter was developed in [7, 41, 42, 47]. We will consider the inverse problem in the more general setting of an anisotropic wave equation in an unbounded domain or on a non-compact manifold. These problems have been studied in detail in [39, 43] also in the case when the measurements are done only on a part of the boundary. In this paper we present a simplified construction method applicable for non-compact manifolds in the case when measurements are done on the whole boundary. We demonstrate these results in the case when we have an isotropic wave speed c(x) in a bounded domain of Euclidean space. For this, we use the fact that in the Euclidean space, the only conformal deformation of a compact domain fixing the boundary is the identity map. This implies that after the abstract manifold structure (M, g) corresponding to the wave speed c(x) in a given domain $\Omega$ is constructed, we can construct in an explicit way the embedding of the manifold M to the domain $\Omega$ and determine c(x) at each point $x \in \Omega$ . We note on the history of this result that using Tataru’s unique continuation result [77], Theorem 2 concerning this case can be proven directly using the boundary control method developed for domains in Euclidean space in [4].

The reconstruction of non-compact manifolds has been considered also in [11, 27] with different kind of data, using iterated time reversal for solutions of the wave equation. We note that the boundary control method can be generalized also for Maxwell and Dirac equations under appropriate geometric conditions [50, 51], and its stability has been analyzed in [1, 45].

Travel Times and Scattering Relation

The problem considered in the previous section of recovering a sound speed from the response operator is highly overdetermined in dimensions n ≥ 2. The Schwartz kernel of the response operator depends on 2n variables, and the sound speed c depends on n variables.

In section “Travel Times and Scattering Relation,” we will show that other types of boundary measurements in wave imaging can be directly obtained from the response operator. One such measurement is the boundary distance function , a function of 2n − 2 variables, which measures the travel times of shortest geodesics between boundary points. The problem of determining a sound speed from the travel times of shortest geodesics is the inverse kinematic problem . The more general problem of determining a Riemannian metric (corresponding to an anisotropic sound speed) up to isometry from the boundary distance function is the boundary rigidity problem . The problem is formally determined if n = 2 but overdetermined for n ≥ 3.

This problem arose in geophysics in an attempt to determine the inner structure of the Earth by measuring the travel times of seismic waves. It goes back to Herglotz [37] and Wiechert and Zoeppritz [84] who considered the case of a radial metric conformal to the Euclidean metric. Although the emphasis has been in the case that the medium is isotropic, the anisotropic case has been of interest in geophysics since the Earth is anisotropic. It has been found that even the inner core of the Earth exhibits anisotropic behavior [24].

To give a proper definition of the boundary distance function, we will consider a bounded domain $\Omega \subset \mathbb{R}^{n}$ with smooth boundary to be equipped with a Riemannian metric g, that is, a family of positive definite symmetric matrices $g(x) =\big (g_{jk}(x)\big)_{j,k=1}^{n}$ depending smoothly on $x \in \overline{\Omega }$ . The length of a smooth curve $\gamma: [a,b] \rightarrow \overline{\Omega }$ is defined to be

$\displaystyle{ L_{g}(\gamma ) =\int _{ a}^{b}\left (\sum _{ j,k=1}^{n}g_{ jk}\big(\gamma (t)\big)\dot{\gamma }^{j}(t)\dot{\gamma }^{k}(t)\right )^{1/2}\,dt. }$

The distance function d _g(x, y) for $x,y \in \overline{\Omega }$ is the infimum of the lengths of all piecewise smooth curves in $\overline{\Omega }$ joining x and y. The boundary distance function is d _g(x, y) for $x,y \in \partial \Omega$ .

In the boundary rigidity problem, one would like to determine a Riemannian metric g from the boundary distance function d _g. In fact, since $d_{g} = d_{\psi ^{{\ast}}g}$ for any diffeomorphism $\psi: \overline{\Omega } \rightarrow \overline{\Omega }$ which fixes each boundary point, we are looking to recover from d _gthe metric g up to such a diffeomorphism. Here, $\psi ^{{\ast}}g(y) = D\psi (y)^{t}g\big(\psi (y)\big)D\psi (y)$ is the pullback of g by ψ.

It is easy to give counterexamples showing that this cannot be done in general; consider, for instance, the closed hemisphere, where boundary distances are given by boundary arcs so making the metric larger in the interior does not change d _g. Michel [55] conjectured that a simple metric g is uniquely determined, up to an action of a diffeomorphism fixing the boundary, by the boundary distance function d _g(x, y) known for all x and y on $\partial \Omega$ . A metric is called simple if for any two points in $\overline{\Omega }$ , there is a unique length minimizing geodesic joining them, and if the boundary is strictly convex.

The conjecture of Michel has been proved for two-dimensional simple manifolds [60]. In higher dimensions, it is open, but several partial results are known, including the recent results of Burago and Ivanov for metrics close to Euclidean [15] and close to hyperbolic [16] (see the survey [40]). Earlier and related works include results for simple metrics conformal to each other [8, 10, 26, 56–58], for flat metrics [34], for locally symmetric spaces of negative curvature [9], for two-dimensional simple metrics with negative curvature [25, 59], a local result [70], a semiglobal solvability result [54], and a result for generic simple metrics [71].

In case the metric is not simple, instead of the boundary distance function, one can consider the more general scattering relation which encodes, for any geodesic starting and ending at the boundary, the start point and direction, the end point and direction, and the length of the geodesic. We will see in section “Travel Times and Scattering Relation” that also this information can be determined directly from the response operator. If the metric is simple, then the scattering relation and boundary distance function are equivalent, and either one is determined by the other.

The lens rigidity problem is to determine a metric up to isometry from the scattering relation. There are counterexamples of manifolds which are trapping, and the conjecture is that on a nontrapping manifold the metric is determined by the scattering relation up to isometry. We refer to [72] and the references therein for known results on this problem.

Curvelets and Wave Equations

In section “Curvelets and Wave Equations,” we describe an alternative approach to the analysis of solutions of wave equations, based on a decomposition of functions into basic elements called curvelets or wave packets . This approach also works for wave speeds of limited smoothness unlike some of the approaches presented earlier. Furthermore, the curvelet decomposition yields efficient representations of functions containing sharp wave fronts along curves or surfaces, thus providing a common framework for representing such data and analyzing wave phenomena and imaging operators. Curvelets and related methods have been proposed as computational tools for wave imaging, and the numerical aspects of the theory are a subject of ongoing research.

A curvelet decomposition was introduced by Smith [67] to construct a solution operator for the wave equation with C ^1, 1 sound speed and to prove Strichartz estimates for such equations. This started a body of research on L ^pestimates for low-regularity wave equations based on curvelet-type methods; see, for instance, Tataru [80–82], Smith [68], and Smith and Sogge [69]. Curvelet decompositions have their roots in harmonic analysis and the theory of Fourier integral operators, where relevant works include Córdoba and Fefferman [23] and Seeger et al. [65] (see also Stein [73]).

In a rather different direction, curvelet decompositions came up in image analysis as an optimally sparse way of representing images with C ² edges; see Candés and Donoho [20] (the name “curvelet” was introduced in [19]). The property that curvelets yield sparse representations for wave propagators was studied in Candés and Demanet [17, 18]. Numerical aspects of curvelet-type methods in wave computation are discussed in [21, 30]. Finally, both theoretical and practical aspects of curvelet methods related to certain seismic imaging applications are studied in [2, 14, 29, 31, 64].

3 Mathematical Modeling and Analysis

Boundary Control Method

Inverse Problems on Riemannian Manifolds

Let $\Omega \subset \mathbb{R}^{n}$ be an open, bounded set with smooth boundary $\partial \Omega$ and let c(x) be a scalar-valued positive function in $C^{\infty }(\overline{\Omega })$ , modeling the wave speed in $\Omega$ . We consider the closure $\overline{\Omega }$ as a differentiable manifold M with a smooth, nonempty boundary. We consider also a more general case, and allow (M, g) to be a possibly non-compact, complete manifold with boundary. This means that the manifold contains its boundary $\partial M$ and M is complete with metric d _gdefined below. Moreover, near each point x ∈ M, there are coordinates (U, X), where U ⊂ M is a neighborhood of x and $X: U \rightarrow \mathbb{R}^{n}$ if x is an interior point, or $X: U \rightarrow \mathbb{R}^{n-1} \times [0,\infty )$ if x is a boundary point such that for any coordinate neighborhoods (U, X) and $(\tilde{U},\tilde{X})$ , the transition functions $X \circ \tilde{ X}^{-1}:\tilde{ X}(U \cap \tilde{ U}) \rightarrow X(U \cap \tilde{ U})$ are C ^∞-smooth. Note that all compact Riemannian manifolds are complete according to this definition. Usually we denote the components of X by $X(y) =\big (x^{1}(y),\ldots,x^{n}(y)\big)$ .

Let u be the solution of the wave equation

$\displaystyle\begin{array}{rcl} u_{tt}(x,t) + Au(x,t)& =& 0\quad \mbox{ in }\quad M \times \mathbb{R}_{+}, \\ u\vert _{t=0}& =& 0,\quad u_{t}\vert _{t=0} = 0, \\ B_{\nu,\eta }u\vert _{\partial M\times \mathbb{R}_{+}}& =& f. {}\end{array}$

(3)

Here, $f \in C_{0}^{\infty }(\partial M \times \mathbb{R}_{+})$ is a real-valued function, and A = A(x, D) is an elliptic partial differential operator of the form

$\displaystyle\begin{array}{rcl} Av = -\sum _{j,k=1}^{n}\mu (x)^{-1}\vert g(x)\vert ^{-\frac{1} {2} } \frac{\partial } {\partial x^{j}}\left (\mu (x)\vert g(x)\vert ^{\frac{1} {2} }g^{jk}(x) \frac{\partial v} {\partial x^{k}}(x)\right ) + q(x)v(x),& &{}\end{array}$

(4)

where $g^{jk}(x)$ is a smooth, symmetric, real, positive definite matrix, $\vert g\vert = \mbox{ det}\big(g^{jk}(x)\big)^{-1}$ , and μ(x) > 0 and q(x) are smooth real-valued functions. On existence and properties of the solutions of Eq. (3), see [52]. The inverse of the matrix $\big(g^{jk}(x)\big)_{j,k=1}^{n}$ , denoted $\big(g_{jk}(x)\big)_{j,k=1}^{n}$ defines a Riemannian metric on M. The tangent space of M at x is denoted by T _xM, and it consists of vectors p which in local coordinates (U, X), $X(y) =\big (x^{1}(y),\ldots,x^{n}(y)\big)$ are written as $p =\sum _{ k=1}^{n}p^{k} \frac{\partial } {\partial x^{k}}$ . Similarly, the cotangent space T _x^∗ M of M at x consists of covectors which are written in the local coordinates as $\xi =\sum _{ k=1}^{n}\xi _{k}dx^{k}$ . The inner product which g determines in the cotangent space T _x^∗ M of M at the point x is denoted by $\langle \xi,\eta \rangle _{g} = g(\xi,\eta ) =\sum _{ j,k=1}^{n}g^{jk}(x)\xi _{j}\eta _{k}$ for $\xi,\eta \in T_{x}^{{\ast}}M$ . We use the same notation for the inner product at the tangent space T _xM, that is, $\langle p,q\rangle _{g} = g(p,q) =\sum _{ j,k=1}^{n}g_{jk}(x)p^{j}q^{k}$ for p, q ∈ T _xM.

The metric defines a distance function, which we call also the travel time function,

$\displaystyle\begin{array}{rcl} d_{g}(x,y) =\inf \vert \mu \vert,\quad \vert \mu \vert =\int _{ 0}^{1}\langle \partial _{ s}\mu (s),\partial _{s}\mu (s)\rangle _{g}^{1/2}ds,& & {}\\ \end{array}$

where | μ | denotes the length of the path μ, and the infimum is taken over all piecewise C ¹-smooth paths μ: [0, 1] → M with μ(0) = x and μ(1) = y.

We define the space L ²(M, dV _μ) with inner product

$\displaystyle\begin{array}{rcl} \langle u,v\rangle _{L^{2}(M,dV _{\mu })} =\int _{M}u(x)v(x)\,dV _{\mu }(x),& & {}\\ \end{array}$

where $dV _{\mu } =\mu (x)\vert g(x)\vert ^{1/2}dx^{1}dx^{2}\ldots dx^{n}$ . By the above assumptions, A is formally self-adjoint, that is,

$\displaystyle\begin{array}{rcl} \langle Au,v\rangle _{L^{2}(M,dV _{\mu })} =\langle u,Av\rangle _{L^{2}(M,dV _{\mu })}\quad \mbox{ for }u,v \in C_{0}^{\infty }(M^{\mathrm{int}}).& & {}\\ \end{array}$

Furthermore, let

$\displaystyle\begin{array}{rcl} B_{\nu,\eta }v = -\partial _{\nu }v +\eta v,& & {}\\ \end{array}$

where $\eta: \partial M \rightarrow \mathbb{R}$ is a smooth function and

$\displaystyle\begin{array}{rcl} \partial _{\nu }v =\sum _{ j,k=1}^{n}\mu (x)g^{jk}(x)\nu _{ k} \frac{\partial } {\partial x^{j}}v(x),& & {}\\ \end{array}$

where $\nu (x) = (\nu _{1},\nu _{2},\ldots,\nu _{m})$ is the interior conormal vector field of $\partial M$ , satisfying $\sum _{j,k=1}^{n}g^{jk}\nu _{j}\xi _{k} = 0$ for all cotangent vectors of the boundary, $\xi \in T^{{\ast}}(\partial M)$ . We assume that ν is normalized, so that $\sum _{j,k=1}^{n}g^{jk}\nu _{j}\nu _{k} = 1$ . If M is compact, then the operator A in the domain $\mathcal{D}(A) =\{ v \in H^{2}(M):\ \partial _{\nu }v\vert _{\partial M} = 0\}$ , where H ^s(M) denotes the Sobolev spaces on M, is an unbounded self-adjoint operator in $L^{2}(M,dV _{\mu })$ .

An important example is the operator

$\displaystyle\begin{array}{rcl} A_{0} = -c^{2}(x)\Delta + q(x)& &{}\end{array}$

(5)

on a bounded smooth domain $\Omega \subset \mathbb{R}^{n}$ with $\partial _{\nu }v = c(x)^{-n+1}\partial _{\mathbf{n}}v$ , where $\partial _{\mathbf{n}}v$ is the Euclidean normal derivative of v.

We denote the solutions of (3) by

$\displaystyle\begin{array}{rcl} u(x,t) = u^{f}(x,t).& & {}\\ \end{array}$

For the initial boundary value problem (3), we define the nonstationary Robin-to-Dirichlet map or the response operator $\Lambda$ by

$\displaystyle{ \Lambda f = u^{f}\vert _{ \partial M\times \mathbb{R}_{+}}. }$

(6)

The finite time response operator $\Lambda ^{T}$ corresponding to the finite observation time T > 0 is given by

$\displaystyle{ \Lambda ^{T}f = u^{f}\vert _{ \partial M\times (0,T)}. }$

(7)

For any set $B \subset \partial M \times \mathbb{R}_{+}$ , we denote $L^{2}(B) =\{ f \in L^{2}(\partial M \times \mathbb{R}_{+}):\ \mathop{\mathrm{supp}}(f) \subset B\}$ . This means that we identify the functions and their zero continuations.

By [78], the map $\Lambda ^{T}$ can be extended to bounded linear map $\Lambda ^{T}: L^{2}(B) \rightarrow H^{1/3}\big(\partial M \times (0,T)\big)$ when $B \subset \partial M \times (0,T)$ is compact. Here, $H^{s}\big(\partial M \times (0,T)\big)$ denotes the Sobolev space on $\partial M \times (0,T)$ . Below we consider $\Lambda ^{T}$ also as a linear operator $\Lambda ^{T}: L_{cpt}^{2}\big(\partial M \times (0,T)\big) \rightarrow L^{2}\big(\partial M \times (0,T)\big)$ , where $L_{cpt}^{2}\big(\partial M \times (0,T)\big)$ denotes the compactly supported functions in $L^{2}\big(\partial M \times (0,T)\big)$ .

For t > 0 and a relatively compact open set $\Gamma \subset \partial M$ , let

$\displaystyle\begin{array}{rcl} M(\Gamma,t) =\{ x \in M\:\ d_{g}(x,\Gamma ) < t\}.& &{}\end{array}$

(8)

This set is called the domain of influence of $\Gamma$ at time t.

When $\Gamma \subset \partial M$ is an open relatively compact set and $f \in C_{0}^{\infty }(\Gamma \times \mathbb{R}_{+})$ , it follows from finite speed of wave propagation (see, e.g., [38]) that the wave $u^{f}(t) = u^{f}(\cdot,t)$ is supported in the domain $M(\Gamma,t)$ , that is,

$\displaystyle\begin{array}{rcl} u^{f}(t) \in L^{2}\big(M(\Gamma,t)\big) =\{ v \in L^{2}(M):\ \mathop{\mathrm{supp}}(v) \subset M(\Gamma,t)\}.& &{}\end{array}$

(9)

We will consider the boundary of the manifold $\partial M$ with the metric $g_{\partial M} =\iota ^{{\ast}}g$ inherited from the embedding $\iota: \partial M \rightarrow M$ . We assume that we are given the boundary data, that is, the collection

$\displaystyle\begin{array}{rcl} (\partial M,g_{\partial M})\mbox{ and }\Lambda,& &{}\end{array}$

(10)

where $(\partial M,g_{\partial M})$ is considered as a smooth Riemannian manifold with a known differentiable and metric structure and $\Lambda$ is the nonstationary Robin-to-Dirichlet map given in (6).

Our goal is to reconstruct the isometry type of the Riemannian manifold (M, g), that is, a Riemannian manifold which is isometric to the manifold (M, g). This is often stated by saying that we reconstruct (M, g) up to an isometry. Our next goal is to prove the following result:

Theorem 1.

Let (M,g) to be a smooth, complete Riemannian manifold with a nonempty boundary. Assume that we are given the boundary data (10). Then it is possible to determine the isometry type of manifold (M,g).

From Boundary Distance Functions to Riemannian Metric

In order to reconstruct (M, g), we use a special representation, the boundary distance representation, R(M), of M and later show that the boundary data (10) determine R(M). We consider next the (possibly unbounded) continuous functions $h: C(\partial M) \rightarrow \mathbb{R}$ . Let us choose a specific point $Q_{0} \in \partial M$ and a constant C ₀ > 0 and using these, endow $C(\partial M)$ with the metric

$\displaystyle\begin{array}{rcl} d_{C}(h_{1},h_{2}) = \vert h_{1}(Q_{0}) - h_{2}(Q_{0})\vert +\sup _{z\in \partial M}\min (C_{0},\vert h_{1}(z) - h_{2}(z)\vert ).& &{}\end{array}$

(11)

Consider a map $R: M \rightarrow C(\partial M)$ ,

$\displaystyle\begin{array}{rcl} R(x) = r_{x}(\cdot );\quad r_{x}(z) = d_{g}(x,z),\,\,z \in \partial M,& &{}\end{array}$

(12)

that is, $r_{x}(\cdot )$ is the distance function from x ∈ M to the points on $\partial M$ . The image $R(M) \subset C(\partial M)$ of R is called the boundary distance representation of M. The set R(M) is a metric space with the distance inherited from $C(\partial M)$ which we denote by d _C, too. The map R, due to the triangular inequality, is Lipschitz,

$\displaystyle\begin{array}{rcl} d_{C}(r_{x},r_{y}) \leq 2d_{g}(x,y).& &{}\end{array}$

(13)

We note that when M is compact and C ₀ = diam (M), the metric $d_{C}: C(\partial M) \rightarrow \mathbb{R}$ is a norm which is equivalent to the standard norm $\|f\|_{\infty } =\max _{x\in \partial M}\vert f(x)\vert$ of $C(\partial M)$ .

We will see below that the map $R: M \rightarrow R(M) \subset C(\partial M)$ is an embedding. Many results of differential geometry, such as Whitney or Nash embedding theorems, concern the question how an abstract manifold can be embedded to some simple space such as a higher dimensional Euclidean space. In the inverse problem, we need to construct a “copy” of the unknown manifold in some known space, and as we assume that the boundary is given, we do this by embedding the manifold M to the known, although infinite dimensional function space $C(\partial M)$ .

Next we recall some basic definitions on Riemannian manifolds; see, for example, [22] for an extensive treatment. A path μ: [a, b] → N is called a geodesic if, for any c ∈ [a, b], there is $$$\varepsilon > 0$$” src=”/wp-content/uploads/2016/04/A183156_2_En_52_Chapter_IEq105.gif”> such that if s, t ∈ [a, b] such that <IMG alt=$ , the path μ([s, t]) is a shortest path between its endpoints, that is,

$\displaystyle{\vert \mu ([s,t])\vert = d_{g}\big(\mu (s),\mu (t)\big).}$

In the future, we will denote a geodesic path μ by γ and parameterize γ with its arclength s, so that $\vert \mu ([s_{1},s_{2}])\vert = d_{g}\big(\mu (s_{1}),\mu (s_{2})\big)$ . Let x(s),

$\displaystyle{x(s) = \left (x^{1}(s),\ldots,x^{n}(s)\right ),}$

be the representation of the geodesic γ in local coordinates (U, X). In the interior of the manifold, that is, for $U \subset M^{\mathrm{int}}$ the path x(s) satisfies the second-order differential equations

$\displaystyle\begin{array}{rcl} \frac{d^{2}x^{k}(s)} {ds^{2}} = -\sum _{i,j=1}^{n}\Gamma _{ ij}^{k}\big(x(s)\big)\frac{dx^{i}(s)} {ds} \frac{dx^{j}(s)} {ds},& &{}\end{array}$

(14)

where $\Gamma _{ij}^{k}$ are the Christoffel symbols, given in local coordinates by the formula

$\displaystyle{\Gamma _{ij}^{k}(x) =\sum _{ p=1}^{n}\frac{1} {2}g^{kp}(x)\left (\frac{\partial g_{jp}} {\partial x^{i}} (x) + \frac{\partial g_{ip}} {\partial x^{j}} (x) -\frac{\partial g_{ij}} {\partial x^{p}} (x)\right ).}$

Let y ∈ M and $\xi \in T_{x}M$ be a unit vector satisfying the condition $$$g\big(\xi,\nu (y)\big) > 0$$” src=”/wp-content/uploads/2016/04/A183156_2_En_52_Chapter_IEq111.gif”> in the case when <IMG alt=$ . Then, we can consider the solution of the initial value problem for the differential equation (14) with the initial data

$\displaystyle\begin{array}{rcl} x(0) = y,\quad \frac{dx} {ds} (0) =\xi.& & {}\\ \end{array}$

This initial value problem has a unique solution x(s) on an interval $[0,s_{0}(y,\xi ))$ such that $$$s_{0}(y,\xi ) > 0$$” src=”/wp-content/uploads/2016/04/A183156_2_En_52_Chapter_IEq114.gif”> is the smallest value s 0 > 0 for which <IMG alt=$ , or $s_{0}(y,\xi ) = \infty$ in case no such s ₀ exists. We will denote $x(s) =\gamma _{y,\xi }(s)$ and say that the geodesic is a normal geodesic starting at y if $y \in \partial M$ and $\xi =\nu (y)$ .

Example 1.

In the case when (M, g) is such a compact manifold that all geodesics are the shortest curves between their endpoints and all geodesics can be continued to geodesics that hit the boundary, we can see that the metric spaces (M, d _g) and $(R(M),\,\|\,\cdot \|_{\infty })$ are isometric. Indeed, for any two points x, y ∈ M, there is a geodesic γ from x to a boundary point z, which is a continuation of the geodesic from x to y. As in the considered case the geodesics are distance minimizing curves, we see that

$\displaystyle\begin{array}{rcl} r_{x}(z) - r_{y}(z) = d_{g}(x,z) - d_{g}(y,z) = d_{g}(x,y),& & {}\\ \end{array}$

and thus $\|r_{x} - r_{y}\|_{\infty }\geq d_{g}(x,y)$ . Combining this with the triangular inequality, we see that $\|r_{x} - r_{y}\|_{\infty } = d_{g}(x,y)$ for x, y ∈ M and R is isometry of (M, d _g) and $(R(M),\,\|\,\cdot \|_{\infty })$ .

Notice that when even M is a compact manifold, the metric spaces (M, d _g) and $(R(M),\|\,\cdot \|_{\infty })$ are not always isometric. As an example, consider a unit sphere in $\mathbb{R}^{3}$ with a small circular hole near the South pole of, say, diameter $\varepsilon$ . Then, for any x, y on the equator and $z \in \partial M$ , $\pi /2-\varepsilon \leq r_{x}(z) \leq \pi /2$ and $\pi /2-\varepsilon \leq r_{y}(z) \leq \pi /2$ . Then $d_{C}(r_{x},r_{y}) \leq \varepsilon$ , while d _g(x, y) may be equal to π.

Next, we introduce the boundary normal coordinates on M. For a normal geodesic $\gamma _{z,\nu }(s)$ starting from $z \in \partial M$ consider $d_{g}(\gamma _{z,\nu }(s),\partial M)$ . For small s,

$\displaystyle\begin{array}{rcl} d_{g}(\gamma _{z,\nu }(s),\partial M) = s,& &{}\end{array}$

(15)

and z is the unique nearest point to γ _z, ν(s) on $\partial M$ . Let τ(z) ∈ (0, ∞] be the largest value for which (15) is valid for all s ∈ [0, τ(z)]. Then for s > τ(z),

$\displaystyle{d_{g}(\gamma _{z,\nu }(s),\partial M) < s,}$

and z is no more the nearest boundary point for γ _z, ν(s). The function $\tau (z) \in C(\partial M)$ is called the cut locus distance function, and the set

$\displaystyle\begin{array}{rcl} \omega =\{\gamma _{z,\nu }\big(\tau (z)\big) \in M:\, z \in \partial M,\ \mbox{ and }\tau (z) < \infty \},& &{}\end{array}$

(16)

is the cut locus of M with respect to $\partial M$ . The set ω is a closed subset of M having zero measure. In particular, M∖ω is dense in M. In the remaining domain M∖ω, we can use the coordinates

$\displaystyle\begin{array}{rcl} x\mapsto \big(z(x),t(x)\big),& &{}\end{array}$

(17)

where $z(x) \in \partial M$ is the unique nearest point to x and $t(x) = d_{g}(x,\partial M)$ . (Strictly speaking, one also has to use some local coordinates of the boundary, $y: z\mapsto \big(y^{1}(z),\ldots,y^{(n-1)}(z)\big)$ and define that

$\displaystyle\begin{array}{rcl} x\mapsto \big(y\big(z(x)\big),t(x)\big) = \left (y^{1}(z(x)\big),\ldots,y^{(n-1)}\big(z(x)\big),t(x)\right ) \in \mathbb{R}^{n},& &{}\end{array}$

(18)

are the boundary normal coordinates.) Using these coordinates, we show that $R: M \rightarrow C(\partial M)$ is an embedding. The result of Lemma 1 is considered in detail for compact manifolds in [42].

Lemma 1.

Let (M,d _g ) be the metric space corresponding to a complete Riemannian manifold (M,g) with a nonempty boundary. The map R: (M,d _g ) → (R(M),d _C ) is a homeomorphism. Moreover, given R(M) as a subset of $C(\partial M)$ , it is possible to construct a distance function d _R on R(M) that makes the metric space (R(M),d _R ) isometric to (M,d _g ).

Proof.

We start by proving that R is a homeomorphism. Recall the following simple result from topology:

Assume that X and Y are Hausdorff spaces, X is compact, and F: X → Y is a continuous, bijective map from X to Y. Then F: X → Y is a homeomorphism.

Let us next extend this principle. Assume that (X, d _X) and (Y, d _Y) are metric spaces and let X _j ⊂ X, $j \in \mathbb{Z}_{+}$ be compact sets such that $\bigcup _{j\in \mathbb{Z}_{+}}X_{j} = X$ . Assume that F: X → Y is a continuous, bijective map. Moreover, let Y _j = F(X _j) and assume that there is a point p ∈ Y such that

$\displaystyle\begin{array}{rcl} a_{j} =\inf _{y\in Y \setminus Y _{j}}d_{Y }(y,p) \rightarrow \infty \mbox{ as }j \rightarrow \infty.& &{}\end{array}$

(19)

Then by the above, the maps $F: \cup _{j=1}^{n}X_{j} \rightarrow \cup _{j=1}^{n}Y _{j}$ are homeomorphisms for all $n \in \mathbb{Z}_{+}$ . Next, consider a sequence y _k ∈ Y such that y _k → y in Y as k → ∞. By removing first elements of the sequence $(y_{k})_{k=1}^{\infty }$ if needed, we can assume that d _Y(y _k, y) ≤ 1. Let now $N \in \mathbb{Z}_{+}$ be such that for j > N, we have $$$a_{j} > b:= d_{Y }(y,p) + 1$$” src=”/wp-content/uploads/2016/04/A183156_2_En_52_Chapter_IEq148.gif”>. Then <IMG alt=$ and as the map $F:\bigcup _{ j=1}^{N}X_{j} \rightarrow \bigcup _{j=1}^{N}Y _{j}$ is a homeomorphism, we see that $F^{-1}(y_{k}) \rightarrow F^{-1}(y)$ in X as k → ∞. This shows that F ⁻¹: Y → X is continuous, and thus F: X → Y is a homeomorphism.

By definition, $R: M \rightarrow R(M)$ is surjective and, by (13), continuous. In order to prove the injectivity, assume the contrary, that is, $r_{x}(\cdot ) = r_{y}(\cdot )$ but x ≠ y. Denote by z ₀ any point where

$\displaystyle\begin{array}{rcl} \min _{z\in \partial M}r_{x}(z) = r_{x}(z_{0}).& & {}\\ \end{array}$

Then

$\displaystyle\begin{array}{rcl} d_{g}(x,\partial M)& =& \min _{z\in \partial M}r_{x}(z) = r_{x}(z_{0}) \\ & =& r_{y}(z_{0}) =\min _{z\in \partial M}r_{y}(z) = d_{g}(y,\partial M),{}\end{array}$

(20)

and $z_{0} \in \partial M$ is a nearest boundary point to x. Let μ _xbe the shortest path from z ₀ to x. Then, the path μ _xis a geodesic from x to z ₀ which intersects $\partial M$ first time at z ₀. By using the first variation on length formula, we see that μ _xhas to hit to z ₀ normally; see [22]. The same considerations are true for the point y with the same point z ₀. Thus, both x and y lie on the normal geodesic $\gamma _{z_{0},\nu }(s)$ to $\partial M$ . As the geodesics are unique solutions of a system of ordinary differential equations (the Hamilton–Jacobi equation (14)), they are uniquely determined by their initial points and directions, that is, the geodesics are non-branching. Thus, we see that

$\displaystyle\begin{array}{rcl} x =\gamma _{z_{0}}(s_{0}) = y,& & {}\\ \end{array}$

where $s_{0} = r_{x}(z_{0}) = r_{y}(z_{0})$ . Hence, $R: M \rightarrow C(\partial M)$ is injective.

Next, we consider the condition (19) for R: M → R(M). Let z ∈ M and consider closed sets $X_{j} =\{ x \in M:\ d_{C}\big(R(x),R(z)\big) \leq j\}$ , $j \in \mathbb{Z}_{+}$ . Then for x ∈ X _j, we have by definition (11) of the metric d _Cthat

$\displaystyle\begin{array}{rcl} d_{g}(x,Q_{0}) \leq j + d_{g}(z,Q_{0}),& & {}\\ \end{array}$

implying that the sets X _j, $j \in \mathbb{Z}_{+}$ are compact. Clearly, $\bigcup _{j\in \mathbb{Z}_{+}}X_{j} = X$ . Let next $Y _{j} = R(X_{j}) \subset Y = R(M)$ and $p = R(Q_{0}) \in R(M)$ . Then for $r_{x} \in Y \setminus Y _{j}$ , we have

$\displaystyle\begin{array}{rcl} d_{C}(r_{x},p)& \geq & r_{x}(Q_{0}) - p(Q_{0}) = d_{g}(x,Q_{0}) {}\\ & \geq & j - d_{g}(z,Q_{0}) - C_{0} \rightarrow \infty \mbox{ as }j \rightarrow \infty {}\\ \end{array}$

and thus the condition (19) is satisfied. As R: M → R(M) is a continuous, bijective map, it implies that R: M → R(M) is a homeomorphism.

Next we introduce a differentiable structure and a metric tensor, g _R, on R(M) to have an isometric diffeomorphism

$\displaystyle\begin{array}{rcl} R: (M,g) \rightarrow (R(M),g_{R}).& &{}\end{array}$

(21)

Such structures clearly exists – the map R pushes the differentiable structure of M and the metric g to some differentiable structure on R(M) and the metric g _R: = R _∗ g which makes the map (21) an isometric diffeomorphism. Next we construct these coordinates and the metric tensor in those on R(M) using the fact that R(M) is known as a subset of $C(\partial M)$ .

We will start by construction of the differentiable and metric structures on $R(M)\setminus R(\omega )$ , where ω is the cut locus of M with respect to $\partial M$ . First, we show that we can identify in the set R(M) all the elements of the form r = r _x ∈ R(M) where $x \in M\setminus \omega$ . To do this, we observe that r = r _xwith $x =\gamma _{z,\nu }(s),\,s <\tau (z)$ if and only if:

$r(\cdot )$ has a unique global minimum at some point $z \in \partial M$ ;

There is $\tilde{r} \in R(M)$ having a unique global minimum at the same z and $r(z) <\tilde{ r}(z)$ . This is equivalent to saying that there is y with $r_{y}(\cdot )$ having a unique global minimum at the same z and r _x(z) < r _y(z).

Thus, we can find $R(M\setminus \omega )$ by choosing all those r ∈ R(M) for which the above conditions (1) and (2) are valid.

Next, we choose a differentiable structure on $R(M\setminus \omega )$ which makes the map $R: M\setminus \omega \rightarrow R(M\setminus \omega )$ a diffeomorphism. This can be done by introducing coordinates near each $r^{0} \in R(M\setminus \omega )$ . In a sufficiently small neighborhood $W \subset R(M)$ of r ⁰ the coordinates

$\displaystyle\begin{array}{rcl} r\mapsto \big(Y (r),T(r)\big) = \left (y(\mbox{ argmin}_{z\in \partial M}r),\min _{z\in \partial M}r\right )& & {}\\ \end{array}$

are well defined. These coordinates have the property that the map $x\mapsto \big(Y (r_{x}),T(r_{x})\big)$ coincides with the boundary normal coordinates (17) and (18). When we choose the differential structure on $R(M\setminus \omega )$ that corresponds to these coordinates, the map

$\displaystyle\begin{array}{rcl} R: M\setminus \omega \rightarrow R(M\setminus \omega )& & {}\\ \end{array}$

is a diffeomorphism.

Next we construct the metric g _Ron R(M). Let $r^{0} \in R(M\setminus \omega )$ . As above, in a sufficiently small neighborhood W ⊂ R(M) of r ⁰, there are coordinates $r\mapsto X(r):=\big (Y (r),T(r)\big)$ that correspond to the boundary normal coordinates. Let $(y^{0},t^{0}) = X(r^{0})$ . We consider next the evaluation function

$\displaystyle\begin{array}{rcl} K_{w}: W \rightarrow \mathbb{R},\quad K_{w}(r) = r(w),& & {}\\ \end{array}$

where $w \in \partial M$ . The inverse of $X: W \rightarrow \mathbb{R}^{n}$ is well defined in a neighborhood $U \subset \mathbb{R}^{n}$ of (y ⁰, t ⁰), and thus we can define the function

$\displaystyle\begin{array}{rcl} E_{w} = K_{w} \circ X^{-1}: U \rightarrow \mathbb{R}& & {}\\ \end{array}$

that satisfies

$\displaystyle\begin{array}{rcl} E_{w}(y,t):= d_{g}\left (w,\gamma _{z(y),\nu (y)}(t)\right ),\quad (y,t) \in U,& &{}\end{array}$

(22)

where $\gamma _{z(y),\nu (y)}(t)$ is the normal geodesic starting from the boundary point z(y) with coordinates $y = (y^{1},\ldots,y^{n-1})$ and ν(y) is the interior unit normal vector at y.

Let now g _R = R _∗ g be the push-forward of g to $R(M\setminus \omega )$ . We denote its representation in X-coordinates by g _jk(y, t). Since X corresponds to the boundary normal coordinates, the metric tensor satisfies

$\displaystyle\begin{array}{rcl} g_{mm} = 1,\quad g_{\alpha m} = 0,\quad \alpha = 1,\ldots,n - 1.& & {}\\ \end{array}$

Consider the function E _w(y, t) as a function of (y, t) with a fixed w. Then its differential, dE _wat point (y, t) defines a covector in $T_{(y,t)}^{{\ast}}(U) = \mathbb{R}^{n}$ . Since the gradient of a distance function is a unit vector field, we see from (22) that

$\displaystyle{\|dE_{w}(y,t)\|_{(g_{jk})}^{2}:= \left ( \frac{\partial } {\partial t}E_{w}(y,t)\right )^{2} +\sum _{ \alpha,\beta =1}^{n-1}(g_{ R})^{\alpha \beta }(y,t)\frac{\partial E_{w}} {\partial y^{\alpha }} (y,t)\frac{\partial E_{w}} {\partial y^{\beta }} (y,t) = 1.}$

Let us next fix a point (y ⁰, t ⁰) ∈ U. Varying the point $w \in \partial M$ , we obtain a set of covectors $dE_{w}(y^{0},t^{0})$ in the unit ball of $\left (T_{(y^{0},t^{0})}^{{\ast}}U,g_{jk}\right )$ which contains an open neighborhood of $(0,\ldots,0,1)$ . This determines uniquely the tensor $g^{jk}(y^{0},t^{0})$ . Thus, we can construct the metric tensor in the boundary normal coordinates at arbitrary $r \in R(M\setminus \omega )$ . This means that we can find the metric g _Ron $R(M\setminus \omega )$ when R(M) is given.

To complete the reconstruction, we need to find the differentiable structure and the metric tensor near R(ω). Let $r^{(0)} \in R(\omega )$ and $x^{(0)} \in M^{\mathrm{int}}$ be such a point that $r^{(0)} = r_{x^{(0)}} = R(x^{(0)})$ . Let z ₀ be some of the closest points of $\partial M$ to the point x ⁽⁰⁾. Then there are points $z_{1},\ldots,z_{n-1}$ on $\partial M$ , given by $z_{j} =\mu _{z_{0},\theta _{j}}(s_{0})$ , where $\mu _{z_{0},\theta _{j}}(s)$ are geodesics of $(\partial M,g_{\partial M})$ and $\theta _{1},\ldots,\theta _{n-1}$ are orthonormal vectors of $T_{z_{0}}(\partial M)$ with respect to metric $g_{\partial M}$ and s ₀ > 0 is sufficiently small, so that the distance functions $y\mapsto d_{g}(z_{i},y)$ , $i = 0,1,2,\ldots,n - 1$ form local coordinates $y\mapsto \big(d_{g}(z_{i},y)\big)_{i=0}^{n-1}$ on M in some neighborhood of the point x ⁽⁰⁾ (we omit here the proof which can be found in [42, Lemma 2.14]).

Let now W ⊂ R(M) be a neighborhood of r ⁽⁰⁾ and let $\tilde{r} \in W$ . Moreover, let $V = R^{-1}(W) \subset M$ and $\tilde{x} = R^{-1}(\tilde{r}) \in V$ . Let us next consider arbitrary points $z_{1},\ldots,z_{n-1}$ on $\partial M$ . Our aim is to verify whether the functions $x\mapsto X^{i}(x) = d_{g}(x,z_{i}),$ $i = 0,1,\ldots,n - 1$ form smooth coordinates in V. As $M\setminus \omega$ is dense on M and we have found topological structure of R(M) and constructed the metric g _Ron $R(M\setminus \omega )$ , we can choose $r^{(j)} \in R(M\setminus \omega )$ such that $\lim _{j\rightarrow \infty }r^{(j)} =\tilde{ r}$ in R(M). Let $x^{(j)} \in M\setminus \omega$ be the points for which $r^{(j)} = R\big(x^{(j)}\big)$ . Now the function $x\mapsto \big(X^{i}(x)\big)_{i=0}^{n-1}$ defines smooth coordinates near $\tilde{x}$ if and only if for functions $Z^{i}(r) = K_{z_{i}}(r)$ , we have

$\displaystyle\begin{array}{rcl} & & \lim _{j\rightarrow \infty }\mbox{ det}\left (\big(g_{R}\big(dZ^{i}(r),dZ^{l}(r)\big)\big)_{ i,l=0}^{n-1}\right )\vert _{ r=r^{(j)}} \\ & & \qquad =\lim _{j\rightarrow \infty }\mbox{ det}\left (\big(g\big(dX^{i}(x),dX^{l}(x)\big)\big)_{ i,l=0}^{n-1}\right )\vert _{ x=x^{(j)}}\not =0.{}\end{array}$

(23)

Thus, for all $\tilde{r} \in W$ we can verify for any points $z_{1},\ldots,z_{n-1} \in \partial M$ whether the condition (23) is valid or not and this condition is valid for all $\tilde{r} \in W$ if and only if the functions $x\mapsto X^{i}(x) = d_{g}(x,z_{i}),$ $i = 0,1,\ldots,n - 1$ form smooth coordinates in V. Moreover, by the above reasoning, we know that any $r^{(0)} \in R(\omega )$ has some neighborhood W and some points $z_{1},\ldots,z_{n-1} \in \partial M$ for which the condition (23) is valid for all $\tilde{r} \in W$ . By choosing such points, we find also near $r^{(0)} \in (\omega )$ smooth coordinates $r\mapsto \big(Z^{i}(r)\big)_{i=0}^{n-1}$ which make the map R: M → R(M) a diffeomorphism near x ⁽⁰⁾.

Summarizing, we have constructed differentiable structure (i.e., local coordinates) on the whole set R(M), and this differentiable structure makes the map R: M → R(M) a diffeomorphism. Moreover, since the metric g _R = R _∗ g is a smooth tensor, and we have found it in a dense subset $R(M\setminus \omega )$ of R(M), we can continue it in the local coordinates. This gives us the metric g _Ron the whole R(M), which makes the map R: M → R(M) an isometric diffeomorphism. $\blacksquare$

In the above proof, the reconstruction of the metric tensor in the boundary normal coordinates can be considered as finding the image of the metric in the travel time coordinates.

Let us next consider the case when we have an unknown isotropic wave speed c(x) in a bounded domain $\Omega \subset \mathbb{R}^{n}$ . We will assume that we are given the set $\Omega$ and an abstract Riemannian manifold (M, g), which is isometric to $\Omega$ endowed with its travel time metric corresponding to the wave speed c(x). Also, we assume that we are given a map $\psi: \partial \Omega \rightarrow \partial M$ , which gives the correspondence between the boundary points of $\Omega$ and M. Next we show that it is then possible to find an embedding from the manifold M to $\Omega$ which gives us the wave speed c(x) at each point $x \in \Omega$ . This construction is presented in detail, e.g., in [42].

For this end, we need first to reconstruct a function σ on M which corresponds to the function c(x)² on $\Omega$ . This is done on the following lemma.