Stable reconstruction of simple Riemannian manifolds from unknown interior sources

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Stable reconstruction of simple Riemannian manifolds from unknown interior sources © 2023 The Author(s). Published by IOP Publishing Ltd Printed in the UK Published version de Hoop, Maarten V.; Ilmavirta, Joonas; Lassas, Matti; Saksala, Teemu de Hoop, M. V., Ilmavirta, J., Lassas, M., & Saksala, T. (2023). Stable reconstruction of simple Riemannian manifolds from unknown interior sources. Inverse Problems, 39(9), Article 095002. https://doi.org/10.1088/1361-6420/ace6c9 2023 Inverse Problems PAPER • OPEN ACCESS Stable reconstruction of simple Riemannian manifolds from unknown interior sources To cite this article: Maarten V de Hoop et al 2023 Inverse Problems 39 095002 View the article online for updates and enhancements. You may also like Connections on a non-symmetric (generalized) Riemannian manifold and gravity Stefan Ivanov and Milan Zlatanovi - Riemannian geometry of fibre bundles A A Borisenko and A L Yampol'skii - Holonomy groups of Lorentzian manifolds A S. Galaev - This content was downloaded from IP address 130.234.90.31 on 28/08/2023 at 09:50 Inverse Problems Inverse Problems 39 (2023) 095002 (43pp) https://doi.org/10.1088/1361-6420/ace6c9 Stable reconstruction of simple Riemannian manifolds from unknown interior sources Maarten V de Hoop1, Joonas Ilmavirta2, Matti Lassas3 and Teemu Saksala4,∗ 1Computational Applied Mathematics & Operations Research, Rice University, 6100 Main MS-134, Houston, TX 77005-1892, United States of America 2Department of Mathematics and Statistics, University of Jyväskylä, PO Box 35 (MaD), FI-40014 Jyväskylä, Finland 3Department of Mathematics and Statistics, University of Helsinki, PO Box 68 (Gustaf Hällströmin katu 2B), FI-00014 Helsinki, Finland 4Department of Mathematics, North Carolina State University, 2311 Stinson Drive, Raleigh, NC 27695-8205, United States of America E-mail: [email protected] Received 23 January 2023; revised 22 June 2023 Accepted for publication 12 July 2023 Published 28 July 2023 Abstract Consider the geometric inverse problem: there is a set of delta-sources in spacetime that emit waves travelling at unit speed. If we know all the arrival times at the boundary cylinder of the spacetime, can we reconstruct the space, a Riemannian manifold with boundary? With a finite set of sources we can only hope to get an approximate reconstruction, and we indeed provide a discrete metric approximation to the manifold with explicit data-driven error bounds when the manifold is simple. This is the geometrization of a seismological inverse problem where we measure the arrival times on the surface of waves from an unknown number of unknown interior microseismic events at unknown times. The closeness of two metric spaces with a marked boundary is measured by a labeled Gromov–Hausdorff distance. If measurements are done for ∗Author to whom any correspondence should be addressed. Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1361-6420/23/095002+43$33.00 © 2023 The Author(s). Published by IOP Publishing Ltd Printed in the UK 1 Inverse Problems 39 (2023) 095002 M V de Hoop et al infinite time and spatially dense sources, our construction produces the true Riemannian manifold and the finite-time approximations converge to it in the metric sense Keywords: inverse problem, Riemannian geometry, distance function, stability, discrete approximation (Some figures may appear in colour only in the online journal) 1. Introduction We study a geometric inverse problem arising from global seismology. Multiple sources go off at unknown times and we measure the arrival times, not knowing which arrivals belong to the same event. The aim is to extract as much information as possible about the planet, modeled here by a simple Riemannian manifold Mof dimension n=dim(M)⩾2, and the sources from this information. Mathematically, the problem boils down to this: Suppose we know the boundary of a Riemannian manifold and for some unknown collection of interior points we know the union of their boundary distance graphs, each shifted with an unknown offset. What can we say about the manifold and these special points from this data? A more detailed description of the data can be found in section 2below. The distances between the source points turn out to be determined exactly by the data (theorem 2), but there is no hope of reconstructing the manifold if the source set is finite. However, the source set can be regarded as a discrete metric space, and this space approximates the Riemannian manifold. It is not only close as a metric space, but we can also assign approximate boundary points with good accuracy (theorem 9). To describe the similarity of two metric spaces ‘with the same boundary’, we define a labeled Gromov–Hausdorff distance. This extends the classical Gromov–Hausdorff distance and compares both the similarity of the metric spaces and the sameness of the boundaries—with a fixed model space for the boundary. It is in the sense of this distance that our approximation is good, and the quality of the approximation can be estimated explicitly directly from data. The sources can be point sources in space time produced by a suitable stochastic process (proposition 12), for example. We present a method to look at the data on the boundary without knowing the manifold and reconstruct an approximate manifold—and get an explicit error bound on our reconstruction. This increases the applicability of geometric inverse problems to the messy cases with real data, but we do not attempt to optimize the bounds. When it comes to explicit geometric estimates, this paper should be seen as a proof of concept. We assume the source events to be discrete in the spacetime, but their projections to the space can be dense. In this case the approximate reconstruction is actually perfect, and even the smooth and Riemannian structure is determined (theorem 1). If measurements are made for finite but increasing time, the approximate reconstructions converge to the true manifold in the labeled Gromov–Hausdorff distance if the sources are spatially dense (theorem 11). Estimation of the quality of the approximation from data alone is lossy. To prove that it is not hopelessly so, we prove that by making the source points dense enough the observed density can be made arbitrarily small. There is an explicit estimate between the true and observed density in both directions (theorems 1and 10). We solve a problem stemming from physics by developing geometry rather than by applying ad hoc tricks. These developments include the introduction of the labeled Gromov–Hausdorff distance and various quantitative descriptions of simplicity of a Riemannian manifold. 2 Inverse Problems 39 (2023) 095002 M V de Hoop et al From a more applied point of view, the key advances are that we introduce of a new and more physically relevant problem and provide quantitative and fully data-driven error estimates. The precise definitions and statements of results are given in section 2. We will discuss the results and their context in section 3. The overall strategy of proofs is given in section 4and the actual proofs follow in subsequent sections. 2. Definitions and theorems 2.1. Arrival time data Let Mbe a Riemannian manifold with boundary. We denote by d:M×M→Rits Riemannian distance. The spacetime M×Rcomes with two natural projections, πonto Mand τonto R. Let S⊂M×Rbe a set of sources. The source s∈Sgoes off at the point π(s)at the time τ(s). The data measured is the arrival times of signals at the boundary, given as the set Q(S) = {(x,τ(s) + d(x,π(s))); x∈∂M,s∈S}. We emphasize that we know the set Q(S) but not the set S. We only observe the data as this point set without any labels to tell which arrivals correspond to the same sources. The arrival time function of a source s∈Sis the function as:∂M→Rgiven by as(x) = τ(s) + d(π(s),x). Let us denote the graph of the function asby G(as)⊂∂M×R. The data set can be rewritten as Q(S) = Ss∈SG(as). Let ϕ:∂M1→∂M2be a diffeomorphism between two boundaries. If A⊂∂M2×R, we define ϕ∗A={(ϕ−1(x),t); (x,t)∈A}⊂∂M1×R. The set Awould typically be the data set Q(S2)or a subset thereof. When there are two manifolds, we decorate all the objects with the subscripts 1 and 2 when clarity requires so. A subset of the spacetime M×Ris said to be discrete when it has no accumulation points. If Mis compact, it follows that any time slice M×[a,b]contains only finitely many points of a discrete set. Our results apply to so-called simple manifolds. A Riemannian manifold is called simple if it is compact with strictly convex boundary (in the sense of definiteness of the second fundamental form) and the exponential map at every point is a diffeomorphism on its maximal domain of definition. Consequently there are no conjugate points and any two points on the manifold, including the boundary, are joined by a unique geodesic depending smoothly on the endpoints. A key concept in our analysis is density. We say that a subset A⊂Xof a metric space Xis ε-dense if the balls B(a,ε)with a∈Acover X. Equivalently, every point in Xhas a point of A at a distance below ε. We set out to study whether and how a Riemannian manifold Mcan be determined by the set Q(S) and prior knowledge of the boundary ∂M. 2.2. Precise determination We begin by studying what can be determined from the data with no error. The source sets Si⊂Mi×Rcan be finite or infinite. However, the discreteness assumption implies that every Si∩(Mi×[−T,T]) is finite and so the source sets are at most countably infinite. Although the sources are not dense in the spacetime, the source points in space can be dense. In that case the same data determines the whole manifold up a diffeomorphism. 3 Inverse Problems 39 (2023) 095002 M V de Hoop et al Theorem 1. Let M1and M2be simple Riemannian manifolds so that the diameter of both manifolds is at most Cdiam >0and the sectional curvature is bounded from above by Csec+>0 so that CdiampCsec+< π. (1) Let ϕ:∂M1→∂M2be a Riemannian isometry. Suppose the two sets Si⊂int(Mi)×Rare discrete in Mi×R. Let Q(Si)⊂∂Mi×Rbe defined as above. If the set π1(S1)⊂M1is dense and Q(S1) = ϕ∗Q(S2), then (1) there is a smooth Riemannian isometry Φ: M1→M2so that ϕ= Φ|∂M1and (2) there is a bijection ξ:S1→S2between the sources so that ξ(p,t) = (Φ(p),t)(2) for all (p,t)∈S1. The condition (1) in most of our theorems always holds when the manifold is a priori known to have non-positive sectional curvature. See section 3.2.4 for details. The theorem immediately implies that the source times coincide: τ1(s) = τ2(ξ(s)) for all s∈S1. The main novelty of theorem 1is in allowing the data set Q(S) to be a union rather than an indexed collection of graphs of the arrival time functions. The rest of our results do not have a similar precedent, and our method of proof for theorem 1is new and constructive. See the discussion in section 3below for details. The next result concerns a situation with a finite amount of sources. Informally, the next theorem states that knowledge of ∂Mas a smooth manifold and the set Q(S) determines the distances between the sources, the time differences, and two kinds of distance differences at the boundary. More formally: Theorem 2. Let M1and M2be simple Riemannian manifolds and ϕ:∂M1→∂M2a diffeomorphism. Suppose the two sets Si⊂int(Mi)×Rare discrete in Mi×Rand #π(S1)⩾2. Let Q(Si)⊂∂Mi×Rbe defined as above. If Q(S1) = ϕ∗Q(S2), then (1) there is a bijection ξ:S1→S2, (2) d1(π1(s),π1(r)) = d2(π2(ξ(s)),π2(ξ(r))) for all s,r∈S1, (3) τ1(s)−τ1(r) = τ2(ξ(s)) −τ2(ξ(r)) for all s,r∈S1, (4) d1(π1(s),x)−d1(π1(r),x) = d2(π2(ξ(s)),ϕ(x)) −d2(π2(ξ(r)),ϕ(x)) for all x ∈∂M1and s,r∈S1, and (5) d1(π1(s),x)−d1(π1(s),y) = d2(π2(ξ(s)),ϕ(x)) −d2(π2(ξ(s)),ϕ(y)) for all x,y∈∂M1 and s ∈S1. While measuring for infinite time gives full uniqueness, it is not realistic. Therefore we consider measurements for a finite but increasing time and show that there is an approximate reconstruction and in the limit of infinite time it tends to the correct one in a suitable sense. To be able to state the results, we need to set up a way to compare the true manifold to a discrete approximation. Before pursuing this direction, we record a result that helps clean up the data by disentangling the union of graphs into separate graphs. The arrival time functions asare easily verified to be smooth when π(s)/∈∂Mand the manifold is simple. 4 Inverse Problems 39 (2023) 095002 M V de Hoop et al Proposition 3. Let M be a simple Riemannian manifold. Suppose S ⊂int(M)×Ris discrete in M ×R. Then the set Q(S) and the Riemannian structure of ∂M determine the set {as;s∈S}, of arrival time functions. Furthermore, if Ω⊂∂M×Ris any open set, then Q(S) determines the set of connected components of the sets G(as)∩Ω, s ∈S. 2.3. Labeled Gromov–Hausdorff distance This subsection is devoted to the metric geometry we will use to compare our approximate reconstructions to the true manifold. Definition 4. Let Xand Ybe compact metric spaces and Lany set (understood as a set of labels). Let α:L→Xand β:L→Ybe any functions. We define the labeled Gromov– Hausdorff distance between (X,α)and (Y,β)to be dL GH(X,α;Y,β) = inf{dZ H(f(X),g(Y)) + sup ℓ∈L dZ(f(α(ℓ)),g(β(ℓ))); Zis a compact metric space, f:X→Zand g:Y→Zare isometric embeddings}. Here dZ His the Hausdorff distance on the metric space Z. If the set of labels is empty, then this distance reduces to the usual Gromov–Hausdorff distance; in this context we understand the supremum of an empty set to be zero. We will be comparing our discrete reconstruction to the manifold Mwith labels given by the known set L=∂M, labeled by the inclusion ι:∂M→M. To make it meaningful to state that the labeled Gromov–Hausdorff distance is small, we must ensure that this is a well-behaved concept of distance. The point of this definition is that in addition to getting a good metric approximation of the manifold in the sense of Gromov–Hausdorff distance, we want to get knowledge of how well the known boundary ∂Msits inside the approximation. This is what the concept is designed to do. Proposition 5. The labeled Gromov–Hausdorff distance is symmetric and satisfies the triangle inequality. Moreover, dL GH(X,α;Y,β) = 0 if and only if there is an isometry h:X→Y so that h ◦α=β. 2.4. Quantitative simplicity We will make approximate reconstructions of the whole manifold, and we can quantify the error of the reconstruction precisely once we have a priori bounds on curvature, diameter, and similar geometric properties. This set of properties might not be minimal in any way, but the following are the estimates we shall need. By proposition 7all simple manifolds do indeed have bounded geometry in the sense of the next definition. Definition 6. We say that a manifold Mwith boundary has bounded geometry with the constants CJF,CSFF,Cdiam,Csec−,Csec+,Cexp,Cdist,CH1, and CH2 (each of them >0) if the following properties hold: 5 Inverse Problems 39 (2023) 095002 M V de Hoop et al (1) The diameter of Mis at most Cdiam. (2) The sectional curvature is always in the range [−Csec−,Csec+]. (3) For any x∈Mand η1,η2∈TxMwe have |η1−η2|⩽Cexpd(expxη1,expxη2) as long as the exponentials are defined. (4) Every Jacobi field Jwith J(0) = 0 along any unit speed geodesic satisfies t 2∂t|J(t)|2⩽CJF |J(t)|2 for all t>0 for which the geodesic is defined. (5) If h1and h2are, respectively, the first and second fundamental forms of the boundary ∂M, then h2⩽CSFFh1. (6) The distances on Mand along ∂Msatisfy dM(x,y)⩽d∂M(x,y)⩽CdistdM(x,y) for all x,y∈∂M. (7) Take any unit speed geodesic γon Mand let ρbe the distance function ρ(x) = d(x,γ(0)). For any t>0, let w∈Tγ(t)Mbe any vector orthogonal to ˙γ(t). Then the Hessian Hρof ρ satisfies hw,Hρwi⩾|w|2(CH1t−1−CH2). We have to point out, however, that bounded geometry does not imply simplicity. This is because the estimates below do not ensure strictly convex boundary, and indeed any smooth domain in a simple manifold satisfies all the estimates. Strict convexity of the boundary is important for our results, but it is used in a qualitative rather than a quantitative way. See section 3.2.2 and proposition 3for the use of convexity. If item (5) is improved from h2≲h1 to h1≲h2≲h1, then strict convexity does indeed follow. We will derive a number of other estimates from those of definition 6with constants depending on the ones listed here. See section 9for details. The constants appearing in the derived estimates are all given in appendix B. We will refer to the appendix whenever a new constant is introduced in a claim. All claims are given with simple constants, so there will be a number of relations between the various constants as is clear from the appendix. It should be noted in part (7) of the definition that wis a vector and Hρwis a covector. The inner product is their duality pairing, which we find most convenient to write as hw,Hρwi instead of (Hρw)(w)or Hρ(w,w). Having these estimates and constants is not an added assumptions, but merely a quantification of the simplicity of the manifold M. Our next result justifies that the constants of definition 6can well be called constants of quantitative simplicity. Proposition 7. Every simple manifold has bounded geometry in the sense of definition 6for some constants. 2.5. Approximate determination Theorems 1and 2state what can be determined exactly from our data. We now turn to studying how well the whole manifold can be reconstructed with finite data. 6 Inverse Problems 39 (2023) 095002 M V de Hoop et al The source set is as before. The measurements at the boundary start at time t=0 and continue up to time t=T. From this finite amount of data we can then draw an approximate conclusion about the geometry of our manifold. The error can be quantified very concretely in terms of the geometric bounds. The only points we know on the manifold are the finite number of sources that have produced a boundary signal on the time interval [0,T]and the a priori known boundary itself. Therefore what we have is a discrete approximation of the smooth Riemannian manifold. In addition to getting a metric approximation of the manifold itself, we need to approximately embed the known boundary ∂Minto the approximation, and so we will use the Gromov– Hausdorff distance labeled by ∂M. Let us denote the set of spatial source points by P=π(S). To each p∈Pis associated a boundary distance function rp:∂M→Rgiven by rp(x) = d(x,p). Let c(p)⊂∂Mbe the set of critical points of this function. By compactness #c(p)⩾2. Consider a simple manifold of bounded geometry in the sense of definition 6. At x∈c(p) let λ(p,x)be the smallest eigenvalue of the Hessian of rpat x. This requires that ∂Mis known as a Riemannian manifold, not just as a smooth one. For y∈c(p)we define E(p,y) = (CJF λ(p,y)−CSFF ,when λ(p,y)>CSFF ∞,otherwise (3) and for all x∈∂M E(x) = inf p∈P,y∈c(p)(E(p,y) + d∂M(x,y)).(4) Finally, we let E=sup x∈∂M E(x).(5) These estimate the following geometric quantities from above: •E(p,y)estimates the distance d(p,y). •E(x) estimates the distance from xto the nearest source point. •Eestimates how far at worst a boundary point can be from a source point. These three types of Es will be all determined by the data and the constants of bounded geometry. Observe that if the Hessians of the distance functions are never large enough, all these Es may well be infinite. For any p,q∈Pand r,s>0 we define the lentil Lp,q r,s=B(p,r)∩B(q,s). The lentil is said to have thickness δp,q r,s=r+s−d(p,q). One can check that on a simple manifold the thickness is the length of the segment of the geodesic γp,qbetween the two points that lies within the lentil, provided that r,s<d(p,q). To see this, observe that a point γp,q(t)is in the lentil if and only if t<rand d(p,q)−t<s. Our inversion procedure starts with preprocessing the data. Proposition 3gives one aspect of it, and the next one says that the derived quantities introduced above are determined by the data. Proposition 8. Let M1and M2be simple Riemannian manifolds. Let ϕ:∂M1→∂M2be a smooth isometry. Suppose that the two sets Si⊂int(Mi)×Rare discrete in Mi×Rand 7 Inverse Problems 39 (2023) 095002 M V de Hoop et al that Randers metrics Fi=Gi+βiindexed with i∈{1,2}with simple and boundary rigid Riemannian norm Gi(x,v) = pgij(x)vivjand closed one-form βi, have the same boundary distance function if and only if G1= Ψ∗G2for some boundary fixing diffeomorphism Ψ: M→ Mand β1−β2=dϕfor some smooth function ϕvanishing on ∂M. It is worth of mentioning that analogous results have been presented earlier on a Riemannian manifold in the presence of a magnetic field [2,10]. The sphere data described in connection to the Riemannian results of [11] above has also been studied on Finsler manifolds [12]. Knowledge of the spheres uniquely determine the fundamental tensor and the curvature operator along any geodesic passing through the known domain, but in contrast to Riemannian geometry this information is insufficient for a full reconstruction of the universal cover. 3.1.3. Geophysical literature. In fact, the inverse problem considered here can be directly related to seismology. The sources correspond with microseismic events with unknown locations and origin times while the metric corresponds with the wave speed. Here, we consider only one wave speed associated with elastic P-waves, but incorporating a second wave speed associated with elastic S-waves is quite straightforward. In the past decade there has been extensive research on the joint recovery of wave speed and event locations, with mining, geothermal and hydraulic-fracturing applications of induced seismicity but also in studies of the crustal structure. Here, we present a comprehensive analysis providing fundamental insight in the feasibility of succeeding in this, while focusing on recovering the wave speed. We use geometric data rather than wave fields; in applications, the extraction of these have been routine, for example, with template matching [18]. Initial empirical studies [5,27] assumed a simple wave speed model varying in one coordinate only and one or two strings of receivers (in boreholes) penetrating the manifold. Strategies have been broadly based on optimization, such as variations of gradient descent [22,50], with intertwining updating [42] and a neighborhood algorithm [55] employing different optimization criteria [52,59] or via intermediate (approximate) interior wave-field recovery [51]. A basic statistical, Bayesian framework has been developed in parallel [61]. 3.2. Discussion of technical assumptions In this subsection we present a discussion of many technical assumptions we made in our theorems. Our focus here is on the key phenomenology of having multiple sources with unknown times and we have chosen to simplify other aspects of the problem. 3.2.1. Density of source points. We assume in all our results that the sources are discrete in the spacetime M×R. This is physically reasonable and also practical. If the sources were dense or even accumulated somewhere, proposition 3would become far more complicated. However, there is little hope of full uniqueness with finitely many source points. The sources can be discrete in spacetime but dense in space, and this is indeed the setting of theorems 1 and 11 in which the Riemannian manifold is determined uniquely either from full time data or asymptotically from increasing time data. 3.2.2. Full data and convexity. All our results concern full data in the sense that arrivals are recorded on the whole boundary ∂M. The full boundary is crucial for proposition 3which disentangles the data into a collection of graphs, and we make so heavy use of differential tools 14 Inverse Problems 39 (2023) 095002 M V de Hoop et al Figure 1. A domain where partial data is insufficient. We split the domain Minto two pieces M1and M2with respect to the line (red dotted line) that is normal to ∂Mat x0∈ ∂M(blue dot). Then we choose a domain Γ⊂∂M1(red arch) so that any minimizing curve joining a point on Γand a point in M2touches the boundary near x0. The curve P⊂Mis any involute of the boundary, meaning that the distance from all points on P to x0is the same. on the boundary throughout the paper that a discrete subset of the boundary is beyond current reach. Let us construct an explicit example of a surface Mand a subset Γ⊂∂Mso that our results fail with data recorded only on Γ. Every pair of points on a smooth compact Riemannian manifold with boundary is connected by a C1-smooth distance minimizing curve [1]. We choose our a manifold to be the horseshoe-shaped domain of figure 1. Because Pis an involute of the boundary as in the figure, d(x,p) = d(x,q)for any x∈Γand p,q∈P. Therefore from the point of view of our data, the set Pappears to collapse to a point. Worse, if Γhappens to be an involute as well, then all the distance functions to points p∈M2are constants in Γ. As the unknown origin times produce unknown constant offsets to the boundary distances, all points in M2look alike when seen from an involutive Γin this sense. The problems of figure 1also illustrate the problems lack of convexity of the boundary may cause. Similar problems may arise in higher dimensions if data is recorded on too small sets. Let the manifold Mbe the closed unit ball in Rnand consider partial data on Γ = ∂M∩H, where H⊂Rnis a hyperplane through the origin. Let p1and p2be two points in int(M)\Hsituated symmetrically about H. Then the boundary distance functions of p1and p2agree on Γ. The sets Mand Γ⊂∂Mhave a reflection symmetry which leaves the data invariant, making it impossible to distinguish the two sides. 3.2.3. Infinitely many lentils. The crucial condition for the estimates in theorem 9was that all of the lentils meet the spatial source set P. There is an infinite family of lentils, so we have a large number of conditions to check from our data. The lentils are open and cover a certain compact subset of M(excluding a layer near the boundary), so in fact using a finite number 15 Inverse Problems 39 (2023) 095002 M V de Hoop et al will suffice. However, it is not easy to identify a covering collection of lentils—for which one would then easily check whether they contain source points—from data. 3.2.4. Conjugate points and boundary sources. If there is a source point on the boundary, the graph of the distance function is singular at that point. The separation of the set Q(S) into graphs in proposition 3relies on smoothness. If several corners happen to coincide, especially if n=2, it may be difficult to choose how to continue the graphs correctly. This problem becomes far worse if there are conjugate points. They also cause the boundary distance function to be non-smooth but can do so for several different points. This makes both disentangling the data into graphs and the analysis of those graphs substantially more complicated. The condition (1) on the constants of definition 6is assumed in theorems 1,9–11. This condition is used in lemma 24 and proposition 30 when the manifold is compared to model manifold with constant sectional curvature. For the comparison to work, we need the model manifold of the same diameter to not have conjugate points, and this is exactly what the condition ensures. The condition (1) always holds in negative curvature. If the manifold is known to have nonpositive sectional curvature and an explicit bound on the diameter, one can simply choose the sectional curvature upper bound Csec+>0 to be small enough to satisfy (1). The condition also holds on all simple manifolds of constant sectional curvature. The constants Caof (37), Cbof (38), and Ccof (39), as well as all constants derived from them, become worse when CdiampCsec+gets close to π. 3.2.5. Constants of quantitative simplicity. The constants of quantitative simplicity have to be the same for both manifolds Miin theorem 9. The set Γiand the approximate boundary inclusion map αdepend on these constants. The distances between the source points are determined irrespective of the constants by theorem 2, but the constants have an effect on the approximate boundary structure and the estimates on the quality of the reconstruction. 3.2.6. Scaling of small quantities. Consider the case when all the various epsilons are very small. In the setting of theorem 9we have δ≈ε2≈ε2and in theorem 10 we have δ≈ε2≲ ε2≈ˆε. Ideally, we would expect to see no second powers so that all concepts of density—the size of lentils, the true density, the observed quantity ε2, the density estimated from data—to be bi-Lipschitz equivalent to each other. The difference in scaling is all due to the lentils being shaped so that outer radius ≈√inner radius. This scaling is easy to verify in Euclidean geometry, and the relevant Riemannian aspects are covered in lemmas 24 and 25, and proposition 30. 3.2.7. Generalization to simple Berwald manifolds. As discussed above, our methods depend on the smoothness of the distance function. Thus these techniques are not applicable for general Riemannian metrics and relaxing the simplicity assumption would likely require different methods. There is however hope to apply these techniques to study analogous questions in the case of certain simple Finsler metrics. 16 Inverse Problems 39 (2023) 095002 M V de Hoop et al Any constant speed geodesic γof a smooth connected and complete Finsler manifold (N,F) is given by the solution of the geodesic equations, in local coordinates, ¨γk(t) + Γk ij(γ, ˙γ)˙γi(t)˙γj(t) = 0.Here Γk ij(x,v)are the coefficients of the Chern connection on the tangent bundle [3]. Due to the dependence of the directional variable v∈TxNa reverse curve of a Finsler geodesic does not need to be a geodesic. If the Chern connection coefficients are directionally independent the Finsler metric Fis called Berwald. Thus Berwald manifolds have a well defined canonical covariant differential operator for tensor fields of any order. Therefore Berwald metrics carry a Gauss type of formula to connect boundary and interior Hessians for the distance function. This is not true for a general Finsler metrics and the lack of Gauss formula would cause issues in several parts of our proofs. It is worth of mentioning that Berwald metrics are not just a theoretical curiosity but have a connection to linear elasticity. It has been observed that elastic Finsler metrics, arising from transverse isotropic medium in weak anisotropy, are actually Berwaldian [60]. 4. Plan of proofs Our data is given in the form of a set, a union of graphs, and we begin by disentangling it into the separate graphs. We prove proposition 3to this effect in section 5. In section 6we turn to finding distances between source points and other information that the data gives exactly. With the help of proposition 3, we may start with the knowledge of the arrival time functions as. Each graph corresponds to a unique source point, giving us the first claim. Computing suitable differences of arrival times, especially along the unique geodesic between two given source points, we cancel out unwanted contributions and find the distance between the two points. Once we have these first two claims, the rest of theorem 2follows straightforwardly. In section 7we depart from differential geometry to the geometry of metric spaces and study the labeled Gromov–Hausdorff distance in more detail, including the proof of proposition 5. The most substantial part of our proofs is in estimating the density of the sources, and we will do this in section 9using the explicit bounds on geometric quantities found in section 8 — including the proof of proposition 7. The first task is to estimate the density of sources near the boundary, which practically amounts to estimating how near the boundary a given source point is. This is done in terms of the curvature of the graph of the arrival time function; when the point is very close to the boundary, the Hessian of the said function blows up. We will prove that d(p,x)⩽E(p,x)for any source point p∈Pand a critical point x∈∂Mof its arrival time function. This gives an estimate in terms of data for how close to any given boundary point must there be a source point. The initial data is first processed into the form of E(p,x)and quantities derived from it, and we verify in proposition 8that the necessary auxiliary quantities are indeed determined by our data. To get a density estimate deep inside the manifold, a different set of tools must be employed. A key tool is the lentil, an intersection of two metric balls. Whether a source point belongs in a lentil defined by two other source points is entirely decided by the distances between the three points, and by theorem 2this information is indeed determined by the data. We will show that with a suitable choice of δ > 0 (as given in theorem 9) the lentils cover the deep interior of the manifold but are small. This requires estimates in two directions for the lentils so that they are large enough to cover the relevant subset of Mwithout gaps but are small enough so that they provide a useful density estimate. If every point is contained in a lentil and that lentil contains 17 Inverse Problems 39 (2023) 095002 M V de Hoop et al a source point, then the density of sources is bounded by the diameters of the lentils. The more complicated side of lentils is to ensure that they cover enough of the manifold, and this relies on a number of estimates based on bounded geometry. These density estimates constitute a proof of theorem 9. Density estimates turn out to be simpler in the reverse direction, showing that the methods used in the proof of theorem 9are fine enough to give a decent estimate of density. These reverse estimates prove theorem 10. In section 10 we turn our attention to matters of convergence. The proof of theorem 11 is mostly based on theorems 9and 10; as Tincreases, the source points become denser and denser both in reality and in terms of data-driven estimates, and thus the quality of the discrete approximation improves. Theorem 1follows a similar path and can be regarded essentially as a corollary of theorem 11, and thus it has a relatively simple proof now that all the tools are available. The metric convergence results show that the two manifolds of theorem 1are isometric as metric spaces, and by the Myers–Steenrod theorem the isometry is in fact smooth. Proposition 12 gets its proof, using similar ideas and basic properties of Poisson point processes, at the very end. 5. Separation of graphs on a bundle Let us denote the scalar second fundamental form by h2and the boundary metric (the first fundamental form) by h1. If p∈int(M),ρ:M→Ris the distance function ρ(x) = d(p,x)at p then the corresponding boundary distance function r:∂M→R, is just the restriction of ρon the boundary ∂M. By the Gauss formula for hypersurfaces (see e.g. [40, theorem 8.13 a]) and the definition of Hessians we have that Hr(x) = Hρ(x) + h2(x).(14) We will make use of this throughout the paper. Our method of proving proposition 3is by means of lifting the data from the boundary ∂M×Rof the spacetime to a suitable bundle where the graphs do not intersect. We start by showing that when s6=s′, then the graphs of asand as′can only have up to first order tangency when they intersect. Therefore 2-jets will separate them. Lemma 13. Let M be a simple Riemannian manifold and s,s′∈int(M)×Rbe two distinct points in the spacetime. Then if as(x) = as′(x)at some point x ∈∂M, then the gradient or the Hessian of the two functions will differ at x. Proof. Take any two s,s′∈Sso that asand as′, their differentials, and their Hessians agree at a point x∈∂M. We aim to show that s=s′. First, observe that the differential das(x)is the covector corresponding to the tangential component of the velocity of the unit speed geodesic from π(s)to x. Therefore das(x) = das′(x)implies that p:= π(s)and p′:= π(s′)are on the same geodesic starting from x. Let γbe this unit speed geodesic starting at x. We have γ(t0) = pand γ(t′ 0) = p′for some t0,t′ 0>0. Consider the Hessian H:= Hρpof the interior distance function ρp:M→Rdefined by ρs(x) = d(x,p)and its counterpart H′corresponding to p′. If a normal Jacobi field Jalong γsatisfies HJ(t) = −DtJ(t)for some t<t0, then J(t0) = 0, (see for instance [40, proposition 11.2]). Similarly, H′J(t) = −DtJ(t)for any t<t′ 0implies J(t′ 0) = 0. 18 Inverse Problems 39 (2023) 095002 M V de Hoop et al Both Hessians satisfy H(′)˙γ(t) = −˙γ(t)for t<min(t0,t′ 0). As the geodesic meets ∂Mtransversally at xand the boundary Hessians of asand as′agree there, we have in fact the equality H(x) = H′(x)for the Hessian operators on TxM. This follows from (14). Now take any nonzero w∈TxMnormal to ˙γ(0)and let Jbe the Jacobi field along γwith J(0) = wand DtJ(0) = Hw=H′w. As observed above, this implies J(t0) = 0 and J(t′ 0) = 0. Due to the lack of conjugate points p=γ(t0) = γ(t′ 0) = p′. As the two source points are equal distance t0=t′ 0from xand the two arrivals are at the same time, we also have τ(s) = τ(s′). This concludes the proof of s=s′. Proof of proposition 3.Let us abbreviate Q:= Q(S)⊂∂M×R. As S⊂M×Ris discrete and Mis compact, there are only finitely many points of Sin any interval M×[a,b]. By simplicity there is no geodesic longer than the diameter of the manifold. Thus Qwritten as Q=[ s∈SG(as) is a locally finite union on ∂M×R. Let σ(Q)be the set of ‘smooth points of Q’, where Qis given locally as a single graph. By definition σ(Q)⊂Qis open. It is also dense, for otherwise there would be two graphs that coincide in an open set, contradicting lemma 13. Let Ebe the bundle over ∂M×Rwith fibers E(x,t)=Tx∂M×T⊗2 x∂M. This Eis the product of the bundles Tx∂Mand Tx∂M⊗Tx∂Mpulled to ∂M×Rover the projection ∂M×R→∂M. For any smooth function h:U→Rdefined in an open set U⊂∂Mwe call the lift of its graph the subset LG(h)⊂Eso that Lx,h(x)G(h)=(∇h(x),∇2h(x)) ∈Tx∂M×T⊗2 x∂M. This defines a smooth submanifold of E. In this way we lift all of σ(Q)into a submanifold Lσ(Q)of E. Since σ(Q)⊂Qis dense, we have Lσ(Q) = [ s∈S LG(as).(15) There are no second order intersections by lemma 13, so the smooth submanifolds L(G(as)) ⊂ Eare pairwise disjoint. Therefore Qdetermines the set {LG(as)⊂E;s∈S}. Projecting from Edown to ∂M×Rgives the graphs G(as)and thus also the functions as, proving the first part of the claim. For the second claim we simply localize (15) to be over Ω⊂∂M×R, and we get a disjoint union of the sets L[G(as)∩Ω]. A single graph G(as)may be cut into several pieces by this procedure and we do not necessarily know which pieces correspond to the same source. Therefore what we obtain is the collection of connected components of graphs restricted to Ω. It was proven in [31] that the boundary distance function rx:∂M→Rof x∈int(M)determines the point xuniquely. The next lemma improves this slightly: the boundary distance function modulo constants is enough for this uniqueness. 19 Inverse Problems 39 (2023) 095002 M V de Hoop et al Lemma 14. Let M be a simple Riemannian manifold and x,y∈int(M)any two distinct points. Denote by rxand rythe boundary distance functions ∂M→Rfrom these points. The difference function rx−rycannot be a constant function. Proof. Suppose that the difference function rx−rytakes the constant value c∈Rdespite x6=y. If c=0, the boundary distance functions coincide and thus by [31]x=y, yields a contradiction. Therefore we may suppose that c6=0. Let γbe the maximal geodesic through the points xand ywith endpoints x′,y′∈∂Mordered so that x′is closer to xand y′to y. By assumption we have d(x′,x) = rx(x′) = ry(x′) + c=d(x′,y) + c. Due to simplicity and the order of the points, we have d(x′,y) = d(x′,x) + d(x,y). Therefore 0=d(x,y) + c. The same calculation with the roles of xand yreversed shows that 0 =d(x,y)− c. These two together imply that d(x,y) = c=0, which is a contradiction. 6. Exact observables Proof of theorem 2.Due to proposition 3the data determines the arrival time functions as(x) = d(π(s),x) + τ(s)at all points x∈∂Mfor every s∈S. It also follows from the same proposition that distinct sources have distinct arrival time functions. Therefore now that the two manifolds M1and M2and their source sets S1and S2have the same data modulo identification by ϕ, we have a bijection between the source sets given by identifying the corresponding arrival time functions which coincide. This proves claim (1). We will prove all subsequent claims by describing how the arrival time functions determine the quantities in question. As these functions coincide on the two manifolds, the reconstructed quantities coincide. We therefore drop the subscripts and work with a single manifold. We can think of the source points being indexed by their arrival time graphs. To simplify notation, we denote the source points as ps=π(s)and source times as ts=τ(s). The arrival time functions now have the form as(x) = d(ps,x) + ts. Two source points ps and prcoincide if and only if as−aris a constant function—this is a straightforward corollary of lemmas 13 and 14. We pass to a subset of Sso that each source point psis only present once; the claims extend easily to the duplicated source points. For any s,r∈Swe define the functions frs :∂M×∂M→Rby frs(x,y) = ar(x)−as(y). These functions are determined by the data. Taking s=rproves claim (5). From now on we suppose that s6=r. The differentials dasand daragree at x∈∂Mif and only if the geodesics joining xto psand prstart in the same direction at x. Therefore the differentials agree at exactly two points, the endpoints of the unique maximal geodesic through the points ps and pr. Let these boundary points be xand y. Depending on how the four points (ps,pr,x,y) are ordered on the geodesic, we have frr(x,y)−fss(x,y) = d(pr,x)−d(ps,x) + d(ps,y)−d(pr,y) =±2d(pr,ps).(16) Thus claim (2) is given by d(pr,ps) = 1 2|frr(x,y)−fss(x,y)|. 20 Inverse Problems 39 (2023) 095002 M V de Hoop et al This determination of distances between source points will play a central role in the proofs of the other theorems. By switching the points xand yif needed (depending on the sign in (16)), we can assume to have d(pr,y) = d(pr,ps) + d(ps,y). Then for any z∈∂M, we have frr(z,y)−fss(z,y) + d(pr,ps) = d(pr,z)−d(ps,z). As everything on the left-hand side is determined by the data, so is then the right-hand side, and we have obtained the function frs :∂M→Rgiven by frs(z) = d(pr,z)−d(ps,z) for all s,r∈S. This proves claim (4). Finally, we note that for any z∈∂Mwe have tr−ts=frs(z,z)−frs(z) and claim (3) follows. As mentioned, the proof of theorem 1will be postponed until we have approximation tools available. 7. The labeled Gromov–Hausdorff distance Our proof of the basic properties of the labeled Gromov–Hausdorff distance is similar to well known proofs of the usual Gromov–Hausdorff distance (see for instance [20, chapter 3] or [6, section 7.3]) but with the labels interwoven into it. For the sake of completeness we record our proof in its entirety, although the aspects unrelated to labels are indeed well known. Proof of proposition 5.The symmetry of the labeled Gromov–Hausdorff distance is evident. The triangle inequality can be verified as in [6, proposition 7.3.16]. Also, if there is has described, then choosing Z=Y,f=h, and g=id shows that the distance is zero. The only nontrivial claim is that if the distance is zero, then such an hexists. Let us prove that. Take any compact metric space Zand isometric embeddings f:X→Zand g:Y→Z. We can define a semimetric (all properties of a metric but d(x,y) = 0 need not imply x=y)don the disjoint union W=XtYby letting d(a,b) =          dZ(f(a),f(b)),when a∈Xand b∈X, dZ(f(a),g(b)),when a∈Xand b∈Y, dZ(g(a),f(b)),when a∈Yand b∈X, dZ(g(a),g(b)),when a∈Yand b∈Y. The embeddings are isometric, so dZ(f(a),f(b)) = dX(a,b)and similarly for Yand g. Take any k∈N. By assumption we have a metric space Zkand thus a semimetric dkon W so that dZk H(X,Y)<1 kand dk(α(ℓ),β(ℓ)) <1 kfor all ℓ∈L. These semimetrics agree with the metrics on Xand Ywhen restricted to either set. By compactness there is a finite set A′ k⊂Xso that the semiballs p∈X;dk(x,p)<1 kx∈A′ k 21 Inverse Problems 39 (2023) 095002 M V de Hoop et al cover X. Similarly, there is a finite A′ ′ k⊂Xso that the semiballs p∈Y;dk(x,p)<1 kx∈A′′ k cover Y. Similarly, there are finite sets B′ kand B′′ kin Yso that the semiballs of radius 1 kcover Y and X. We define Ak=A′ k∪A′′ kand Bk=B′ k∪B′′ k. We can make all of these choices so that Ak⊂Ak+1and Bk⊂Bk+1. For any k∈Non the finite set AktBkthe sequence of semimetrics (dj)∞ j=kis bounded pointwise, since dj|Ak×Ak=dj′|Ak×Akand dj|Bk×Bk=dj′|Bk×Bkfor j,j′⩾k. For x∈Akand y∈Bk the assumption dZk H(X,Y)<1 kyields dj(x,y)⩽max(diam(X),diam(Y)) + 2 k for all j⩾k. Thus (dj)∞ j=khas a converging subsequence and the pointwise limit of semimetrics is a semimetric on AktBk. Constructing a diagonal sequence gives us a subsequence of (dj)which converges pointwise on W′=Sk∈N(AktBk)to a semimetric δ. Because each dk agrees with the original metrics on Xand Y, the only thing to inspect are the ‘cross-distances’ between Xand Y. We can extend δto all of Was follows. When xand yare both in Xor both in Y, we use the metrics on these spaces. When x∈Xand y∈Y, we pick for each k∈Npoints xk∈ BdX(x,1 k)∩Ak⊂Xand yk∈BdY(y,1 k)∩Bk⊂Y, where we have indicated the metric defining the balls as a superscript. We then let δ(x,y) = lim k→∞ δ(xk,yk). A simple argument shows that this limit is independent of the choice of the approximating sequences. It is also straightforward to check that the extended δis indeed a semimetric on W. We want to show that for every x∈Xthere is a unique y∈Yso that δ(x,y) = 0. The triangle inequality and δagreeing with the metric on Yshows that the point is unique. For existence, there is a point in xk∈Akthat is 1 k-close to x, and there is yk∈Bkthat is 1 k-close to xk. We can then set y=limk→∞ yk. Similarly, each y∈Yhas a unique x∈Xso that δ(x,y) = 0. This gives rise to a bijection h:X→Ythat satisfies δ(x,h(x)) = 0. This can be checked to be an isometry. Take any ℓ∈Land k∈N. We have dk(α(ℓ),β(ℓ)) <1 kby construction of dk, so in the limit of the subsequence that we got we find δ(α(ℓ),β(ℓ)) = 0. This means that h(α(ℓ)) = β(ℓ). Therefore h◦α=β. We record two additional propositions. The proofs are immediate and we omit them. Proposition 15. Let X and Y be compact metric spaces and L any set. Let α:L→X and β:L→Y be any two functions. If L′⊂L, then dL′ GH(X,α|L′;Y,β|L′)⩽dL GH(X,α;Y,β). As d∅ GH(X,∅;Y,∅)is the usual Gromov–Hausdorff distance between X and Y, the labeled kind of convergence implies the usual kind of convergence. Proposition 16. Let X be a compact metric space and Y ⊂X a subset. Let α:L→X and β:L→Y be any two functions on a set L. If Y ⊂X is ε1-dense and supℓ∈L|d(α(ℓ),β(ℓ))|⩽ε2, then dL GH(X,α;Y,β)⩽ε1+ε2. 22 Inverse Problems 39 (2023) 095002 M V de Hoop et al 8. Bounded geometry We begin our study of bounded geometry by proving proposition 7with the help of two lemmas, and then we move on to finding further estimates based on these basic bounds. Lemma 17. Let M be a simple manifold, let Ux⊂TxM be the maximal domain of definition of the exponential map, and let U ⊂TM be the subset with fibers Ux. The map exp ×π:U→ M×M that maps U3(x,v)7→(expx(v),x)∈M×M is a diffeomorphism and so has a smooth inverse θ:M×M→U⊂TM. Proof. As each expx:Ux→Mis a diffeomorphism on a simple manifold, the map exp ×π is clearly smooth and bijective. What remains to check is the invertibility of the differential at every point (x,v)∈U. The differential has a convenient block structure due to π(x,v) = x being independent of v, and so d(exp ×π)is invertible at (x,v) if and only if dexpxis invertible at v. The invertibility of dexpxis true by assumption. Lemma 18. Let W1,...,Wnbe vector fields on a simple manifold M constituting a global orthonormal frame. Let Γi jk(x,y)be the Christoffel symbol at x ∈M written in the normal coordinates centered at y ∈M. Let W1(y),...,Wn(y)be the coordinate basis for TyM. Then sup x,y∈MΓi jk(x,y)<∞ for all indices i,j,k. There is a global orthonormal frame on every simple manifold given by vector fields Wi as stated. Such a frame can be produced by the Gram–Schmidt method from a general frame coming from the trivializability of the tangent bundle. Proof of lemma 18.The vectors W1(y),...,Wn(y)constitute an orthonormal basis for TyM. With the help of lemma 17 we see that the basis vectors for TxMin the normal coordinates about yare given by wi(x,y) := d[x]expy(θ(x,y))Wi(y). Here and later in this proof we indicate the variable of differentiation by superscript in [square brackets] when needed for clarity. The invariant coordinate map for the normal coordinates of yis θy:M→TyMgiven by θy(x) = θ(x,y). Using the basis given by the frame, we get a proper coordinate map θy:M→ Rngiven by θy i=Wi(y)♭◦θy, which means θy(x)=(hW1(y),θy(x)i,...,hWn(y),θy(x)i). Now we have a concrete description of the normal coordinates in terms of the frame. Let us denote wi(x,y) = d[x]θy i(x)∈T∗ xM. This is the differential of a coordinate map and thus a basis covector for T∗ xMinduced by the normal coordinates about y. The corresponding basis vectors on TxMare wi(x,y). In terms of the coordinate map z=θy:M→Rnthe Christoffel symbols are given by Γi jk(x,y) = dzi∇∂/∂zj ∂ ∂zk =wi(x,y)(∇[x] wj(x,y)wk(x,y)). 23 Inverse Problems 39 (2023) 095002 M V de Hoop et al (1) Take any p ∈P. If x ∈c(p), then d(p,x)⩽E(p,x). (2) For any x ∈∂M there is p ∈P with d(p,x)⩽E(x) + ε1. (3) For any x ∈∂M there is p ∈P so that d(p,x)⩽E+ε1. Proof. Part (1): take any p∈int(M). Let ρ:M→Rbe the distance function ρ(x) = d(p,x)and r=ρ|∂M. The Hessians Hρand Hrof these functions are symmetric quadratic forms on TM and T∂M, respectively. Let γbe any unit speed geodesic with γ(0) = pand Ja Jacobi field along γso that J(0) = 0. Suppose γmeets the boundary orthogonally, which is equivalent with the exit point being a critical point of r. Then the Hessian of the distance function has the property (cf proof of lemma 13) DtJ(t) = HρJ(t) for all t>0 and thus by part (4) of definition 6we also have hJ(t),HρJ(t)i⩽CJFt−1|J(t)|2. As there are no conjugate points, this amounts to Hρ⩽CJFt−1g in the sense of quadratic forms on ˙γ⊥⊂Tγ(t)M. If λ(x)denotes the smallest eigenvalue of Hr(x)at a boundary point x∈∂M, then (14) yields λ(x)h1(x)⩽Hr(x) = Hρ(x) + h2(x)⩽[CJFr(x)−1+CSFF]h1(x) as quadratic forms on the boundary. If λ(x)>CSFF, this estimate gives r(x)⩽CJF λ(x)−CSFF =: E(p,x). If λ(x)⩽CSFF, then E(p,x) = ∞. The claimed estimate thus holds in both cases. Part (2): let x∈∂M. By the triangle inequality, the estimate d(y,x)⩽d∂M(y,x)and Part (1) we have d(p,x)⩽d(p,y) + d∂M(y,x)⩽E(p,y) + d∂M(y,x) for all p∈Pand y∈c(p). Thus by the definition of E(x) in (4) there is a point p∈Pwith a distance to xless than E(x) + ε1to x. Part (3): follows immediately from the previous one. With the aid of lemma 26 we can now prove proposition 8. Proof of proposition 8.Proposition 3implies that the data determines the arrival time functions asfor all s∈S, although we do not have a description of the index set Syet. As these functions differ from the boundary distance functions rs:∂M→R,rs(x) = d(x,π(s)), only by a constant, we thus know the differential and the Hessian of each rson all of ∂M. Therefore the data determines the critical points of these functions and of the function E(p,y), for p=π(s) and y∈c(p), defined in (3). From this one can easily compute E(x) and Efrom (4) and (5). The set Γis defined in terms of these quantities, so it is uniquely determined by the data as well. Theorem 2indicates that the data determines the pointwise spatial distances of the source points, and the last part of the claim follows. 30 Inverse Problems 39 (2023) 095002 M V de Hoop et al 9.3. Interior density estimates For the density of sources deep in the manifold M, we use two estimates. The first one concerns the density of the set of geodesics connecting near-boundary source points in Γ. The second one ensures that the lentils cover enough of the manifold. Lemma 27. Let M be a simple manifold of bounded geometry and let ε2>0. Let Γ⊂M be such that ∂M⊂B(Γ,ε2). Then for every z ∈int(M)there are points x,y∈Γso that d(z,γx,y([0,d(x,y)])) ⩽Cdε2. Proof. Fix any x∈Γand take any z∈int(M). Let ˆγ=γx,zbe a constant speed geodesic with ˆγ(0) = xand ˆγ(t) = zfor some t>0. We extend this geodesic beyond zso that it meets ∂Mat some time t′>t. We scale the constant speed so that t′=1 and we denote ˆ y:= ˆγ(1)∈∂M. There is y∈Γso that d(y,ˆ y)< ε2. Let γ: [0,1]→Mbe the constant speed geodesic for which γ(0) = xand γ(1) = y. By lemma 23 we have d(ˆγ(t),γ(t)) ⩽Cdtd(ˆ y,y). As t<1, we have thus d(z,γ([0,1])) ⩽d(ˆγ(t),γ(t)) <Cdε2 as claimed. Lemma 28. Take any ε1>0and ε2>0, and let M be a simple manifold with bounded geometry and Γ⊂B(∂M,ε1)so that ∂M⊂B(Γ,ε2). Let z ∈M. Whenever d(z,∂M)⩾ε1+Cgε2, then there is a lentil Lx,y r,swith x,y∈Γand thickness δ=C1ε2<min(r,s)containing z. The constants are given by (40) and (51) in appendix B. Proof. We showed in lemma 27 above that for any z∈int(M)there are x,y∈Γso that the distance from zto the trace of the geodesic γx,yis at most Cdε2. We will choose h>0 later (see (24) for an explicit expression) and require that d(x,z)⩾hand d(y,z)⩾h. These additional requirements may not hold for the pair of points x,y∈Γwithout further assumptions for z. As both xand yare ε1-close to the boundary, the additional requirements d(x,z)⩾h and d(y,z)⩾hare certainly satisfied if d(z,∂M)⩾ε1+h. This is why we placed a boundary distance assumption on zin the statement of this lemma. Let z′∈γx,ybe the nearest point on this geodesic to the point z. After choosing hlarge enough we obtain by the reverse triangle inequality d(x,z′),d(y,z′)⩾h−Cdε2>0.(21) This implies that z′is an interior point of γx,y, and the geodesics γx,yand γz,z′meet orthogonally at z′. To make sure that zis contained in a lentil Lx,y r,s, with thickness δand midpoint z′, we require that the transversal radius of that lentil is more than d(z,z′). By lemma 25 the aforementioned condition of the transversal radius is satisfied when Cdε2⩽1 2δ. (22) The condition for z′to be a midpoint of a lentil of thickness δis that d(x,z′)>1 2δand d(y,z′)> 1 2δ, which by equation (21) are satisfied when h−Cdε2>1 2δ. (23) 31 Inverse Problems 39 (2023) 095002 M V de Hoop et al Combining the two conditions (22) and (23) gives Cdε2<h−Cdε2 and therefore h>2Cdε2. We choose h=3Cdε2.(24) The conditions (22) and (23) then become 2Cdε2⩽δ < 4Cdε2, so we may choose δ=2Cdε2.(25) This is why we chose C1in (40) as we did. Finally we note that by the choice of z′∈γx,yand due to (21), (24), and (25) we have d(x,y)> δ. If we set r=d(x,z′) + 1 2δand s=δ+d(x,y)−r=d(z′,y) + 1 2δ, then the lentil Lx,y r,sis of width δand has z′as a midpoint. Moreover, by the previous argument this lentil contains the point z, if d(z,∂M)⩾ε1+Cgε2, as required in the claim of this lemma. The choices of hand δin the proof above are optimal up to a constant and the explicit choice simplifies our estimates. 9.4. Global density estimates With all the estimates we have collected, we are ready to prove both parts of theorem 9. Claim (1) states that if all lentils of thickness δmeet a source point (which is verifiable from data), then the source set P1⊂M1has a concrete density estimate. The thickness δis chosen carefully depending on ε2. Claim (2) states that under the assumptions of the previous claim the measurements define a discrete metric space which is a quantitatively good approximation for the true space. Proof of theorem 9.For claim (1), take any z∈M. We wish to show that there is p∈Pso that d(z,p)< ε with εchosen as in (8). We split the proof in two cases: near the boundary and deep within the manifold, and use different tools in either of these cases to verify the validity of the density claim. The boundary case: if d(z,∂M)< ε1+Cgε2, then there is x∈∂Mso that d(z,x)< ε1+ Cgε2. By claim (3) of lemma 26, there is p∈Pso that d(x,p)<E+ε1. Recall that ε2=E+ ε1> ε1. By the triangle inequality we have d(z,P)⩽d(z,p)⩽d(z,x) + d(x,p)<(Cg+2)ε2.(26) This is a sufficient estimate near the boundary. The deep interior case: by definition of the set Γin (6) and lemma 26 we have Γ⊂B(∂M,ε1) and ∂M⊂B(Γ,ε2). If d(z,∂M)⩾ε1+Cgε2, then lemma 28 shows that there is a lentil L=Lx,y r,s of the correct thickness δ=C1ε2<r,sso that z∈L. By the crucial assumption of the theorem, there is p∈P∩L. Lemma 24 then gives us 32 Inverse Problems 39 (2023) 095002 M V de Hoop et al d(z,P)⩽d(z,p)⩽δ+Ce√δ=C1ε2+CepC1ε2.(27) This is a sufficient estimate for deep interior points. Now it remains to combine the estimates for points in the deep interior and near the boundary. Combining (26) and (27) gives d(z,P)⩽max(Cg+2)ε2,C1ε2+CepC1ε2 ⩽(Cg+2+C1)ε2+CepC1ε2 for all z∈M. This proves that P⊂Mis ε-dense with the choice of (8) and concludes the proof of claim (1). The labeled Gromov–Hausdorff distance is now straightforward to estimate for claim (2).By the definition of E(x) in (4) we can choose a function α:∂M→Pthat satisfies (9); heuristically, the point α(x)is an ε1-approximate minimizer of the right-hand side of (4). This function can be actually identified from the data Q(S) as we know the first fundamental form of ∂Mand the value of E(p,y)for all p∈π(S)and y∈c(p), given by (3), is obtained by the computing the gradient and the Hessian of the corresponding arrival time function.Then by lemma 26 it follows that for each x∈∂Mthere exists y∈c(α(x)) that satisfies ε2⩾E(x) + ε1>E(α(x),y) + d∂M(x,y) ⩾d(α(x),y) + d(x,y) ⩾d(α(x),x). Since P⊂Mis ε-dense, proposition 16 finally gives d∂M GH(P,α;M,ι)⩽ε+ε2. By the choice of εin (8) this proves estimate (10) and completes the proof. 9.5. Reverse density estimates We do not need any more preparations before the next proof. Proof of theorem 10.Consider a unit speed geodesic γstarting at a source point p=γ(0)∈ P. Take any t>0 and any w∈Tγ(t)Morthogonal to ˙γ(t). By bounded geometry we have hw,Hρwi⩾|w|2(CH1t−1−CH2). Now suppose that y=γ(t)is on ∂Mand ˙γ(t)is normal to the boundary. We recall the earlier notation that ρ:M→Ris the distance function to the point pand r=ρ|∂M. The second fundamental form h2is positive definite by simplicity and by (14) we deduce that the Hessians on ∂Mand Mevaluated at yare related by Hr(y) = Hρ(y) + h2(y)⩾(CH1t−1−CH2)h1(y), where h1is the first fundamental form. Therefore the smallest eigenvalue λ(p,y)of Hr(y) satisfies λ(p,y)⩾CH1d(p,y)−1−CH2.(28) This enables us to estimate E(p,y). 33 Inverse Problems 39 (2023) 095002 M V de Hoop et al To utilize the assumption that P⊂Mis ˆε-dense with ˆεchosen by (11), suppose that d(∂M,p)<ˆεand let zp∈∂Mbe the nearest boundary point to it. Estimate (28) then implies that E(p,zp)⩽CJF CH1 ˆε−1−CH2 −CSFF when ˆε < CH1/(CH2 +CSFF). Due to our choice of ˆε, we have ˆε⩽CH1 2(CH2 +CSFF),(29) and so E(p,zp)⩽2CJFC−1 H1 ˆε. (30) Consider then an arbitrary point x∈∂M. There is a source point q∈Pwith d(x,q)<ˆε. Let zq∈∂Mbe the closest boundary point to q. We haved(zq,q)<ˆε, and by bounded geometry (part (6) of definition 6) d∂M(x,zq)⩽CdistdM(x,zq) ⩽Cdist(dM(x,q) + dM(zq,q)) <2Cdist ˆε. By (30) we get E(x)⩽E(q,zq) + d∂M(x,zq)⩽2CJFC−1 H1 ˆε+2Cdist ˆε. As this bound it is independent of x∈∂M, we have the same bound for Efrom (5): E⩽2(CJFC−1 H1 +Cdist)ˆε, (31) under the assumption (29). Combining our choices ε1=C5ˆεand ε2=ε1+Ewith (31), we get ε2=ε1+E⩽(C5+2CJFC−1 H1 +2Cdist)ˆε. (32) All other estimates for ε2follow from this one. To satisfy (12), we require that C4ε2<1 2εand C3√ε2<1 2ε. By (32) both follow from requiring C4(C5+2CJFC−1 H1 +2Cdist)ˆε < 1 2ε and C3q(C5+2CJFC−1 H1 +2Cdist)ˆε < 1 2ε, or equivalently ˆε < min ε 2C4(C5+2CJFC−1 H1 +2Cdist),ε2 4C2 3(C5+2CJFC−1 H1 +2Cdist)!.(33) Thus (29) and (33) are satisfies due to (11), and by the previous remark the estimate (12) follows. It remains to show that each lentil L=Lx,y r,d(x,y)−r+δwith x,y∈Pcontain source points as claimed. The thickness of the lentil Lis δ=C1ε2, so by lemma 25 B(m,1 2δ)⊂L, where mis 34 Inverse Problems 39 (2023) 095002 M V de Hoop et al the midpoint of the lentil L. Therefore each lentil of thickness δcontains a ball of radius 1 2δ. It follows from ˆε-density of P⊂Mthat that each such lentil contains a source point because ˆε=C−1 5ε1=1 2C1ε1⩽1 2C1ε2=1 2δ. This is why we chose C5so that C−1 5=1 2C1. 10. Convergence of discrete approximations 10.1. Deterministic convergence With proposition 8, theorems 9and 10 the proof of theorem 11 is straightforward. Proof of theorem 11.We begin with proposition 3with Ω = ∂M×(0,T). Out of the parts of graphs we only choose the ones that are full graphs of a function as:∂M→R. It follows from the diameter bound that τ(s)⩽as(x)⩽τ(s) + Cdiam for all x∈∂M. Therefore the set Q(S,T)of (13) contains the full graph of asfor all the sources s∈Swith 0 < τ(s)<T−Cdiam. Let us denote the set of source points with their complete graphs contained in ∂M×(0,T)by PT⊂π(S). If Tis too small, we may have PT=∅. In this case we choose the metric approximation MT to be a set of one point and αT:∂M→MTthe constant map. When there are sources, we set MT=PT. Take any ε > 0 and let ˆεbe given by (11). As ST>0PTis dense in Mby assumption, by compactness there is T(ε)>0 so that PTis ˆε-dense for all T⩾T(ε). With the choices ε1=C5ˆε and ε2=ε1+E(T), where E(T) is defined as in (5), and δ=C1ε2we have C4ε2+C3√ε2< ε, (34) by theorem 10. By the last claim of theorem 10 the assumption of theorem 9is satisfied. Thus the set PTis also ε-dense as required in claim ((2)) of theorem 9. We choose a map αT:∂M→ PTso that (9) is satisfied when α=αTand P=PT. Finally the estimate (10) with (34) yields d∂M GH(PT,αT;M,ι)< ε. It should be noted that due to proposition 8, the auxiliary quantity E(T) and the map αTused above are defined by the data. 10.2. Dense sources Now we are finally able to prove theorem 1. We will make use of a version of the Myers– Steenrod Theorem [46,47] which states that a metric isometry between smooth Riemannian manifolds is necessarily smooth. In the case of simple manifolds it is straightforward to prove and the boundary causes no technical trouble. We record the proof here for the sake of completeness. Lemma 29 (Myers–Steenrod Theorem). If Φ: M1→M2is a metric isometry between simple Riemannian manifolds, it is smooth up to the boundary and an isometry also in the sense that Φ∗g2=g1. 35 Inverse Problems 39 (2023) 095002 M V de Hoop et al Proof. It suffices to show that Φis smooth. Then preserving distances quickly implies Φ∗g2= g1via differentiating distance functions. Take any x1∈int(M1)and let x2= Φ(x1). If γis any constant speed geodesic with γ(0) = x1, then Φ◦γis a constant speed geodesic starting at x2with the same speed. Let Φ∗:Tx1M1→ Tx2M2be the map that maps the velocity vectors of these geodesics to each other. That is, if ∂tγ(t)|t=0=v, then ∂tΦ(γ(t))|t=0= Φ∗(v). We have Φ◦expx1=expx2◦Φ∗. The map Φ∗is clearly homogeneous in positive scalings, maps unit vectors to unit vectors, and is bijective. By considering reversed geodesics we find Φ∗(−v) = −Φ∗(v)for all v∈Tx1M1. For any two v,w∈Tx1M1we have [28, equation (2.15)] d(expx1(tv),expx1(tw))2=t2(|v|2+|w|2−2hv,wi) + O(t4) for small t>0. Due to Φ∗preserving the norm, the intertwining property Φ◦expx1=expx2◦ Φ∗, and the isometric nature of Φ, this implies that hv,wi=hΦ∗v,Φ∗wi for all v,w∈Tx1M1. The map Φ∗preserves norms and inner products, and therefore it preserves the (squared) distance between any pair of vectors in Tx1M1. Thus Φ∗is an isometry between finitedimensional inner product spaces, and as a Euclidean isometry fixing the origin it is linear and thus smooth. As the exponential maps are diffeomorphisms on their maximal domain of definition, the map Φ = expx2◦Φ∗◦exp−1 x1is smooth. Proof of theorem 1.Let ι1:∂M1→M1be the usual inclusion map and let ι2=ϕ:∂M1→ M2(with extended codomain here for convenience) so that we have boundary inclusions ιi:∂M1→Miwith the same domain. For any T>0 consider data on the bounded set ∂M×(−T,T). This time the time interval extends into both the past and the future. As in the proof of theorem 11 above, we get a finite metric space MTand a map αT:∂M1→MTso that d∂M1 GH (MT,αT;Mi,ιi)→0 as T→∞ for both i=1,2. By proposition 8the approximating metric space MTconstructed from data is the same for the two manifolds M1and M2. We then apply proposition 5. By the triangle inequality d∂M1 GH (M1,ι1;M2,ι2)→0 as T→∞and so d∂M1 GH (M1,ι1;M2,ι2) = 0. Now proposition 5gives an isometry Φ: M1→M2 so that Φ◦ι1=ι2. The condition Φ◦ι1=ι2simply means that Φ|∂M1=ϕ. By lemma 29 the isometry Φis actually smooth. It remains to check that the sources s= (p,t)∈S1⊂int(M1)×Rare mapped correctly. By theorem 2we know that there is a bijection, but we have to verify that it corresponds to the isometry Φas claimed. In light of proposition 3, the data can be seen as a collection of graphs. Let us denote Ai={as;s∈Si}. The two manifolds having equivalent data means that the map ϕ∗:A2→A1that takes a7→ϕ∗a=a◦ϕis bijective. We will show that the maps bi:Si→Ai with bi(s) = asare bijective, and therefore the natural bijection between the sources S1→S2 is ξ=b−1 2◦(ϕ∗)−1◦b1. We need this map to satisfy (2), which now amounts to a(p,t)=a(Φ(p),t)◦ϕ 36 Inverse Problems 39 (2023) 095002 M V de Hoop et al as functions ∂M1→Rfor all (p,t)∈S1. This follows straightforwardly from the definitions and Φbeing an isometry. Let us then show that biis bijective for both i. To this end, take any two distinct sources s,ˆ s∈Si. The two arrival time functions are as=rs+τ(s)and aˆ s=rˆ s+τ(ˆ s). If the spatial source points are different, π(s)6=π(ˆ s), then rs−rˆ scannot be a constant function on ∂Mdue to lemma 14. Therefore asand aˆ sdo not coincide. If the spatial points are the same, then rs−rˆ sis the constant function 0. As s6=ˆ sbut π(s) = π(ˆ s), we must have τ(s)6=τ(ˆ s)and thus as−aˆ sis a non-zero constant function. We have thus proven that s6=ˆ simplies bi(s)6=bi(ˆ s). This concludes the proof that biis injective and thus the theorem is proven. 10.3. Stochastic convergence The proof of the stochastic result only requires checking that the point process almost surely produces the correct kind of source set. Proof of proposition 12.Countability of the set S⊂M×Rgiven by the homogeneous Poisson point process is certain. We need to prove that the following properties hold almost surely: (1) π(S)∩∂M=∅. (2) S⊂M×Ris discrete. (3) π(S∩(M×[0,∞))) ⊂Mis dense. The probability that a measurable set A⊂M×Rcontains k∈Npoints of the source set Sis P[#(A∩S) = k] = e−λµ(A)(λµ(A))k k!.(35) For this and other basic properties of homogeneous Poisson point processes, see e.g. [23,57]. The first property follows simply from µ(∂M×R) = 0. The second property follows from A∩Sbeing almost certainly finite when µ(A)<∞. For the third property, let (xi)∞ i=1be a dense sequence in Mand define for any i,j⩾1 the set Ai,j=BM(xi,j−1)×[0,∞)⊂M×R. Comparability with the natural product measure gives µ(Ai,j) = ∞. Therefore it follows from (35) that P[#(Ai,j∩S) = 0] = 0. The projection π(S∩(M×[0,∞))) is dense in Mif #(Ai,j∩S)>0 for all iand j. As this event for individual indices has probability 1, the event for density is a countable intersection of events of full probability. Therefore density is almost certain. The conditions of theorem 11 are thus met almost surely. Data availability statement No new data were created or analysed in this study. 37 Inverse Problems 39 (2023) 095002 M V de Hoop et al Acknowledgments M V d H was supported by the Simons Foundation under the MATH +X program, the National Science Foundation under Grant DMS-1815143, and the corporate members of the Geo-Mathematical Imaging Group at Rice University. J I was supported by the Academy of Finland (Projects 332890 and 336254). M L was supported by Academy of Finland (Projects 284715 and 303754). T S was supported by the Simons Foundation under the MATH +X program and the corporate members of the Geo-Mathematical Imaging Group at Rice University. The authors want to thank Peter Caday and Vitaly Katsnelson for useful discussions and the anonymous referees for valuable comments. Appendix A. An improved estimate on transversal radius As mentioned in connection to lemma 25, the estimate can be improved although it is not necessary for the proof of our theorems. The estimate R≳δof lemma 25 can be improved to R≳√δas follows. Proposition 30. Consider a simple manifold of bounded geometry with the constants satisfying the condition (1). Take any x,y∈M and consider the lentil Lx,y r,swith r,s∈(0,d(x,y)). Suppose δx,y r,s=r+s−d(x,y)satisfies δx,y r,s⩽min(r,s). The transversal radius of any lentil satisfies Rx,y r,s>minChqδx,y r,smin(r,s),1 2Cdiam, where the constant is given by (52) in appendix B. Most of our constants are used for estimates from above, but Chis used for estimating from below. Therefore, unlike most of our constants, it can be made smaller but not larger if needed. It is also natural that the transversal radius cannot exceed half of the diameter. Lemma 25 is more convenient to use and we do not benefit significantly from the improved exponent, so we will not employ proposition 30 in the proofs of our main results. Proof of proposition 30.Consider a lentil L=Lx,y r,swith thickness δ=δx,y r,sand any point z∈ Mfor which the closest point on γx,yis the midpoint m=mx,y r,s. We want to find conditions on the distance w:= d(z,m)which ensure that z∈L. We will do so by making explicit estimates on a comparison manifold of constant sectional curvature and then translating the resulting estimate to the actual manifold. We assume that w⩽1 2Cdiam. Due to (1) this implies w2Csec+< π2/4, and so 1 − 1 2w2Csec+∈(−1,1]. This will ensure that the argument of arccos: [−1,1]→[0,π]stays within the domain in the following treatment. We may freely assume w>0 when convenient, as the case w=0 can be given a trivial separate treatment. Let us denote R:= d(x,m)and d:= d(x,z). The points x,z, and mform a triangle with a right angle at m. Consider the corresponding triangle on the MCsec+with constant sectional curvature Csec+with with the right angle at ˜ mwith the same lengths Rand wof the catheti. The length of the hypotenuse is ˜ d=C−1/2 sec+arccos[cos(RpCsec+)cos(wpCsec+)]. 38 Inverse Problems 39 (2023) 095002 M V de Hoop et al The cosine satisfies cos(x)⩾1−1 2x2, so ˜ d⩽C−1/2 sec+arccos[cos(RpCsec+)[1−1 2w2Csec+]]. The function t7→arccos(1−t2)is Lipschitz-continuous on [0,T]with the Lipschitz constant 2 √2−T2 whenever T<√2. In our setting T=√Cc(with the constant of (39)) and the two values of t where we compare arccos(1−t2)are q1−cos(RpCsec+)[1−1 2w2Csec+] and q1−cos(RpCsec+). The condition (1) implies that indeed T<√2, so the Lipschitz constant is Cias given in (53). This Lipschitz-continuity now gives ˜ d⩽R+C−1/2 sec+Ciq1−cos(RpCsec+)[1−1 2w2Csec+] −q1−cos(RpCsec+). Using the estimate √1−a+ab ⩽√1−a+1 2|ab|√1−a(which is valid for a<1 and b∈R) leads to ˜ d⩽R+C−1/2 sec+Ci Csec+cos(RpCsec+) 4q1−cos(RpCsec+) w2. A simple calculation gives |cos(x)|[1−cos(x)]−1/2⩽4/xfor all x∈(0,π). Thus we find ˜ d⩽R+CiR−1w2, completing our estimate of the comparison length ˜ d. By lemma 20 we have d⩽˜ d, and so d⩽R+CiR−1w2.(36) Let us use this to see when z∈L. To ensure z∈L, we need d<r. As r=R+1 2δ, estimate (36) says that d<rwhenever CiR−1w2<1 2δ. Thus we get the condition w2<(2Ci)−1δ(r−1 2δ)and the rougher bound 1 2Cdiam. A similar treatment of the condition d(y,z)<sleads us to w2<(2Ci)−1δ(s−1 2δ). We have thus shown that Rx,y r,s>min(2Ci)−1/2qδ(r−1 2δ),(2Ci)−1/2qδ(s−1 2δ),1 2Cdiam. As δ < r,s, this implies the claimed inequality. 39