Lucie Baudouin
email: baudouin@laas.fr

Maya De Buhan
email: maya.de-buhan@parisdescartes.fr

Sylvain Ervedoza
email: ervedoza@math.u-bordeaux.fr

Axel Osses
email: axosses@dim.uchile.cl

Carleman-based reconstruction algorithm for the waves

Keywords: hyperbolic equation, inverse problem, reconstruction algorithm, Carleman estimates AMS subject classifications: 93B07, 35R30

1 Introduction and main resuts.

Setting

We consider a smooth bounded domain of R d (d 1) denoted by Ω and a time horizon T > 0. This article focuses on the reconstruction of the velocity in a wave equation, that is for the following inverse problem:

Given the distributed and boundary source terms f and f ∂ and the initial data (w 0 , w 1 ), considering the solution of the hyperbolic equation

   ∂ 2 t w * -∇ • (a * ∇w * ) = f, in (0, T ) × Ω, w * = f ∂ ,
on (0, T ) × ∂Ω, w * (0) = w 0 , ∂ t w * (0) = w 1 , in Ω, (1.1) can we determine the unknown coefficient a * = a * (x), assumed to depend on the spacevariable x only, from the additional knowledge of the flux of the solution through a part Γ 0 of the boundary ∂Ω, namely

a * ∂ ν w * = a * ∇w * • ν, on (0, T ) × Γ 0 , (1.2)
where ν ∈ R d is the unit outward normal vector to ∂Ω?

Before going any further, let us recall that if a ∈ W 1,∞ (Ω) with inf Ω a α 0 > 0, f ∈ L 1 (0, T ; L 2 (Ω)), f ∂ ∈ H 1 ((0, T ) × ∂Ω), w 0 ∈ H 1 (Ω) and w 1 ∈ L 2 (Ω) and assuming the compatibility condition f ∂ (0, x) = w 0 (x) for all x ∈ ∂Ω, the initial and boundary value problem

   ∂ 2 t w -∇ • (a∇w) = f, in (0, T ) × Ω, w = f ∂ ,
on (0, T ) × ∂Ω, w(0) = w 0 , ∂ t w(0) = w 1 , in Ω, (1.3) is well-posed, and has a unique solution

w ∈ C 0 ([0, T ]; H 1 (Ω)) ∩ C 1 ([0, T ]; L 2 (Ω)),
whose normal derivative a∂ ν w is well-defined as an element of L 2 ((0, T )×∂Ω), see [START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF]Theorem 4.1].

We emphasize that we are interested in the effective reconstruction of the coefficient a * (x) in (1.1) from the extra information given by the outgoing flux. As a matter of fact, this coefficient a * corresponds to the speed of propagation of waves up to a square root (if c * denotes the speed of propagation of waves in the medium, then a * (x) = c 2 * (x)), and we will sometimes call it the speed coefficient with a slight abuse of denomination. In particular, we will simply recall the results obtained concerning uniqueness and stability for this inverse problem. These questions are already well documented, for instance the Lipschitz stability from large geometrical observations has been studied in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF][START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF], among others.

The approach we propose to reconstruct the coefficient a * , inspired by our previous works [START_REF] Baudouin | Global carleman estimates for waves and applications[END_REF][START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF] for the reconstruction of an unknown potential in a wave equation through a measurement of the flux, is based on an algorithm that involves Carleman estimates, and whose convergence will be ensured when taking a tuning parameter, corresponding to the Carleman parameter, large enough. The goal of this article is thus to fully develop this strategy, from the theoretical aspects to the numerical ones.

The paper is organized as follows. The end of this introductory section describes the stability results already known for the inverse problem under consideration and will present the reconstruction algorithm of the speed coefficient built to be implemented numerically. Section 2 presents the proof of the Lipschitz stability in the recovery of the speed from the flux, since this will enlighten the algorithms' design. The proof of the convergence of the algorithm will be detailed in Section 3. Section 4 presents the numerical implementation, Section 5 then provides the reader with some further comments, and the appendix gives the proofs of the Carleman estimates we shall use.

Stability results

The results we will introduce require the following quite usual geometric conditions (sometimes called "multiplier condition" or "Γ-condition"):

Ω Γ 0 L ν(x) x x 0
Geometric condition:

∃ x 0 ∈ Ω, such that Γ 0 ⊃ {x ∈ ∂Ω, (x -x 0 ) • ν(x) 0}, (1.4)
The geometric condition (1.4) appears naturally when deriving observability results for the wave equation (1.3) using the multiplier (x -x 0 ) • ∇w + (d -1)w/2 for constant wave speed in large times, developed in references [START_REF] Ho | Observabilité frontière de l'équation des ondes[END_REF][START_REF] Lions | Contrôlabilité exacte, Stabilisation et Perturbations de Systèmes Distribués. Tome 1. Contrôlabilité exacte[END_REF]. In fact, when considering variable wave speed a ∈ W 1,∞ (Ω), following the multiplier technique, one easily checks the observability of the wave equation (1.3) under a condition of the form ∃ρ ∈]0, 1], such that ∇a • (x -x 0 ) 2(1 -ρ)a in D (Ω).

(

Indeed, this result can be achieved by multiplying the equation (1.3) by (x -x 0 ) • ∇w + (d -1)w/2, and the time T for observability through Γ 0 provided by that approach satisfies ρT > 2L √

inf Ω a , where we have set L = sup Ω {|x -x 0 |}. (1.6) Details are left to the reader (similar arguments can be found for instance in [START_REF] Dehman | Observability estimates for the wave equation with rough coefficients[END_REF]). Since we will develop a strategy based on Carleman estimates, which can be seen as a fancy version of multiplier techniques to some extent, we start by making some assumptions on the speed a * that we aim at recovering linked to the previous discussion. In particular, we suppose that we have some a priori informations on a * and some a priori bounds m > 0, 0 < α 0 < α 1 and ρ ∈]0, 1] such that

a * ∈ C 1 Ω ∩ H 2 (Ω) with ∇a * L ∞ (Ω) m, 0 < α 0 inf Ω a * sup Ω a * α 1 , ∇a * • (x -x 0 ) 2(1 -ρ)a * in Ω, (1.7)
where x 0 is the one in (1.4). We shall also assume that a * is known on the boundary, i. Another set of assumptions concerns the solution w * of (1.1), namely

w * ∈ H 2 (0, T ; W 2,∞ (Ω)), (1.9)
and the initial conditions:

w 0 ∈ W 3,∞ (Ω) and inf Ω |∇w 0 (x) • (x -x 0 )| r 0 > 0, (1.10)
w 1 = 0.

(1.11)

Then we can prove uniqueness and stability of our inverse problem in the following class:

V a,aν ,α0,α1,ρ,m := a ∈ C 1 (Ω) ∩ H 2 (Ω), ∇a L ∞ (Ω) m, 0 < α 0 inf Ω a sup Ω a α 1 and ∇a • (x -x 0 ) 2(1 -ρ)a in Ω, a = a and ∂ ν a = a ν on ∂Ω , (1.12)

for which we will use the shorthand notation V a,aν = V a,aν ,α0,α1,ρ,m .

In the sake of clarity, we will prove the specific following Lipschitz stability result:

Theorem 1.1. [Stability of the inverse problem] Let Γ 0 satisfy (1.4) for some x 0 / ∈ Ω and assume that a * satisfies (1.8) and (1.7) for some m > 0, 0 < α 0 α 1 and ρ ∈ (0, 1], w * satisfies (1.9), (w 0 , w 1 ) satisfy (1.10)- (1.11) and let T be such that

T > L α 0 ρ 0 (Lm/α 0 ) , (1.13)
where ρ 0 = ρ 0 (r) is the function given by

ρ 0 (r) = ρ + r 2 4 - r 2 2
.

(

1.14)

There exists a constant C > 0 such that for all a ∈ V a,aν (defined in (1.12)), denoting by w the solution of equation (1.3), we have

a -a * H 1 0 (Ω) C ∂ ν w -∂ ν w * H 2 (0,T ;L 2 (Γ0)) .
(1.15)

Remark 1.2. According to (1.9), ∂ ν w * belongs to H 2 (0, T ; L 2 (∂Ω)). Besides, as we will see in (2.3), for w solution of (1.3) with a ∈ V a,aν , one can check that ∂ ν (w -w * ) belongs to H 2 (0, T ; L 2 (∂Ω)), and thus ∂ ν w also belongs to H 2 (0, T ; L 2 (∂Ω)) and the right hand-side of (1.15) is always finite.

Remark 1.3. The time condition (1.13) differs from the time condition (1.6). First, it differs from a factor 2 which is due to the fact that for the inverse problem under consideration, we use a time symmetrization argument, and thus observing the flux on (0, T ) is equivalent to observe the flux on (-T, T ), see the proof of Theorem 1.1. Formula (1.13) also differs from the time condition (1.6) by the fact that ρ 0 (Lm/α 0 ) appears instead of ρ. This is probably due only to the method of proof. Also note that ρ 0 (0) = ρ, so that when m = 0, i.e. when a is constant, we recover from (1.13) a time condition of the form T > L/ √ α 0 ρ, which looks very similar to (1.6) except for the square root, and the time L/ √ α 0 ρ is slightly larger than the time L/( √ α 0 ρ) -which should be the expected time in view of (1.6) -for ρ ∈ (0, 1).

Although Theorem 1.1 is not entirely new (it should be compared at least with the results in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF] and [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF] that we shall comment later), its proof is at the origin of the construction of the algorithm that we will propose next. We will therefore give it in Section 2.

More precisely, as in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF][START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF], it will require two Carleman estimates, one for the wave equation, and one for a transport equation. The one for the wave equation reads as follows:

Theorem 1.4. [Carleman estimates for a wave equation] Let x 0 ∈ Ω ⊂ R d and assume condition (1.4) for Γ 0 . Let m > 0, 0 < α 0 α 1 and ρ ∈ (0, 1], and define the admissible set

V = {a ∈ W 1,∞ (Ω), ∇a L ∞ (Ω) m, 0 < α 0 inf Ω a sup Ω a α 1 and ∇a • (x -x 0 ) 2(1 -ρ)a a.e. in Ω}. (1.16)
There exists a parameter ρ 0 = ρ 0 (Lm/α 0 ) > 0 such that for all β ∈ (0, α 0 ρ 0 ), there exist s 0 > 0 and a positive constant C such that for all s s 0 and for all a ∈ V,

Ω e 2sϕ(0) ∂ t v(0) 2 + |∇v(0)| 2 + s 2 v(0) 2 dx + s T -T Ω e 2sϕ (∂ t v) 2 + |∇v| 2 + s 2 v 2 dxdt C T -T Ω e 2sϕ ( a v) 2 dxdt + Cs T -T Γ0 e 2sϕ |∂ ν v| 2 dσdt + Cs {ϕ<0} e 2sϕ (∂ t v) 2 + |∇v| 2 + s 2 v 2 dxdt + Cs Ω e 2sϕ(±T ) ∂ t v(±T ) 2 + |∇v(±T )| 2 + s 2 v(±T ) 2 dx, (1.17) for all v ∈ L 2 ((-T, T ); H 1 0 (Ω)) satisfying a v = ∂ 2 t v -∇ • (a∇v) ∈ L 2 ((-T, T ) × Ω) and ∂ ν v ∈ L 2 ((-T, T ) × ∂Ω)
, where ϕ denotes the weight function

ϕ(t, x) = |x -x 0 | 2 -βt 2 ,
(1.18)

and {ϕ < 0} = {(t, x) ∈ (-T, T ) × Ω, ϕ(t, x) < 0}.
As such, Theorem 1.4 is not new, even if we did not find it as stated here in the literature. In fact, it can be deduced from [22, Theorem 2.2.4 p.45], but we make it slightly more precise in particular on the estimate on ρ 0 . This is why the proof is given in appendix, Section A.1.

We shall also need a Carleman estimate for a transport equation, which can be found for instance in [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF]Lemma 4.1] (or in a weaker form in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF]Proposition 2.2]). We give the proof in appendix, Section A.2 in the sake of completeness and readability. Theorem 1.5. [Carleman estimates for a transport equation] Let x 0 ∈ Ω and X be a vector field such that

X ∈ W 2,∞ (Ω; R d ) ∩ C 0 (Ω; R d ),
and

inf x∈Ω {|X(x) • (x -x 0 )|} > 0,
and set

γ X = sign(X(x) • (x -x 0 )) and Γ X = x ∈ ∂Ω, (X(x) • ν(x))γ X > 0 . (1.19)
(Note that, under the above assumptions, x → X(x) • (x -x 0 ) does not change sign in Ω, so the choice of x ∈ Ω in the definition of γ X is irrelevant.) Then there exist s 0 > 0 and a positive constant C such that for all s s 0 ,

Ω e 2s|x-x0| 2 |∇(∇ • (bX))| 2 + s 2 |∇b| 2 + s 4 b 2 dx C Ω e 2s|x-x0| 2 |∇ (∇ • (bX)) | 2 + s 2 |∇ • (bX)| 2 dx (1.20) for any b ∈ H 1 X (Ω) satisfying ∇ • (bX) ∈ H 1 X (Ω) where H 1 X (Ω) = b ∈ H 1 (Ω), b = 0 on Γ X .
Let us also point out that similar strategies were used in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF][START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF] to deduce stability results for the inverse problem under consideration. On the one hand, Hölder stability results have been obtained for large time and from large geometrical boundary observations in [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF] for the L 2 -norm distance of the coefficients, assuming that the speed is known at the boundary and in [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF], Lipschitz stability results have been obtained for the H 2 -norm distance between the coefficients, assuming that the speed and some of its normal derivatives are known at the boundary, but assuming C 3 a priori bounds on the coefficients.

The result of Theorem 1.1 uses global Carleman estimates, and thus requires geometric conditions which are very close to the multiplier conditions (1.4), (1.13). But one can in fact go further provided that we use a foliation condition as it is explained in [START_REF] Stefanov | Recovery of a source term or a speed with one measurement and applications[END_REF], which states, roughly speaking, that the whole set Ω can be foliated into strictly convex surfaces; in our case, these strictly convex surfaces simply correspond to the level sets of |x -x 0 | 2 . However, this strategy seems hard to implement numerically, since this condition is rather abstract, and we should thus stick to the strategy based on Carleman estimates for the wave equation with weights of the form (1.18).

We also recall that logarithmic stability results are proved in large times for arbitrarily small boundary observations for recovering a potential [START_REF] Bellassoued | Global logarithmic stability in inverse hyperbolic problem by arbitrary boundary observation[END_REF] or the speed [START_REF] Bellassoued | Logarithmic stability in determination of a coefficient in an acoustic equation by arbitrary boundary observation[END_REF], but the corresponding stability results are logarithmic and therefore numerical issues seem more intricate and are postponed to future works.

A Carleman based Reconstruction Algorithm (C-bRec)

We now focus on reconstructing the coefficient a * and we give an algorithm inspired by the proof of Theorem 1.1. In particular, we shall work under all the assumptions of the previous section. Namely, we assume the geometric condition (1.4), the a priori informations of (1.7) given by m, α 0 , α 1 and ρ, and that the corresponding solution w * of (1.1) satisfies (1.9), while (w 0 , w 1 ) follow (1.10)- (1.11). We shall also assume the knowledge of a * on its boundary, but we will parametrize it slightly differently from (1.8):

a = a * | ∂Ω and b = ∇a * • ∇w 0 | ∂Ω are known.
(

Under these conditions, we assume that T satisfies (1.13) and we choose β ∈ (0, α 0 ρ 0 (Lm/α 0 )) such that

T > L √ β , (1.22)
and then choose ϕ = ϕ(t, x) = |x -x 0 | 2 -βt 2 as in (1.18).

Of course, we also assume that we know the measurement a * ∂ ν w * on (0, T ) × Γ 0 . Since we know a * on the boundary, this is equivalent to know

M * := ∂ ν w * on (0, T ) × Γ 0 , (1.23)
which is the quantity we will work with in the following. The algorithm we propose to reconstruct a * will produce a sequence of iterates that belong to the following set:

V * a,b,α0,α1,ρ,m,w0 := a ∈ W 1,∞ (Ω), ∇ • (a∇w 0 ) ∈ H 1 (Ω), ∇a L ∞ (Ω) m, and
0 < α 0 inf Ω a sup Ω a α 1 , ∇a • (x -x 0 ) 2(1 -ρ)a in Ω,
and a = a and ∇a • ∇w 0 = b on ∂Ω , (1.24)

which will be simply denoted by V * a,b in the following. Note that the trace of ∇a • ∇w 0 makes sense when a ∈ W 1,∞ (Ω) and ∇ • (a∇w 0 ) ∈ H 1 (Ω), since then ∇a • ∇w 0 ∈ H 1 (Ω) as well. This set is a closed variant of V a,aν in (1.12), but it differs mainly due to regularity issues (a ∈ C 1 (Ω) ∩ H 2 (Ω) in (1.12)). In particular, V a,aν is included in V * a,b . Also note that Theorem 1.1 also applied within the class V * a,b instead of V a,aν through easy adaptations of its proof which are left to the reader. Still, we have chosen to present Theorem 1.1 within the class V a,aν since the class V * a,b is more delicate to work with, and depends on the datum w 0 .

To motivate the introduction of this class V * a,b , let us point out that we will consider a projection operator depending on a parameter s 1 on the set V * a,b as follows:

P a,b,s is the projection on V * a,b in {a ∈ H 1 (Ω), ∇ • (a∇w 0 ) ∈ H 1 (Ω)} (1.25)
for the norm

a 2 s = Ω e 2sϕ(0) s 2 |∇a| 2 + s 4 a 2 + |∇ (∇ • (a∇w 0 )) | 2 dx. (1.26)
The important point is that the set V * a,b is closed and convex for the topology induced by the norm • s defined in (1.26), while the set V a,aν is not closed for this topology.

For later use, we choose ε ∈ (0, inf

Ω |x -x 0 | 2 ) and introduce a cut-off function η ∈ C ∞ (R) such that 0 η 1 and η(r) = 1, if r ε, η(r) = 0, if r 0. (1.27)
We also define the constant γ ∇w0 = sign(∇w 0 • (x -x 0 )) (here, the choice of x in Ω is irrelevant due to (1.10)), which shall be used hereafter.

A first reconstruction algorithm for the coefficient a * from the knowledge of ∂ ν w * on (0, T ) × Γ 0 is as follows, and depends on a tuning parameter s 1:

Algorithm 1.
Initialization: Let a 0 be any guess in the set of admissible coefficients

V * a,b . Iteration, from k to k + 1: We suppose that a k ∈ V * a,b is known. • Step 1 -We set μk = η(ϕ)∂ 2 t ∂ ν w[a k ] -∂ ν w * on (0, T ) × Γ 0 , where ∂ ν w * is the measurement and w[a k ] denotes the solution of    ∂ 2 t w -∇ • (a k ∇w) = f, in (0, T ) × Ω, w = f ∂ , on (0, T ) × ∂Ω, w(0) = w 0 , ∂ t w(0) = 0, in Ω, (1.28) corresponding to (1.3) with the coefficient a k . • Step 2 -We minimize the functional J s,a k (z) = 1 2 T 0 Ω e 2sϕ |∂ 2 t z -∇ • (a k ∇z)| 2 dxdt + s 2 T 0 Γ0 e 2sϕ |∂ ν z -μk | 2 dσdt + s 2 {ϕ<0} e 2sϕ (|∂ t z| 2 + |∇z| 2 + s 2 |z| 2 ) dxdt + s 2 Ω e 2sϕ(T ) ∂ t z(T ) 2 + |∇z(T )| 2 + s 2 z(T ) 2 dx (1.29) on T a k = {z ∈ L 2 (0, T ; H 1 0 (Ω)) such that ∂ 2 t z-∇•(a k ∇z) ∈ L 2 ((0, T )×Ω), ∂ t z(0, •) = 0 in Ω, and ∂ ν z ∈ L 2 ((0, T ) × Γ 0 )},
and we denote its unique minimizer by Z k , for which in particular, Z k (0,

•) ∈ H 1 0 (Ω). • Step 3 -We solve ∇ • (δa k (x)∇w 0 (x)) = -Z k (0, x), for x ∈ Ω, δa k = 0, on Γ ∇w0 = {x ∈ ∂Ω, (∇w 0 (x) • ν(x))γ ∇w0 > 0} .
(1.30)

Then set ãk+1 = a k + δa k .

• Step 4 -Finally, using the projection P a,b,s introduced in (1.25), we set

a k+1 = P a,b,s (ã k+1 ). (1.31) ν Ω ∇w 0 Γ ∇w 0 x 0 Figure 1: Illustration of Γ ∇w0 in the case ∇w 0 constant.
An alternative algorithm under the exact same setting can be given as follows:

Algorithm 2.

Initialization: Let a 0 be any guess in the set of admissible coefficients

V * a,b . Iteration, from k to k + 1: We suppose that a k ∈ V * a,b is known. • Step 1 -identical to Algorithm 1. • Step 2 -identical to Algorithm 1. • Step 3 bis -We minimize the functional K s,k (δa) = 1 2 Ω e 2sϕ(0,•) |∇(∇ • (δa∇w 0 )) -∇ Z k (0, •)| 2 dx (1.32)
on the vector space δa ∈ H 1 0 (Ω), ∇δa • ∇w 0 ∈ H 1 0 (Ω) , and denote its unique minimizer by δa k . Then we set ãk+1 = a k + δa k .

• Step 4 -identical to Algorithm 1.

From a numerical point of view, the projection of Step 4 in both algorithm is complicated to perform. Step 3 bis of Algorithm 2 has the advantage over Step 3 of Algorithm 1 of corresponding to a kind of partial projection step. In particular, the boundary conditions contained in V * a,b are satisfied by the minimizer. Consequently, for the numerical implementation in Section 4, we will actually work from Algorithm 2 and Step 4 will be essentially ignored.

We prove that these algorithms both converge, provided the parameter s is chosen large enough: Theorem 1.6. Under the above assumptions, there exist constants C > 0 and s 0 1 such that for all s s 0 , for all a 0 ∈ V * a,b , both Algorithm 1 and 2 are well-defined and the iterates a k constructed in each algorithm satisfy, for all k ∈ N,

a k+1 -a * 2 s C inf{s 2 , e 2s infΩ(ϕ(0)-ε) } a k -a * 2 s . (1.33)
In particular, for s large enough, the sequence a k strongly converges to a * as k → ∞ in the norm • s defined in (1.26).

Theorem 1.6 is the main result of our article, and establishes a way to reconstruct the speed coefficient a * from the knowledge of the measurement ∂ ν w * (and other informations, in particular a priori bounds, geometric conditions, regularity conditions and knowledge of the velocity on the boundary) without an a priori guess. In this sense, using the wording of [START_REF] Beilina | Approximate global convergence and adaptivity for coefficient inverse problems[END_REF], this procedure is globally convergent.

We also want to emphasize that this algorithm can be implemented numerically, and this is presented in Section 4. Let us now point out some of the numerical issues of interest regarding the numerical implementation of these algorithms.

First, the minimization of J s,a k corresponds to the minimization of a strictly convex and coercive quadratic functional. Still, there are some issues in minimizing it due to the fact that the algorithm requires s to be large to be convergent, and therefore one should handle numerically very large exponential terms. This issue is in fact comparable to the one we dealt with in our previous work [START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF] and can be partially solved by working on the conjugate variable and using multiple cut-offs to decompose the minimization of J s,a k in areas in which the exponential terms are of the same order, see Section 4.3.

Second, when considering Algorithm 1, one has to solve the transport equation (1.30) with velocity ∇w 0 , and third, one should compute the projection on the set V * a,b for the topology induced by • s in (1.26). Since the set V * a,b is a convex closed set for this topology, the projection P a,b,s is well-defined. However, computing it in practice seems to be a challenging issue, due to the fact that the elements of V * a,b are delicate to enumerate. We shall therefore prefer to use Algorithm 2 and simply discard the projection step. In fact, we will check numerically that the iterates constructed by Algorithm 2 stay in the class V * a,b even when discarding the projection step. Related references. We shall not give here extensive references to the inverse problems at hand. As we have said, the stability result stated in Theorem 1.1 is closely related to the works [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF] and [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF], in which the assumptions are slightly stronger. In fact, these results, as ours, are based on the original Bukhgeim-Klibanov method [START_REF] Bukhgeım | Uniqueness in the large of a class of multidimensional inverse problems[END_REF], see also [START_REF] Klibanov | Inverse problems and Carleman estimates[END_REF], which relies on the use of Carleman estimates. A survey article with many applications of this strategy can be found in [START_REF] Klibanov | Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems[END_REF].

The main feature of our approach is that it provides a numerical method to reconstruct the unknown coefficient. Although this is a natural question, it turns out that there are very few works addressing these issues in a thorough theoretical setting. Indeed, the problem is strongly non-linear and thus there are no a priori reason to expect the convergence of naive methods such as least squares approach, in particular when there are no good guess of the coefficient to be recovered, due to the existence of possible local minima. Thus, one of the best feature of our algorithm is that it is globally convergent, in the sense that we do not need the initial guess a 0 to be close to the coefficient a * : Indeed, note that, even if we require some partial knowledge on a * , the initial guess a 0 can be very far from a * .

Globally convergent methods are quite rare, but decisive works in this direction have been done since the nineties, in particular based on the idea that Carleman estimates may help in the construction of globally convergent methods, see [START_REF] Klibanov | Global convexity in a three-dimensional inverse acoustic problem[END_REF], for which solving the inverse problem consists in minimizing a strictly convex functional. For a quite recent account on these questions and how it can be used on several models, we refer to the book [START_REF] Beilina | Approximate global convergence and adaptivity for coefficient inverse problems[END_REF]. Since then, several other works have been proposed by Klibanov and coauthors on several models, among which we shall quote in particular the closely related results [START_REF] Beilina | Globally strongly convex cost functional for a coefficient inverse problem[END_REF] for the recovery of a density in a wave equation, and more recently [START_REF] Klibanov | Convexification for the inversion of a time dependent wave front in a heterogeneous medium[END_REF] in the more delicate case in which the initial datum is a Dirac mass (in other context, we also refer to [START_REF] Klibanov | Convexification of electrical impedance tomography with restricted Dirichlet-to-Neumann map data[END_REF] for the recovery of an electrical conductivity in electrical impedance tomography, or [START_REF] Smirnov | On an inverse source problem for the full radiative transfer equation with incomplete data[END_REF] for the recovery of a source term in a radiative transfer equation). In parallel, we have implemented similar ideas in the works [START_REF] Baudouin | Global carleman estimates for waves and applications[END_REF][START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF] for recovering potentials in a wave equation (see also [START_REF] Boulakia | Numerical reconstruction based on Carleman estimates of a source term in a reaction-diffusion equation[END_REF] for the recovery of a source term in a parabolic context). The main novelty in our approach compared to [START_REF] Beilina | Globally strongly convex cost functional for a coefficient inverse problem[END_REF] is that the algorithm we propose consists essentially in minimizing a sequence of quadratic coercive functionals, for which efficient minimization algorithms are available, while [START_REF] Beilina | Globally strongly convex cost functional for a coefficient inverse problem[END_REF] considers one strictly convex functional which is not quadratic and is thus mainly only suitable to classical descent strategies.

Notations. In the rest of the article, C will denote a generic constant which may change from line to line but is always independent of both the Carleman parameter s and the iteration number k.

Lipschitz stability estimate for the inverse problem

Proof of Theorem 1.1. The structure of the proof is basically the same as [START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF] but takes advantage of Carleman estimates (1.20) and (1.17) that are slighlty less demanding than the ones used previously in the litterature.

We work within the setting of Theorem 1.1, having in particular w a solution of (1.3) corresponding to a ∈ V a,aν , and w * solution of (1.1) satisfying the regularity assumption (1.9). We set

u = ∂ 2 t (w -w * ) that satisfies the system    ∂ 2 t u -∇ • (a∇u) = ∇ • ((a -a * )∂ 2 t ∇w * ), in (0, T ) × Ω, u = 0, on (0, T ) × ∂Ω, u(0) = ∇ • ((a -a * )∇w 0 ), ∂ t u(0) = 0, in Ω (2.1)
and we underline that the compatibility of boundary and initial conditions is assured by the assumption that a and a * belong to V a,aν , implying in particular that on ∂Ω, a = a * = a and

∂ ν a = ∂ ν a * = a ν , and therefore u(0, •) ∈ H 1 0 (Ω). (2.2)
Besides, from the assumptions

w * ∈ H 2 (0, T ; W 2,∞ (Ω)) and a, a * ∈ V a,aν , the source term ∇ • ((a - a * )∂ 2 t ∇w * ) belongs to L 2 (0, T ; L 2 (Ω)). Hence it is clear that u ∈ C 0 ([0, T ]; H 1 0 (Ω)) ∩ C 1 ([0, T ]; L 2 (Ω)) and ∂ ν u ∈ L 2 (0, T ; L 2 (∂Ω)), (2.3)
and classical estimates for the wave equation yield

u C 0 ([0,T ];H 1 0 (Ω))∩C 1 ([0,T ];L 2 (Ω)) C( u(0) H 1 0 (Ω) + ∇ • ((a -a * )∂ 2 t ∇w * L 1 (0,T ;L 2 (Ω)) ) C u(0) H 1 0 (Ω) + C a -a * H 1 0 (Ω) . (2.4)
Next, we extend the solution u to negative times by u(t, x) = u(-t, x) for t < 0 so that the wave equation (2.1) is satisfied in (-T, T ).

The proof then consists into two main steps, the first one being the estimation of the difference a -a * by the initial state u(0) in the last line of (2.1) and the second one relying on the estimation of the initial weighed energy of the solution u of system (2.1). Both steps stem from appropriate Carleman estimates.

Estimate of the difference between two coefficients. In order to estimate δa = a -a * from the initial condition in (2.1), we need a Carleman estimate for the following first order operator:

δa → ∇ • (δa∇w 0 ).
It is given by Theorem 1.5 for the variable b = δa and the particular vector field X(x) = ∇w 0 (x) (that satisfies the initial sensitivity condition (1.10)). Since w 0 ∈ W 3,∞ (Ω), a, a * ∈ V a,aν and ∇ • (δa∇w 0 ) = u(0) where u satisfies (2.1), we have δa ∈ H 1 0 (Ω) and u(0) = ∇ • (δa∇w 0 ) ∈ H 1 0 (Ω) and we can write

s 2 Ω e 2sϕ(0) |∇(δa)| 2 + s 2 (δa) 2 dx + Ω e 2sϕ(0) |∇ (∇ • (δa∇w 0 )) | 2 C Ω e 2sϕ(0) |∇u(0)| 2 + s 2 (u(0)) 2 dx, (2.5)
so that the next step is to estimate the right hand side, which is basically the (weighted) initial energy of the solution of (2.1).

Estimate of the initial energy. We assume the geometric condition (1.4) on Γ 0 , the time condition (1.13) on T , and since

x 0 ∈ Ω, we chose 0 < ε < inf x∈Ω |x -x 0 | 2 . Moreover, we use the smooth cut-off function η ∈ C ∞ (R) introduced in (1.27) such that 0 η 1, η(r) = 0 if r 0 and η(r) = 1 if r ε. We then choose β such that L 2 T 2 < β < α 0 ρ 0 Lm α 0 ,
which can be done according to (1.13).

Using then the weight function ϕ defined in (1.18), we set

v = η(ϕ)u = η(ϕ)∂ 2 t (w -w * ) .
Thanks to the support of η • ϕ (see figure 2), v = η(ϕ)u has the property to vanish with its derivatives in {ϕ < 0}. It satisfies

0 r ε 1 η 0 ϕ = 0 ϕ = ε x 0 T η • ϕ = 0 η • ϕ = 1 L Ω
       ∂ 2 t v -∇ • (a∇v) = η(ϕ)∇ • ((a -a * )∂ 2 t ∇w * ) + [∂ 2 t -∇ • (a∇), η(ϕ)]u, in (-T, T ) × Ω, v = 0, on (-T, T ) × ∂Ω, v(0) = ∇ • ((a -a * )∇w 0 ), ∂ t v(0) = 0, in Ω, v(±T ) = 0, ∂ t v(±T ) = 0, in Ω, (2.6
) where [A, B] is the commutator of the operators A and B. Applying the Carleman estimate (1.17) of Theorem 1.4, we get

Ω e 2sϕ(0) (∂ t v(0)) 2 + |∇v(0)| 2 + s 2 (v(0)) 2 dx C T -T Ω e 2sϕ (η(ϕ)∇ • ((a -a * )∂ 2 t ∇w * )) 2 dxdt +C T -T Ω e 2sϕ ([∂ 2 t -∇ • (a∇•), η(ϕ)]u) 2 dxdt + Cs T -T Γ0 e 2sϕ |∂ ν v| 2 dσdt.
Concerning the left hand side of this estimate, we can use that for all x ∈ Ω, v(0, x) = u(0, x). For the boundary term, we have

∂ ν v = η(ϕ)∂ ν ∂ 2 t (w -w * ) in (-T, T
) × ∂Ω and we can bound η(ϕ) by 1. In order to deal with the two other terms of the right hand side, for the first one we use w * ∈ H 2 (-T, T ; W 2,∞ (Ω)) and ϕ(t, x) ϕ(0, x) for all (t, x) ∈ (-T, T ) × Ω, and for the second one, a simple calculation allows to notice that the support of the commutator [∂ 2 t -∇ • (a∇•), η(ϕ)] is confined where ϕ ε (see figure 2). Therefore one can write that

Ω e 2sϕ(0) |∇u(0)| 2 + s 2 (u(0)) 2 dx C Ω e 2sϕ(0) |∇(δa)| 2 + |δa| 2 dx + Ce 2sε T -T Ω (∂ t u) 2 + |∇u| 2 + u 2 dxdt + Cs T -T Γ0 e 2sϕ ∂ ν ∂ 2 t (w -w * ) 2 dσdt (2.7)
and we come back to the time domain (0, T ) by recalling u(-t, •) = u(t, •).

Final estimates. Now, using (2.4), we have

T 0 Ω (∂ t u) 2 + |∇u| 2 + u 2 dxdt C ∇u(0) 2 L 2 (Ω) + C δa 2 H 1 0 (Ω) .
Noticing then that inf Ω ϕ(0, •) > ε, we get e 2sε δa 2

H 1 0 (Ω) Ω e 2sϕ(0) |∇(δa)| 2 + |δa| 2 dx.
Therefore, (2.7) becomes

Ω e 2sϕ(0) |∇u(0)| 2 + s 2 (u(0)) 2 dx C Ω e 2sϕ(0) |∇(δa)| 2 + |δa| 2 dx + Ce 2sε Ω |∇u(0)| 2 dx + Cs T 0 Γ0 e 2sϕ ∂ ν ∂ 2 t (w -w * ) 2 dσdt. (2.8)
Recalling u(0) = ∇ • (δa∇w 0 ), and gathering (2.8) and (2.5) now gives

s 2 Ω e 2sϕ(0) |∇(δa)| 2 + s 2 (δa) 2 dx + Ω e 2sϕ(0) |∇ (∇ • (δa∇w 0 )) | 2 dx C Ω e 2sϕ(0) |∇(δa)| 2 + |δa| 2 dx + Ce 2sε Ω |∇ (∇ • (δa∇w 0 )) | 2 dx + Cs T 0 Γ0 e 2sϕ ∂ ν ∂ 2 t (w -w * ) 2 dσdt. (2.9)
Using then again inf Ω ϕ(0, •) > ε and taking s large enough, we can absorb the terms in δa in the right hand side by the ones on the left hand side and we get

s 2 Ω e 2sϕ(0) |∇(δa)| 2 + s 2 (δa) 2 dx + Ω e 2sϕ(0) |∇ (∇ • (δa∇w 0 )) | 2 dx C T 0 Γ0 e 2sϕ ∂ ν ∂ 2 t (w -w * ) 2 dσdt.
Then, bounding the exponential terms, we obtain the desired Lipschitz stability result

a -a * H 1 0 (Ω) C ∂ ν w -∂ ν w * H 2 (0,T ;L 2 (Γ0))
that concludes the proof of Theorem 1.1.

Convergence of the Algorithm

The core idea of our work is to take advantage of the structure of the proof of Theorem 1.1 to conceive a reconstruction algorithm of the coefficient a * from the knowledge of the flux of the solution. In this section, we will explain how we can prove the convergence of Algorithm 1 using the same tools (Carleman estimates of Theorem 1.4 and Theorem 1.5).

Let us first give the main ideas we develop later in the proper proof of Theorem 1.6. Basically, at the step k, Algorithm 1 uses the fact that if w * is the solution of equation (1.1) and w[a k ] solves (1.28), then

z k = ∂ 2 t w[a k ] -w * (3.1) solves    ∂ 2 t z k -∇ • (a k ∇z k ) = g k , in (0, T ) × Ω, z k = 0, on (0, T ) × ∂Ω, z k (0, •) = z k 0 , ∂ t z k (0, •) = 0, in Ω, (3.2)
where

g k = ∇ • ((a k -a * )∇∂ 2 t w * ), z k 0 = ∇ • ((a k -a * )∇w 0 ), (3.3)
and

µ k = ∂ ν z k on (0, T ) × Γ 0 .
First of all, we can a priori estimate z k with respect to the data z k 0 and g k as follows :

z k L ∞ (0,T ;H 1 0 (Ω)) + ∂ t z k L ∞ (0,T ;L 2 (Ω)) C z k 0 H 1 0 (Ω) + g k L 1 (0,T ;L 2 (Ω)) . (3.4) Indeed, g k ∈ L 1 (0, T ; L 2 (Ω)) and z k 0 ∈ H 1 0 (Ω) thanks to w 0 ∈ H 3 (Ω), w * ∈ H 2 (0, T ; W 2,∞ (Ω)) and a k , a * ∈ V a,aν .
Using the hidden regularity properties of the wave equation [START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF], this also implies that ∂ ν z k is well-defined as an element of L 2 (0, T ; L 2 (∂Ω)), and thus

µ k ∈ L 2 (0, T ; L 2 (∂Ω)). Then we set zk = η(ϕ)z k = η(ϕ)∂ 2 t w[a k ] -w * , (3.5)
which satisfies:

   ∂ 2 t zk -∇ • (a k ∇z k ) = η(ϕ)g k + η(ϕ), ∂ 2 t -∇ • (a k ∇) z k := gk , in (0, T ) × Ω, zk = 0, on (0, T ) × ∂Ω, zk (0, •) = η(ϕ(0))z k 0 , ∂ t zk (0, •) = 0, in Ω, (3.6)
and

μk = ∂ ν zk = η(ϕ)µ k on (0, T ) × Γ 0 . (3.7)
We also have zk (T, [START_REF] Klibanov | Carleman estimates for coefficient inverse problems and numerical applications[END_REF], thanks to the definitions (1.27) of η and (1.18) of ϕ (see the illustration of η • ϕ in Figure 2).

•) = 0 and ∂ t zk (T, •) = 0 in Ω since T √ β > sup x∈Ω |x -x 0 | from (1.
Let us now explain the real insight of the construction of the algorithm. In system (3.2), both the source g k and the initial data z k 0 are unknown. We are ultimately interested in finding a good approximation of z k 0 , which contains informations on a k -a * . In order to do so, we will try to fit "at best" the flux ∂ ν z with µ k on the boundary, approximating the unknown source term g k by 0. This approach works since we can prove that the source term g k brings less information than µ k . This idea is the cornerstone of the construction of the algorithms and convergence results presented in [START_REF] Baudouin | Global carleman estimates for waves and applications[END_REF] and [START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF], and is in fact the main ingredient of the proof of Theorem 1.1. Algorithm 1 is constructed on the minimization of the functional J s,a k defined in (1.29), through which we actually try to approximate zk = η(ϕ)z k , that satisfies the following crucial properties:

• zk (0, •) = η(ϕ(0, •))z k 0 = ∇ • ((a k -a * )∇w 0 ) encodes the information on a k -a * ; • zk = η(ϕ)z k vanishes in {ϕ < 0}; • zk (T, •) = 0 because T > L/ √ β; • ∂ ν zk = μk in (0, T ) × Γ 0 .
Then, as one can understand from the definition of z k 0 in (3.3), we need to access specifically the information on a k -a * . This is the third step of both algorithms. In Algorithm 1, it takes spontaneously the form of studying the first order differential equation (1.30) that encapsulates a k -a * . In Algorithm 2, we propose instead the minimization of a functional, K s,k to approximate a * at best, ensuring at the same time that we stay in an admissible subclass of V * a,b (in particular concerning the boundary conditions). Finally, Step 4 is a very formal final step to guarantee that each iterative candidate a k remains in the admissible space V * a,b .

Most of these ideas take actually their roots in the proofs of the stability of similar inverse problems by compactness uniqueness arguments [START_REF] Puel | On a global estimate in a linear inverse hyperbolic problem[END_REF][START_REF] Puel | Generic well-posedness in a multidimensional hyperbolic inverse problem[END_REF][START_REF] Yamamoto | Uniqueness and stability in multidimensional hyperbolic inverse problems[END_REF] or by Carleman estimates [START_REF] Imanuvilov | Global Lipschitz stability in an inverse hyperbolic problem by interior observations[END_REF][START_REF] Imanuvilov | Global uniqueness and stability in determining coefficients of wave equations[END_REF][START_REF] Imanuvilov | Determination of a coefficient in an acoustic equation with a single measurement[END_REF][START_REF] Baudouin | Lipschitz stability in an inverse problem for the wave equation[END_REF]. Of course, Section 2, with the detailed proof of the stability result has already given a good idea of the technical points that allow to construct the proof of the algorithms convergence result.

Proof of Theorem 1.6. We place ourselves in the setting of Theorem 1.6. The proof below will be detailed first in the context of Algorithm 1. The analysis of Algorithm 2 will be given afterwards.

Stage 1. Well-posedness of the algorithms. For T satisfying (1.22) and a ∈ V * a,b , we consider the space

T a = z ∈ C 0 ([0, T ]; H 1 0 (Ω)) ∩ C 1 ([0, T ]; L 2 (Ω)), with ∂ 2 t z -∇ • (a∇z) ∈ L 2 ((0, T ) × Ω) and ∂ t z(0, •) = 0 in Ω , (3.8)
endowed with the norm

z 2 obs,s,a = T 0 Ω e 2sϕ |∂ 2 t z -∇ • (a∇z)| 2 dxdt + s T 0 Γ0 e 2sϕ |∂ ν z| 2 dσdt + {ϕ<0} e 2sϕ (|∂ t z| 2 + |∇z| 2 + s 2 |z| 2 ) dxdt + s Ω e 2sϕ(T ) (∂ t z(T )) 2 + |∇z(T )| 2 + s 2 z(T ) 2 dx.
This semi-norm is actually a norm thanks to the Carleman estimate of Theorem 1.4, for all s large enough. Since the weight function ϕ is bounded from above and below on [0, T ] × Ω, it is also a norm for all s > 0. One can read details and remarks on that matter in reference [START_REF] Baudouin | Global carleman estimates for waves and applications[END_REF]. Therefore, if a k ∈ V * a,b , the functional J s,a k defined in (1.29) is clearly continuous, strictly convex and coercive on the space T a k and consequently admits a unique minimizer. What remains of an iteration of Algorithm 1 is clearly well defined, and the last step guarantees that a k+1 ∈ V * a,b , so that Algorithm 1 is well-defined for any choice of a 0 ∈ V * a,b and s > 0.

Stage 2. Dependence of the minimizer Zk with respect to the source term. To show that the minimizer Zk of J s,a k is close to zk in (3.5), we shall use the fact that they correspond to minimizers of similar functionals. Thus, for generic a ∈ V (defined in (1.16)), µ ∈ L 2 (0, T ; L 2 (Γ 0 ))

and g ∈ L 2 ((0, T ) × Ω), we introduce the functional

J s,a [µ, g](z) = 1 2 T 0 Ω e 2sϕ |∂ 2 t z -∇ • (a∇z) -g| 2 dxdt + s 2 T 0 Γ0 e 2sϕ |∂ ν z -μ| 2 dσdt + s 2 {ϕ<0} e 2sϕ (|∂ t z| 2 + |∇z| 2 + s 2 |z| 2 ) dxdt + s 2 Ω e 2sϕ(T ) ∂ t z(T ) 2 + |∇z(T )| 2 + s 2 z(T ) 2 dx (3.9)
defined on T a (defined in (3.8)).

The following result describes how the minimizer of the functional J s,a [µ, g] depends on the source term g in the wave equation. Proposition 3.1. Assume that a ∈ V (defined in (1.16)), µ ∈ L 2 ((0, T ) × Γ 0 ) and both g and ḡ belong to L 2 ((0, T ) × Ω). Let Z (resp. Z) be the unique minimizer of the functional J s,a [µ, g] (resp. J s,a [µ, ḡ]). Then there exist positive constants s 0 and C independent of a ∈ V such that for s s 0 we have:

Ω e 2sϕ(0) |∇Z(0) -∇ Z(0)| 2 dx + s 2 Ω e 2sϕ(0) (Z(0) -Z(0)) 2 dx C T 0 Ω e 2sϕ |g -ḡ| 2 dxdt, (3.10)
where ϕ and s 0 are chosen as in Theorem 1.4.

Proof of Proposition 3.1. If Z is the unique minimizer of the functional J s,a [µ, g], then the Euler Lagrange equation satisfied by Z writes: For all z ∈ T a ,

T 0 Ω e 2sϕ (∂ 2 t Z -∇ • (a∇Z) -g)(∂ 2 t z -∇ • (a∇z)) dxdt + s T 0 Γ0 e 2sϕ (∂ ν Z -µ)∂ ν z dσdt + s {ϕ<0} e 2sϕ (∂ t Z∂ t z + ∇Z • ∇z + s 2 Zz) dxdt + s Ω e 2sϕ(T ) ∂ t Z(T )∂ t z(T ) + ∇Z(T ) • ∇z(T ) + s 2 Z(T )z(T ) dx = 0. (3.11)
A similar Euler Lagrange equation can be derived for Z minimizer of J s,a [µ, ḡ]. We apply both identities to z = Z -Z and subtract one from another. We obtain:

T 0 Ω e 2sϕ |∂ 2 t z -∇ • (a∇z)| 2 dxdt + s T 0 Γ0 e 2sϕ |∂ ν z| 2 dσdt + s {ϕ<0} e 2sϕ (|∂ t z| 2 + |∇z| 2 + s 2 |z| 2 ) dxdt + s Ω e 2sϕ(T ) (∂ t z(T )) 2 + |∇z(T )| 2 + s 2 z(T ) 2 dx = T 0 Ω e 2sϕ (g -ḡ)(∂ 2 t z -∇ • (a∇z)) dxdt.
Using Young's inequality, it implies that

1 2 T 0 Ω e 2sϕ |∂ 2 t z -∇ • (a∇z)| 2 dxdt + s T 0 Γ0 e 2sϕ |∂ ν z| 2 dσdt + s {ϕ<0} e 2sϕ (|∂ t z| 2 + |∇z| 2 + s 2 |z| 2 ) dxdt + s Ω e 2sϕ(T ) (∂ t z(T )) 2 + |∇z(T )| 2 + s 2 z(T ) 2 dx 1 2 T 0 Ω e 2sϕ |g -ḡ| 2 dxdt. (3.12)
But the left hand side of (3.12) is precisely the right hand side of the Carleman estimate (1.17). Hence, applying Theorem 1.4, which applies under the assumptions of Proposition 3.1, to z = Z -Z which satisfies ∂ t z(0, •) = 0 (to be perfectly accurate, we apply Theorem 1.4 to the extension of z to (-T, T ) × Ω by z(-t, •) = z(t, •)), we immediately deduce estimate (3.10).

Applying Proposition 3.1, we deduce that

Ω e 2sϕ(0) |∇ Z k (0)-∇z k (0)| 2 dx+s 2 Ω e 2sϕ(0) ( Z k (0)-zk (0)) 2 dx C T 0 Ω e 2sϕ |g k | 2 dxdt, (3.13)
where we used that Zk , respectively zk in (3.5), is the minimizer of

J s,a k [μ k , 0], respectively of J s,a k [μ k , gk ],
with gk as in (3.7).

Stage 3. Estimation through the first order differential operator. Now, by design in Algorithm 1, ãk+1 satisfies

∇ • ((ã k+1 -a k )∇w 0 ) = -Z k (0, •), in Ω, ãk+1 -a k = 0, on Γ ∇w0 = {x ∈ ∂Ω, (∇w 0 (x) • ν(x))γ ∇w0 > 0} , (3.14)
while zk (0,

•) = η(ϕ(0, •))z k 0 = ∇ • ((a k -a * )∇w 0 ) in Ω.
Therefore, the estimate (3.13) is in fact an estimate on zk (0,

•) -Z k (0, •) = ∇ • ((ã k+1 -a * )∇w 0 ). Now, turning to the definition of gk in (3.6), since η(ϕ), ∂ 2 t -∇ • (a k ∇•)
has support in a region where ϕ ε, we obtain

Ω e 2sϕ(0) |∇(∇ • ((ã k+1 -a * )∇w 0 )| 2 dx + s 2 Ω e 2sϕ(0) (∇ • ((ã k+1 -a * )∇w 0 ) 2 dx C T 0 Ω e 2sϕ |g k | 2 dxdt C T 0 Ω e 2sϕ η(ϕ)∇ • ((a k -a * )∇∂ 2 t w * ) + η(ϕ), ∂ 2 t -∇ • (a k ∇) z k 2 dxdt C w * 2 H 2 (0,T ;W 2,∞ (Ω)) Ω e 2sϕ(0) |a k -a * | 2 + ∇(a k -a * ) 2 dx + Ce 2sε T 0 Ω |z k | 2 + |∇z k | 2 + |∂ t z k | 2 dxdt.
Notice here that the bound of η(ϕ),

∂ 2 t -∇ • (a k ∇) independently of k stems from the a priori estimates a L ∞ (Ω) α 1 and ∇a k L ∞ (Ω) m in the set V * a,b .
Using the classical energy estimate recalled in (3.4), and since inf Ω ϕ(0,

•) > ε in Ω, we get Ω e 2sϕ(0) |∇(∇ • ((ã k+1 -a * )∇w 0 )| 2 dx + s 2 Ω e 2sϕ(0) (∇ • ((ã k+1 -a * )∇w 0 ) 2 dx C w * 2 H 2 (0,T ;W 2,∞ (Ω)) Ω e 2sϕ(0) |a k -a * | 2 + ∇(a k -a * ) 2 dx + Ce 2sε ∇ • ((a k -a * )∇w 0 ) 2 H 1 0 Ω) + ∇ • ((a k -a * )∇∂ 2 t w * ) 2 L 1 (0,T ;L 2 (Ω)) C w 0 W 2,∞ (Ω) + w * H 2 (0,T ;W 2,∞ (Ω)) 2 Ω e 2sϕ(0) |a k -a * | 2 + ∇(a k -a * ) 2 dx + Ce 2sε Ω ∇(∇ • ((a k -a * )∇w 0 )) 2 dx.
Thus, for a constant C independent of k and s, we can write

Ω e 2sϕ(0) |∇(∇ • ((ã k+1 -a * )∇w 0 )| 2 dx + s 2 Ω e 2sϕ(0) (∇ • ((ã k+1 -a * )∇w 0 ) 2 dx C Ω e 2sϕ(0) |a k -a * | 2 + ∇(a k -a * ) 2 dx + Ce 2sε Ω ∇(∇ • ((a k -a * )∇w 0 )) 2 dx. ( 3
b ∈ H 1 ∇w0 (Ω) and ∇ • (b∇w 0 ) ∈ H 1 ∇w0 (Ω). Indeed, ∇•(b∇w 0 ) = -Z k (0)+∇((a k -a * )∇w 0 ) actually belongs to H 1 0 (Ω).
The function b is also vanishing on the boundary Γ ∇w0 since δa k = ãk+1 -a k = 0 on Γ ∇w0 while a k -a * = 0 on ∂Ω since a k and a belong to V * a,b . Therefore, using Theorem 1.5 and (3.15), we obtain

s 2 Ω e 2sϕ(0) |∇(ã k+1 -a * )| 2 + s 2 |ã k+1 -a * | 2 dx + Ω e 2sϕ(0) |∇ ∇ • ((ã k+1 -a * )∇w 0 ) | 2 dx C Ω e 2sϕ(0) |a k -a| 2 + ∇(a k -a * ) 2 dx + Ce 2sε Ω ∇(∇ • ((a k -a * )∇w 0 )) 2 dx, which in particular implies that ãk+1 ∈ {a ∈ H 1 (Ω), ∇ • (a∇w 0 ) ∈ H 1 (Ω)}.
Using here the norm defined in (1.26) and since we have ϕ(0) > ε, this can be rewritten as

ãk+1 -a * 2 s C inf{s 2 , e 2s infΩ(ϕ(0)-ε) } a k -a * 2 s
Finally, since a k+1 = P a,b,s (ã k+1 ), where the projection P a,b,s is associated to the norm • s and a * ∈ V * a,b , we have

a k+1 -a * 2 s ãk+1 -a * 2
s , so that it yields

a k+1 -a * 2 s C inf{s 2 , e 2s infΩ(ϕ(0)-ε) } a k -a * 2 s , (3.16)
where the constant C is independent of k and s s 0 .

In particular, if we chose s large enough such that C/inf{s 2 , e 2s infΩ(ϕ(0)-ε) } = q < 1, then we obtain

∀k ∈ N, a k+1 -a * 2 s q a k -a 2 s , and then a k -a * 2 s q k a 0 -a * 2 s , (3.17)
which clearly goes to 0 as k → ∞. This concludes the proof of the convergence of Algorithm 1.

Proof of convergence for algorithm 2. Most of Algorithm 2 is similar to Algorithm 1, and we thus only focus on the differences between the two.

Stage 1 bis. Estimation through another minimization. The well-posedness of Algorithm 2 follows from the well-posedness of the minimization of the functional J s,a k on T a k as before, and from the well-posedness of the minimization of the functional

K s,k in (1.32) on C = δa ∈ H 1 0 (Ω), ∇δa • ∇w 0 ∈ H 1 0 (Ω) endowed with δa 2 C = Ω e 2s|x-x0| 2 |∇(∇ • (δa∇w 0 ))| 2 dx.
This is indeed a norm thanks to the Carleman estimate of Theorem 1.5, for all s large enough, since X = ∇w 0 is given under assumption (1.10) and δa ∈ C. Therefore, the functional K s,k defined in (1.32) is continuous, strictly convex and coercive on the space C and consequently admits a unique minimizer δa k .

Stage 2 bis. Dependence of the minimizer Zk with respect to the source term. This step is exactly as in Algorithm 1.

Stage 3 bis. Estimation through another minimization. In Algorithm 2, the minimization of K s,k by δa k ensures that K s,k (δa k ) K s,k (a k -a * ) and since zk (0,

•) = ∇ • ((a k -a * )∇w 0 ), this rewrites Ω e 2sϕ(0,•) |∇(∇ • (δa k ∇w 0 )) -∇ Z k (0, •)| 2 dx Ω e 2sϕ(0) |∇ Z k (0) -∇z k (0)| 2 dx
From the weighted Poincaré estimate of Lemma A.2, this yields

Ω e 2sϕ(0,•) |∇(∇ • (δa k ∇w 0 )) -∇ Z k (0, •)| 2 dx + s 2 Ω e 2sϕ(0,•) |∇ • (δa k ∇w 0 ) -Z k (0, •)| 2 dx C Ω e 2sϕ(0) |∇ Z k (0) -∇z k (0)| 2 dx
Hence, from (3.13) and the estimations of gk in the previous stage, we get (3.15) immediately.

Conclusion.

We can now conclude the proof of the convergence of Algorithm 2 as before, following Stage 4 of the proof of the convergence of Algorithm 1.

Numerical issues

In this subsection, we present the discretization schemes and the technical solutions we have developed to implement numerically the algorithm. In order to simplify the presentation, from now on, we focus on the one-dimensional case Ω = (0, L x ) and Γ 0 = {x = L x }. Note that in this case, our problem can probably be solved in an easier way using the specific features of the 1-d wave equation (see e.g. [18, Section 8.1]), but our goal is to propose an approach that is generalizable to the multi-dimensional setting.

Generation of the data and noise

In this article, we work with synthetic data. To discretize the wave equation (1.1) in (0, T )×(0, L x ), we use a finite differences scheme in space and time domains. The space and time steps are denoted by h and τ respectively. We introduce N x and N t such that L x = (N x + 1)h and T = (N t + 1)τ , and we define, for 0 j N x + 1 and 0 n N t + 1, (w * ) n j a numerical approximation of the solution w * (t n , x j ) with t n = nτ and x j = jh. It satisfies the following system:

         (w * ) n+1 j -2(w * ) n j + (w * ) n-1 j τ 2 -(∆ a,h (w * ) n ) j = f (t n , x j ), 1 j N x , 1 n N t , (w * ) 0 j = w 0 (x j ) and (w * ) 1 j = w 0 (x j ) + τ 2 2 ((∆ a,h w 0 )(x j ) + f (0, x j )) , 1 j N x , (w * ) n 0 = f ∂ (t n , 0) and (w * ) n Nx+1 = f ∂ (t n , L x ), 0 n N t + 1, (4.1)
with

(∆ a,h w) j = 1 h a j+1 + a j 2 w j+1 -w j h - a j + a j-1 2
w j -w j-1 h .
Then, we compute the discrete counterpart of the continuous measurement M * given in (1.23) as follows:

M n * = (w * ) n Nx+1 -(w * ) n Nx h , 0 n N t + 1.
In the applications we have in mind, the entry data for the inverse problem are directly the second derivative in time of the flux M * , measured for instance with an accelerometer. To compute this discrete time derivative, instead of the classical finite difference formulae that generates instabilities, we use a Savitzki-Golay formula [START_REF] Savitzky | Smoothing and Differentiation of Data by Simplified Least Squares Procedures[END_REF] associated with a cubic polynomial and a window size of 9 points to guarantee the smoothness of the result. We set:

A n * = SavGol(M * ) n , 0 n N t + 1.
that is an approximation of the acceleration A * := ∂ 2 t M * at time t n = nτ . On these measured data, we may add a Gaussian noise:

A n * ←-A n * + θ(max A * )N (0, 1), 0 n N t + 1 (4.2)
where N (0, 1) satisfies a centered normal law with deviation 1 and α is the level of noise. Note that the chosen model of noise is additive. In the presence of noise, the first step of the inverse problem is to regularize the data in a way that depends on the nature of the noise and the physical phenomena. Here, we use a 3-order low-pass Butterworth filter [START_REF] Butterworth | On the Theory of Filter Amplifiers[END_REF] associated to a cutoff frequency that should be known a priori.

Remark 4.1. In order to avoid the inverse crime, we don't use the same schemes for the direct and the inverse problems. Hence, we solve (4.1) thanks to an explicit scheme and we use an implicit scheme for equation (4.4), see afterwards. This implies that even in the absence of additional noise the reconstruction has to deal with numerical errors.

Discretization of Algorithm 2

Of course, our goal is to mimic, as much as possible, Algorithm 1 or Algorithm 2. Still, we will see that each step requires to be adapted to the discrete settings below, even if the whole structure of the algorithm will be preserved. As a matter of fact, from numerical experiments, it seems that Algorithm 2 is more suitable numerically, and we shall thus discuss only the discretization of Algorithm 2.

Initialization process.

This is the easiest step to be discretized. Somehow, it simply consists in exhibiting a coefficient a 0 lying in the class V a,bν . But even in the continuous case, this is not so simple. To handle the task, we introduce a functional

K init (a) = Ω |∇(∇ • (a∇w 0 ))| 2 dx (4.3)
among all a ∈ C a,bν defined by

C a,bν = {a ∈ H 1 (Ω), ∇ • (a∇w 0 ) ∈ H 1 (Ω)
, with a = a and ∇a • ∇w 0 = b ν on ∂Ω}, and we let a 0 be the minimizer of K init . Note that nothing guarantees that a 0 belongs to the class V * a,b , in particular regarding condition (1.5). Still, when testing our numerical algorithms, it did not seem to bring any problem. The discretization can be done classically, for instance using finite differences or finite element methods. We decided to use finite differences in our numerical simulations.

Remark 4.2. In 1d, the minimization of K init can even be done by hand. Indeed, if Ω = (0, L x ), since ∂ x w 0 is assumed to be non-vanishing on (0, L x ), we easily get that the Euler Lagrange equation satisfied by the minimizer a 0 yields ∂ 4

x (a 0 (x)∂ x w 0 (x)) = 0 for all x ∈ (0, L x ), so that a 0 (x)∂ x w 0 (x) is a polynomial of degree 3 which is fully determined by the four boundary conditions a 0 | x=0,Lx , and

∂ x a 0 ∂ x w 0 | x=0,Lx .
Iterative steps.

Discretization of Step 1. It simply consists in computing a discrete version of w in (1.28):

         w n+1 j -2w n j + w n-1 j τ 2 -(∆ a k ,h w n+1 ) j = f (t n+1 , x j ), 1 j N x , 1 n N t , w 0 j = w 0 (x j ) and w 1 j - τ 2 2 (∆ a k ,h w 1 ) j = w 0 (x j ) + τ 2 2 f (dt, x j ), 1 j N x , w n 0 = f ∂ (t n , 0) and w n Nx+1 = f ∂ (t n , L x ), 0 n N t + 1. (4.4)
Then, in order to compute a discrete version of μk , we define the cut-off function η as:

η(τ ) =        1, if τ ε, 0, if τ 0, 1 - ε τ f (t)dt ε 0 f (t)dt , if 0 < τ < ε, with f (t) = exp -1 t(ε -t) , (4.5)
and we set

μk,n = η(ϕ(t n , L x )) SavGol w Nx+1 -w Nx h n -A n * , 0 n N t + 1.
Discretization of Step 2. The minimization process for J s,a k [μ k , 0] is equivalent to the resolution of the following variational formulation: Find Z ∈ T a k such that for all z ∈ T a k ,

T 0 Ω e 2sϕ ∂ 2 t Z -∇ • (a k ∇Z) ∂ 2 t z -∇ • (a k ∇z) dxdt + s T 0 Γ0 e 2sϕ ∂ ν Z∂ ν z dσdt + s Ω e 2sϕ(T ) ∂ t Z(T )∂ t z(T ) + ∇Z(T ) • ∇z(T ) + s 2 Z(T )z(T ) dx + s {ϕ<0} e 2sϕ ∂ t Z∂ t z + ∇Z • ∇z + s 2 Zz dxdt = s T 0 Γ0 e 2sϕ μk ∂ ν z dσdt. (4.6)
To discretize (4.6), we approximate the integrals using rectangle quadrature rules and we use standard centered finite differences for the various integrands. We do not detail here these developments.

As far as the discrete setting is concerned, we have to add new terms to the cost functional to guaranty its coercivity property uniformly with respect to the discretization parameters τ and h. These additional terms (sometimes called viscosity terms) that help handling high frequency spurious waves generated by the numerical scheme take the form

shτ Nt n=0 Nx j=0 e 2sϕ(t n ,xj ) z n+1 i+1 -z n+1 i -z n i+1 + z n i 2 . (4.7)
For explanations and proofs, we refer to the works [START_REF] Baudouin | Convergence of an inverse problem for a 1-D discrete wave equation[END_REF][START_REF] Baudouin | Stability of an inverse problem for the discrete wave equation and convergence results[END_REF] on discrete Carleman estimates for the wave operator. Finally the minimizer Z k is obtained as the unique solution of a linear system of the form AZ = b (4.8) associated with a sparse matrix A of size (N t + 2)(N x + 2) corresponding to the discretization of the variational problem (4.6), that is solved using a Conjugate Gradient approach.

Discretization of Step 3. Here, we discretize the functional K s,k defined by (1.32) using finite difference as for the rest of the algorithm. This will allow to compute a minimizer δa k efficiently.

Remark 4.3. Similarly as in Remark 4.2, in 1-d, i.e. when Ω = (0, L x ), the minimizer δa k can be computed explicitly. Indeed, the Euler-Lagrange equation yields that

∂ xx (e 2ϕ(0,x) (∂ x (∂ x (δa k (x)∂ x w 0 (x)) -Z k (0, x)))) = 0 for all x ∈ (0, L x ).
It follows that for some coefficients γ 0 , γ 1 , and for all x ∈ (0, L x ),

∂ xx (δa k (x)∂ x w 0 (x)) = ∂ x Z k (0, x) + e -2sϕ(0,x) (γ 0 + γ 1 x).
Integrating, we get

∂ x (δa k (x)∂ x w 0 (x)) = Z k (0, x) + γ 0 x 0 e -2sϕ(0,x ) dx + γ 1 x 0
x e -2sϕ(0,x ) dx , and

δa k (x)∂ x w 0 (x) = x 0 Z k (0, x ) dx + γ 0 x 0 (x -x )e -2sϕ(0,x ) dx + γ 1 x 0 (x -x )x e -2sϕ(0,x dx .
Therefore, the solution δa k can be computed explicitly, by choosing γ 0 and γ 1 such that

γ 0 Lx 0 e -2sϕ(0,x) dx + γ 1 Lx 0 xe -2sϕ(0,x) dx = 0, γ 0 Lx 0 (L x -x)e -2sϕ(0,x) dx + γ 1 Lx 0 (L x -x)xe -2sϕ(0,x dx = - Lx 0 Z k (0, x) dx,
which will guarantee that δa k satisfies the appropriate boundary conditions.

Discretization of Step 4. Until now we did not manage to build an efficient way to compute the projection on the set of admissible coefficients V * a,b endowed with the norm • s , and so we simply discard this step. In fact, this is what makes us choose Algorithm 2 over Algorithm 1. Indeed, the step 3 of Algorithm 1 only guarantees that ãk+1 = ãk in Γ ∇w0 , while the step 3 of Algorithm 2 guarantees that ãk+1 = ãk on the whole boundary ∂Ω, and is in some sense closer to be in the correct class V * a,b . Currently, we only assert that at each step a k+1 = ãk+1 remains in V a,aν by checking numerically all the corresponding conditions. Fortunately, in all the examples we considered, these conditions are satisfied at each iteration, if the exact coefficient a * satisfies (1.7).

Stopping criterion.

The iterative loop is stopped when one of the three following criteria is satisfied:

μn 2 A * 2 max( , α) or a k+1 -a k ∞ α 1 , (4.9)
where is a fixed tolerance and α is the level of noise (recall (4.2)), or when the maximal number of iterations is reached. If the exact solution a * is known, we can compute the final error

err a = a ∞ -a * 2 a * 2 , (4.10)
with a ∞ the converged numerical solution. This error ratio is reported in Table 2.

Numerical challenges

The main drawback of the approach is related to the presence of large exponential factors in the functional J s,a k [μ] that leads to severe numerical difficulties when performing the minimization for s large. As in [START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF], this difficulty is solved thanks to the three following ideas:

• We prove and use Carleman estimates with a single weight of the form e 2sϕ instead of the most frequently used double weight e 2se λϕ for λ > 0;

• At the discrete level, we work on the conjugate variable Y n j = Z n j e sϕ(t n ,xj ) , for 0 n N t + 1 and 0 j N x + 1. This change of unknown acts as a preconditioner of A presented in (4.8) that decreases its condition number. Indeed, it allows to remove all the exponential terms appearing in the matrix A. Unfortunately, there are still exponentials in the right hand side vector b.

• The idea to tackle the exponential factor in the right hand side b is to develop a progressive process to compute the solution Y as the aggregation of several problems localized in subdomains in which the exponential factors are all of the same order. In this objective, from the smooth cut-off function η equal to 1 for r ε defined in (1.27), we introduce N s cut-off functions {η i } 1 i Ns (these ones are not necessarily smooth) such that ∀r ∈ R,

Ns i=1 η i (r) = η(r), (4.11)
as illustrated in Figure 3.

0 r ε 1 η 3 η 2 η 1 η Figure 3: Example of cut-off functions η i for 1 i 3.
One immediately gets that, if for each i ∈ {1, • • • , N s }, we denote by Y i the solution of the linear system associated to μi = η i (ϕ)µ, then the solution Y associated to μ = η(ϕ)µ is simply given by

Y = Ns i=1 Y i .
The interest of this approach is that the target flux η i (ϕ)µ involves exponential terms in ϕ on the support of η i (ϕ(t, x)). This becomes particularly interesting if we impose that for each

i ∈ {1, • • • , N s }, Supp η i ⊂ [c i , b i ] with b i -c i C, (4.12)
for some given constant C > 0. Indeed, in that case, we get

sup Supp ηi(ϕ) e sϕ inf Supp ηi(ϕ)
e sϕ e sC , so that if C 1/s, all the exponentials are of the same order when computing μi . This approach can also be used, at least theoretically, to parallelize the minimization of J s,a k [μ, 0].

Numerical results

This subsection is devoted to the presentation of some numerical examples to illustrate the properties of the reconstruction algorithm and its efficiency. All simulations are executed with Python and the source codes are available on request. Table 1 and Figure 4 gather the numerical values used for all the following examples, unless specified otherwise where appropriate. In all the figures showing reconstruction results, the exact coefficient that we want to recover is plotted by a red line, whereas the numerical coefficient recovered by the algorithm is represented by a dotted black line. The convergence informations (number of iterations, running time, convergence error) are reported in Table 2. For the explicit scheme used in (4.1) to remain stable, one has to impose the following CFL condition on the time and space steps:

L x N x w 0 f ∂ f 1 200L x x x
√ α 1 τ h 1,
where √ α 1 is the maximum value of the wave speed. Thus, we set

N t = T (Nx+1) √ α1
Lx where • stands for the integer part. Table 2: Convergence results of the test cases. One can observe that in all cases the convergence criteria (4.9) for = 1 • 10 -3 is met in less than 4 iterations. The algorithm is running on a personal laptop in a few minutes. The final error err a defined in (4.10) shows the effectiveness of the reconstruction.

* = ∂ 2 t w * ,    ∂ 2 t u * -∇ • (a * ∇u * ) = ∂ 2 t f, in (0, T ) × Ω, u * = ∂ 2 t f ∂ , on (0, T ) × ∂Ω, u * (0) = f (0) + ∇ • (a * ∇w 0 ), ∂ t u * (0) = ∂ t f (0), in Ω. (4.13) (a) w 0 (x) = x (b) f ∂ (t) = x cos(πt) (c) f (0, x)
Using the regularity result [26, Theorem 4.1], we observe that a compatibility condition

∂ 2 t f ∂ (0) = f (0) + ∇ • (a * ∇w 0 ), on ∂Ω, (4.14)
has to be imposed to guarantee the regularity condition (1.9) on w * . Numerically, given f ∂ , w 0 , a and b, we compute f 0 as a lifting of the compatibility condition (4.14) defined on ∂Ω by solving

-∆f 0 (x) = 0, x ∈ Ω, f 0 (x) = ∂ 2 t f ∂ (0, x) -a(x)∆w 0 (x) -b(x), x ∈ ∂Ω,
and choose f (t, x) = f 0 (x), for all (t, x) ∈ (0, T ) × (0, L x ).

Simulations from data without noise

In this subsection, several results of reconstruction obtained in the absence of noise are given.

Figure 5 shows the evolution of the recovered coefficient a k at each iteration k of the convergence process and below the final result a ∞ together with the exact solution a * . Figure 6 illustrates the progressive process (4.11) on the first iteration of the algorithm for the recovery of a * (x) = 6 + exp

-0.3 1-8|x-0.4| 2 1 1>8|x-0.4| 2 .
Starting from a 0 0 = a 0 , we represent successively where δa 0 i is the solution returned by Step 3 of Algorithm 2 associated to the measurements η i (ϕ)µ. In this example, T is given by (1.13), s = 10 and N s = 100 and we observe that only the 8 first cut-off functions bring information for the reconstruction, which corresponds approximately to the time condition (1.6). As we will see below, see Figure 9, numerical simulations seem to indicate that the time given by (1.6) is the one needed to guarantee the convergence of the algorithm, thus supporting the idea that the time condition in (1.13), much larger, is mainly technical. 2 remain proportional to θ.

a 0 i = a 0 i-1 + δa 0 i , 1 i N s , (4.15)

Wrong choices of the parameters

The first two columns Figure 9 present several reconstruction results when the assumption on the minimal observation time T given in (1.13) is not satisfied. In the first column, the observation time is taken equal to the intuitive physical time T = Lx √ α0 , corresponding to the time needed for the slowest wave (the one travelling at speed √ α 0 ) to reach the observation zone {x = L x } even if it starts from the opposite side of the domain {x = 0}. In that case, the reconstruction is as accurate as the one obtained with the much larger time given in (1.13). This observation is also illustrated in Figure 6 when only the first cut-off functions are needed to recover the solution. Then, in the second column, we try the same reconstructions but with a shorter time T = 0.8 Lx √ α0 . In that case, the result begins to deteriorate near the boundary {x = 0}. This illustrates the fact that the minimal time to get the convergence of the algorithm is very likely Lx √ α0 , instead of the more restrictive time condition (1.13). a,b because it does not satisfy the condition (1.5) for any ρ > 0. The time T for this example is computed using (1.13) formula with ρ 0 = 0.01.

In the last column of Figure 9 (column (c)), we also plot the results given by the algorithm in the 1-dimensional setting when the exact coefficient to be recovered does not satisfy the multiplier condition (1.5) for any ρ > 0. Note that in this case, of course, it is important to discard the projection step 4 of Algorithm 2. We also choose a very large time T = L/0.01. The numerical experiments seem to indicate good convergence of the algorithm in this 1-dimensional case, although it is used in a setting in which we do not know how to provide a proof of convergence. However, in the 2-dimensional case, we have found some examples of velocity (typically of the form a * (x) = |x -x 0 | 4 ã * (θ)), which do not satisfy (1.5) and for which the algorithm does not converge. This dichotomy between the 1-dimensional case and the 2-dimensional case is somehow expected as the geometry is much more complex in dimension larger than two: in particular, there might be rays of geometric optics that will never meet the observation set, thus violating the observability inequality in this case.

Simulations in two dimensions in space

We also performed some reconstructions procedures in Ω = [0, 1]×[0, 1.1], discretized with 40 points in each direction, x 0 = (-0.3, 0.55), s = 5 and Γ 0 = {x 2 = 0} ∪ {x 1 = 1} ∪ {x 2 = 1.1}. Figures 10 and 11 present the results obtained for two different coefficients a * in the absence of noise. In each case, the gray scales are identical for the exact and the recovered graphics. The convergence informations are reported in Table 2.

Conclusion

The reconstruction approach we presented theoretically and implemented numerically in this article has the interest of being constructed on a quite well-spread methodology for the proof of Lipschitz stability of coefficient inverse problems. It has every chance to be generalizable to situations where the so called Bukhgeim-Klibanov method is used to prove stability via Carleman estimates. It should therefore bring adapted globally convergent algorithms for several other reconstruction of partial differential equations coefficients.

Still, some numerical difficulties remain, essentially linked to the 2D setting and the size of the Carleman parameter s with respect to the mesh size h. The management of the level of noise is also an issue when s gets large. The next step of our work is now to seek for more practical applications and the challenge of dealing with the corresponding alleged constraints.

Nevertheless, we believe that the generalizability of the approach, the fact that the algorithm is proved to be globally convergent from any initial guess, and the extend of our study from theory to numerics are strong features of our work. A Main tools : Carleman estimates.

A.1 Carleman estimates for the wave equation

In order to deal with the effective reconstruction of the speed in a wave equation, we need to prove a specific one parameter Carleman estimate for the wave equation with non constant speed, which will keep track of the multiplier type condition (1.5). In order to do that, we follow the proof of [22, Theorem 2.24 p.45] in our case, in a spirit close to the one in [START_REF] Baudouin | Convergent algorithm based on carleman estimates for the recovery of a potential in the wave equation[END_REF].

Theorem A.1. [Carleman estimate for a wave operator without terms at t = 0] Let x 0 ∈ Ω ⊂ R d and assume the multiplier condition (1.4) for Γ 0 . Let m > 0, 0 < α 0 < α 1 and ρ ∈ (0, 1], and V as in (1.16). There exists a parameter ρ 0 = ρ 0 (Lm/α 0 ) > 0 such that for all β ∈ (0, α 0 ρ 0 ), there exist s 0 > 0 and a positive constant C such that for all s s 0 and for all a ∈ V,

s T -T Ω e 2sϕ (∂ t v) 2 + |∇v| 2 + s 2 v 2 dxdt + T -T Ω (P 1 (e sϕ v)) 2 + (P 2 (e sϕ v)) 2 dxdt C T -T Ω e 2sϕ ( a v) 2 dxdt + Cs T -T Γ0 e 2sϕ |∂ ν v| 2 dσdt + Cs {ϕ<0} e 2sϕ (∂ t v) 2 + |∇v| 2 + s 2 v 2 dxdt + Cs Ω e 2sϕ(±T ) (∂ t v(±T )) 2 + |∇v(±T )| 2 + s 2 v(±T ) 2 dx, (A.1)
where v ∈ L 2 ((-T, T );

H 1 0 (Ω)) satisfies a v = ∂ 2 t v -∇ • (a∇v) ∈ L 2 ((-T, T ) × Ω) and ∂ ν v ∈ L 2 ((-T, T ) × ∂Ω)
, and where ϕ denotes the weight function (1.18) and

{ϕ < 0} = {(t, x) ∈ (-T, T ) × Ω, ϕ(t, x) < 0}.
Besides, ρ 0 can be taken as in (1.14).

Proof. Using the weight ϕ defined by (1.18), we set, for s > 0, w(t, x) = v(t, x)e sϕ(t,x) for all (t, x) ∈ (-T, T ) × Ω. Then, we introduce the conjugate operator P defined by P w = e sϕ a (e -sϕ w). Some easy computations give P w = P 1 w + P 2 w + Rw, with

P 1 w = ∂ 2 t w -a∆w + s 2 ((∂ t ϕ) 2 -a|∇ϕ| 2 )w, P 2 w = -2s∂ t w∂ t ϕ + 2sa∇w • ∇ϕ + αsaw, Rw = -s(∂ 2 t ϕ -∇ • (a∇ϕ))w -∇a • ∇w -αsaw,
where α is a constant which will be chosen later, in (A.6). Since we have the main part of the proof is then to bound from below the cross-term in

P 1 w P 2 w dxdt a 2
by positive and dominant terms, similar to the one of the left hand side of (A.1), and a negative boundary term, that will be moved to the right hand side of the estimate. For the sake of clarity, we divide the proof in several steps. All the computations below are done for smooth velocities a and smooth functions v (equivalently, w). Then, by a classical density argument, we can extend the results to any a ∈ V and to any v satisfying a v ∈ L 2 ((-T, T ) × Ω) and ∂ ν v ∈ L 2 ((-T, T ) × ∂Ω).

Step 1. Explicit calculations

We set

P 1 w, P 2 w L 2 ((-T,T )×Ω,dt dx/a 2 ) = 3 i,k=1 I i,k
where I i,k is the L 2 scalar product of the ith-term in P 1 w by the kth-term in P 2 w with the metric dxdt/a 2 (x). We will mainly use integrations by parts and the boundary properties of w (such as w = 0 on (-T, T ) × ∂Ω). We shall also persistently use the fact that ∂ t ϕ is independent of x and ∇ϕ is independent of t, so that ∂ t ∇ϕ = ∂ 2 t ∇ϕ = 0 and ∂ t ∆ϕ = ∂ t (|∇ϕ| 2 ) = 0, and we emphasize that a depends only on the space variable and therefore ∂ t a = 0. Integrations by part in time give easily

I 11 = T -T Ω ∂ 2 t w (-2s∂ t w ∂ t ϕ) dxdt a 2 = -s T -T Ω ∂ t ((∂ t w) 2 )∂ t ϕ dxdt a 2 = s T -T Ω (∂ t w) 2 ∂ 2 t ϕ dxdt a 2 -s Ω (∂ t w) 2 ∂ t ϕ dx a 2 T -T
.

Similarly, one has

I 12 = T -T Ω ∂ 2 t w (2sa -1 ∇w • ∇ϕ) dxdt = -2s T -T Ω ∂ t w (∇∂ t w • (a -1 ∇ϕ)) dxdt + 2s Ω a -1 ∂ t w (∇w • ∇ϕ) dx T -T = -s T -T Ω ∇(∂ t w) 2 • (a -1 ∇ϕ) dxdt + 2s Ω a -1 ∂ t w (∇w • ∇ϕ) dx T -T = s T -T Ω (∂ t w) 2 ∇ • (a -1 ∇ϕ) dxdt + 2s Ω a -1 ∂ t w (∇w • ∇ϕ) dx T -T
, and

I 13 = T -T Ω ∂ 2 t w (αsa -1 w) dxdt = -s T -T Ω αa -1 (∂ t w) 2 dxdt + s Ω αa -1 ∂ t w w dx T -T .
We compute in the same way

I 21 = T -T Ω (-a -1 ∆w) (-2s∂ t w ∂ t ϕ) dxdt = -2s T -T Ω ∇w • ∇(a -1 ∂ t w ∂ t ϕ) dxdt = -s T -T Ω a -1 ∂ t (|∇w| 2 ) ∂ t ϕ dxdt + 2s T -T Ω ∇w • ∇a a 2 ∂ t w ∂ t ϕ dxdt = s T -T Ω a -1 |∇w| 2 ∂ 2 t ϕ dxdt -s Ω a -1 |∇w| 2 ∂ t ϕ dx T -T + 2s T -T Ω ∇w • ∇a a 2 ∂ t w ∂ t ϕ dxdt, and
I 23 = T -T Ω (-∆w)(αsw) dxdt = αs T -T Ω |∇w| 2 dxdt.
Using that w| (-T,T )×∂Ω = 0 implies ∇w = (∂ ν w)ν and |∇w| 2 = |∂ ν w| 2 on (-T, T ) × ∂Ω, we can compute (since ∇ 2 ϕ = 2I):

I 22 = T -T Ω (-∆w) (2s∇w • ∇ϕ) dxdt = 2s T -T Ω ∇w • ∇(∇w • ∇ϕ) dxdt -2s T -T ∂Ω (∇w • ν)(∇w • ∇ϕ)dσdt = s T -T Ω ∇(|∇w| 2 ) • ∇ϕ dxdt + 2s T -T Ω ∇ 2 ϕ(∇w, ∇w) dxdt -2s T -T ∂Ω |∂ ν w| 2 ∂ ν ϕdσdt = -s T -T Ω |∇w| 2 ∆ϕ dxdt + 4s T -T Ω |∇w| 2 dxdt -s T -T ∂Ω |∂ ν w| 2 ∂ ν ϕdσdt.
Finally, some integrations by part enable to obtain

I 31 = T -T Ω (s 2 ((∂ t ϕ) 2 -a|∇ϕ| 2 )w)(-2s∂ t w ∂ t ϕ) dxdt a 2 = -s 3 T -T Ω ((∂ t ϕ) 2 -a|∇ϕ| 2 )∂ t ϕ∂ t (w 2 ) dxdt a 2 = s 3 T -T Ω ∂ t ((∂ t ϕ) 2 -a|∇ϕ| 2 )∂ t ϕ w 2 dxdt a 2 -s 3 Ω ((∂ t ϕ) 2 -a|∇ϕ| 2 )∂ t ϕ w 2 dx a 2 T -T , I 32 = T -T Ω (s 2 ((∂ t ϕ) 2 -a|∇ϕ| 2 )w)(2sa -1 ∇w • ∇ϕ) dxdt = s 3 T -T Ω a -1 ((∂ t ϕ) 2 -a|∇ϕ| 2 )∇(w 2 ) • ∇ϕ dxdt = -s 3 T -T Ω ∇ • a -1 ∇ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) w 2 dxdt, and
I 33 = T -T Ω (s 2 ((∂ t ϕ) 2 -a|∇ϕ| 2 )w)(αsa -1 w) dxdt = αs 3 T -T Ω a -1 ((∂ t ϕ) 2 -a|∇ϕ| 2 )w 2 dxdt.
Gathering all the terms that have been computed, we get

T -T Ω P 1 wP 2 w dxdt = s T -T Ω (∂ t w) 2 a -2 ∂ 2 t ϕ + ∇ • (a -1 ∇ϕ)) -αa -1 dxdt +s T -T Ω |∇w| 2 a -1 ∂ 2 t ϕ -∆ϕ + 4 + α dxdt + 2s T -T Ω a -2 ∇w • ∇a ∂ t w ∂ t ϕ dxdt -s T -T ∂Ω |∂ ν w| 2 ∂ ν ϕdσdt +s 3 T -T Ω w 2 a -2 ∂ t ∂ t ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) + αa -1 ((∂ t ϕ) 2 -a|∇ϕ| 2 ) -∇ • a -1 ∇ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) dxdt -s Ω (∂ t w) 2 + a|∇w| 2 a -2 ∂ t ϕ dx T -T -s 3 Ω ((∂ t ϕ) 2 -a|∇ϕ| 2 )∂ t ϕ w 2 a -2 dx T -T +2s Ω a -1 ∂ t w (∇w • ∇ϕ) dx T -T + αs Ω a -1 ∂ t w w dx T -T . (A.3)
Step 2. Dominating terms

Step 2a. Terms involving derivatives • Focusing on the terms in s(∂ t w) 2 and s|∇w| 2 that we want to be strictly positive, we need to have, for all (t, x) ∈ (0, T ) × Ω, [START_REF] Imanuvilov | Global uniqueness and stability in determining coefficients of wave equations[END_REF], and recalling ϕ(t, x) = |x-x 0 | 2 -βt 2 , multiplying the first and second conditions respectively by a 2 and a, the two conditions will be satisfied if for all x ∈ Ω,

a -2 ∂ 2 t ϕ + ∇ • (a -1 ∇ϕ)) -αa -1 > 0 and (A.4) a -1 ∂ 2 t ϕ -∆ϕ + 4 + α > 0. (A.5) Using (x-x 0 )•∇a 2(1-ρ)a in Ω from a ∈ V defined in (1.
-2β + (2d -4 + 4ρ -α)a > 0 and -2β -(2d -4 -α)a > 0.
We choose to take

α = 2d -4 + 2ρ (A.6) so that -2β + (2d -4 + 4ρ -α)a = -2β -(2d -4 -α)a = 2aρ -2β.
Thus, for β ∈ (0, ρα 0 ), (A.7) the integrals of (A.3) involving |∂ t w| 2 and |∇w| 2 satisfy

s T -T Ω (∂ t w) 2 (a -2 ∂ 2 t ϕ + ∇ • (a -1 ∇ϕ)) -αa -1 ) dxdt + s T -T Ω |∇w| 2 (a -1 ∂ 2 t ϕ -∆ϕ + 4 + α) dxdt (A.8) 2s(ρα 0 -β) T -T Ω ((∂ t w) 2 + a|∇w| 2 )a -2 dxdt.
• The term of (A.3) involving products of ∂ t w and ∇w is the following one

2s T -T Ω a -2 ∇w∇a∂ t w ∂ t ϕ dxdt
We shall bound this term differently in {ϕ > 0} and in {ϕ < 0}.

In the set {ϕ < 0}, we bound this term as follows:

2s {ϕ<0} a -2 ∇w∇a∂ t w ∂ t ϕ dxdt Cs {ϕ<0} ((∂ t w) 2 + a|∇w| 2 ) dxdt, (A.9)
using that a ∈ W 1,∞ (Ω) with explicit bounds within the class V in (1.16).

In the set {ϕ > 0}, we use that βt

2 |x -x 0 | 2 , hence ∀(t, x) with ϕ(t, x) > 0, |∂ t ϕ| = 2β|t| 2 β sup Ω |x -x 0 | = 2 βL.
This implies

2s {ϕ>0} a -2 ∇w∇a∂ t w ∂ t ϕ dxdt 4s √ β √ α 0 L ∇a L ∞ {ϕ>0} |a -1 ∂ t w| |a -1/2 ∇w|) dxdt. 2s √ β √ α 0 L ∇a L ∞ {ϕ>0} ((∂ t w) 2 + a|∇w| 2 )a -2 dxdt. (A.10)
• Using (A.8) and (A.9)-(A.10), we obtain that if β > 0 is such that (A.7) holds and

√ β √ α 0 L ∇a L ∞ < ρα 0 -β, (A.11)
then there exist c * > 0 and C > 0 (depending in particular of a W 1,∞ ) such that for all s > 0,

s T -T Ω (∂ t w) 2 (a -2 ∂ 2 t ϕ + ∇ • (a -1 ∇ϕ)) -αa -1 ) dxdt + s T -T Ω |∇w| 2 (a -1 ∂ 2 t ϕ -∆ϕ + 4 + α) dxdt + 2s T -T Ω a -2 ∇w∇a∂ t w ∂ t ϕ dxdt (A.12) sc * T -T Ω ((∂ t w) 2 + |∇w| 2 ) dxdt -Cs {ϕ<0} ((∂ t w) 2 + a|∇w| 2 ) dxdt.
Step 2b. Terms in |w| 2 .

We write the term in s 3 |w| 2 as s 3 a -2 |w| 2 c w dxdt, where c w = c w (t, x) is given by

c w = ∂ t ∂ t ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) + αa((∂ t ϕ) 2 -a|∇ϕ| 2 ) -a 2 ∇ • a -1 ∇ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) . (A.13)
While it is clear that there exists a constant C depending in particular of a W 1,∞ such that

sup (t,x)∈(-T,T )×Ω |c w (t, x)| C,
we will show that this coefficient is actually positive in the set {ϕ > 0}. We first rewrite c w as:

c w = (∂ tt ϕ + αa -a∆ϕ)((∂ t ϕ) 2 -a|∇ϕ| 2 ) + 2∂ t ϕ 2 ∂ tt ϕ + 2a 2 D 2 ϕ(∇ϕ, ∇ϕ) + ∇ϕ • ∇a(∂ t ϕ) 2 = 2a 2 D 2 ϕ(∇ϕ, ∇ϕ) + 2∂ tt ϕa|∇ϕ| 2 + a|∇ϕ| 2 (-3∂ tt ϕ + a∆ϕ -αa) + (∂ t ϕ) 2 (3∂ tt ϕ -a∆ϕ + αa + ∇ϕ • ∇a) 16a 2 |x -x 0 | 2 -16βa|x -x 0 | 2 + 4a|x -x 0 | 2 (6β + a(2d -α)) + 4β 2 t 2 (-6β -a(2d -α)) -8β 2 t 2 |x -x 0 ||∇a(x)|.
From (A.6), α was chosen such that 2d -α = 4 -2ρ, implying

c w 16a|x -x 0 | 2 (a -β) + 4a|x -x 0 | 2 (6β + (4 -2ρ)a) -4β 2 t 2 (6β + (4 -2ρ)a)) -8β 2 t 2 |x -x 0 ||∇a(x)| 8|x -x 0 | 2 a 2 (4 -ρ) + 8a|x -x 0 | 2 β -4β 2 t 2 (6β + (4 -2ρ)a)) -8β 2 t 2 |x -x 0 ||∇a(x)|. Now, for (t, x) ∈ (-T, T )×Ω with ϕ(t, x) > 0, one has βt 2 |x-x 0 | 2 , so that for all (t, x) ∈ {ϕ > 0}, c w 8|x -x 0 | 2 a 2 (4 -ρ) + 8a|x -x 0 | 2 β -8β|x -x 0 | 2 (3β + (2 -ρ)a)) -8β|x -x 0 | 3 |∇a(x)| 8|x -x 0 | 2 a 2 (4 -ρ) + aβ -2βa + βρa -3β 2 -β|x -x 0 | ∇a L ∞ = 8|x -x 0 | 2 a [(4 -ρ)a -β(1 -ρ)] -3β 2 -β|x -x 0 | ∇a L ∞ . Assuming ρ 1 one can prove that a [(4 -ρ)a -β(1 -ρ)] α 0 [(4 -ρ)α 0 -β(1 -ρ)] . Indeed, (4 -ρ)α 0 -β(1 -ρ) 0 is equivalent to β α 0 4-ρ 1-ρ which is implied by condition (A.7) (0 < β < ρα 0 ) since ρ(1 -ρ) 4 -ρ is always satisfied. Therefore, c w 8|x -x 0 | 2 α 0 [(4 -ρ)α 0 -β(1 -ρ)] -3β 2 -βL ∇a L ∞
and if we impose that β satisfies (A.7) and

3β 2 + β (L ∇a L ∞ + α 0 (1 -ρ)) -α 2 0 (4 -ρ) < 0, (A.14)
which is satisfied for β small enough, then c w is positive on {ϕ > 0}.

• Summary: Provided β > 0 satisfies (A.7), (A.11) and (A.14), we get positive constants c * and C such that for all s > 0,

s T -T Ω (∂ t w) 2 (a -2 ∂ 2 t ϕ + ∇ • (a -1 ∇ϕ)) -αa -1 ) dxdt +s T -T Ω |∇w| 2 (a -1 ∂ 2 t ϕ -∆ϕ + 4 + α) dxdt + 2s T -T Ω a -2 ∇w∇a∂ t w ∂ t ϕ dxdt +s 3 T -T Ω w 2 a -2 ∂ t ∂ t ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) + αa -1 ((∂ t ϕ) 2 -a|∇ϕ| 2 ) -∇ • a -1 ∇ϕ((∂ t ϕ) 2 -a|∇ϕ| 2 ) dxdt sc * T -T Ω ((∂ t w) 2 + |∇w| 2 + s 2 |w| 2 ) dxdt -Cs {ϕ<0} ((∂ t w) 2 + a|∇w| 2 + s 2 |w| 2 ) dxdt.
From now on, we fix β satisfying (A.7), (A.11) and (A.14). Denoting r = L ∇a L ∞ /α 0 , it is not very difficult to check that this is equivalent to impose i.e. β ∈ (0, α 0 ρ 0 (L ∇a L ∞ /α 0 )) with ρ 0 = ρ 0 (r) given by

0 < β < α 0 ρ + r 2 4 - r 2
ρ 0 (r) = min    ρ + r 2 4 - r 2 2 , 1 6 (r + 1 -ρ) 2 + 12(4 -ρ) -(r + 1 -ρ)    . (A.16)
Tedious but easy computations show that the first term is in fact always smaller than the second one, so that we can simply take ρ 0 = ρ 0 (r) as in (1.14).

Step 3. Remaining right hand side terms |∂ ν w| 2 (x -x 0 ) • ν(x)dσdt.

Since we also have We take now s 0 large enough so that the second term of the second line is absorbed by the dominant term in s|∇w| 2 + s 3 w 2 as soon as s s 0 and we obtain Step 4. Back to the variable v Since w = ve sϕ and using P w = e sϕ a v, we can go back to the variable v in the previous Carleman estimate on w and obtain that (A.1) holds. This concludes the proof of Theorem A.1.

Proof of Theorem 1.4. This theorem is a corollary of the Carleman estimate of Theorem A.1. Indeed, let E s be the weighted energy of the waves defined by

E s [v](t) = 1
2 Ω e 2sϕ(t) (∂ t v(t)) 2 + a|∇v(t)| 2 dx.

As usual, C will be a positive constant independent of s that may vary from line to line.

On the one hand, recall first the weighted Poincaré inequality proved in [START_REF] Baudouin | Global carleman estimates for waves and applications[END_REF]:

Lemma A.2. Let φ ∈ C 2 (Ω) and assume that the weight φ defined on Ω satisfies inf Ω |∇ φ| δ > 0. Then there exist s 0 > 0 and C > 0 such that, for all s s 0 and for all z ∈ H 1 0 (Ω),

A.2 A Carleman estimate for a transport equation

Proof of Theorem 1.5. We denote by ϕ 0 the weight function ϕ 0 (x) = ϕ(0, x) = |x -x 0 | 2 where ϕ is defined in (1.18). We start with the following lemma: Proof of Lemma A.3. Without loss of generality, by replacing X by -X if needed, we assume that the vector field X satisfies ∀x ∈ Ω, X(x) • ∇ϕ 0 (x) -r 0 < 0. (A.20)

We multiply (A.18) by be 2sϕ0 and integrate over Ω: With X ∈ W 1,∞ (Ω) (in fact, ∇ • X ∈ L ∞ (Ω) is sufficient), choosing s 0 large enough, we obtain that there exists a constant C > 0 depending on r 0 and s 0 such that (A.19) holds for s s 0 .

To prove Theorem 1.5, we consider the following first order transport equation with source term g ∈ H We then focus on estimates of b in H 1 norm. In order to do that, we use that g ∈ H 1 X (Ω) and b ∈ H 1 X , so that the equation (A.21) yields X • ∇b = 0 on Γ X . On the other hand, b = 0 on Γ X , so ∇b = ∂ ν bν on Γ X . Hence we should have (X • ν)∂ ν b = 0 on Γ X . Since X • ν = 0 on Γ X , this implies that ∂ ν b = 0 on Γ X , and therefore ∇b = 0 on Γ X . Then one easily checks that for all i ∈ {1, • • • , d}, ∂ i b satisfies (A.18) with f i = ∂ i (g -b∇ • X). Applying Lemma A.3, we then get that for s s 0 In particular, summing over i ∈ {1, • • • , d} and using X ∈ W 1,∞ (Ω), there exists s 0 1 such that for all s s 0 , Using then X ∈ W 2,∞ (Ω), the last term can be absorbed by taking s sufficiently large, leading Theorem 1.5.

  e. a = a * | ∂Ω and a ν = ∂ ν a * | ∂Ω are known.(1.8)

Figure 2 :

 2 Figure 2: Definition and application of the cut-off function η, for T > L/ √ β.

. 15 ) 4 .

 154 Stage Conclusion using the first order Carleman estimate. The objective now is to estimate from below the left hand side of inequality (3.15) using the first order equation Carleman estimate presented in Theorem 1.5.The vector field used here is X = ∇w 0 ∈ W 2,∞ (Ω; R d ) and satisfies (1.10). We can thus apply estimate (1.20) to b = ãk+1 -a * , which by construction satisfies

Remark 4 . 5 .

 45 Compatibility condition.We consider the equations satisfied by u

Figure 4 :Figure 5

 45 Figure 4: Data used in the 1d numerical examples according to equation (4.4). (a) Initial data. (b) Boundary source term. (c) Distributed source term satisfying the compatibility condition (4.14) for a * (x) = 6 + sin(2πx). Example # of iterations Running time in sec err a Figure 5 (a) 4 405 2 • 10 -4 Figure 5 (b) 4 320 6 • 10 -4 Figure 5 (c) 4 264 1 • 10 -3 Figure 8 (a) 4 422 3 • 10 -3 Figure 8 (b) 4 498 5 • 10 -3 Figure 8 (c) 4 478 9 • 10 -3 Figure 10 3 1514 4 • 10 -3 Figure 11 4 1785 3 • 10 -3

2 Figure 5 :

 25 Figure 5: Top line: Convergence history of the reconstruction process. The blue line is the initial guess a 0 computed by the minimization of K init (4.3). In the case (c), it is equal to the constant coefficient a 0 = 6 and does not contain any information regarding the position of the bump. We notice that at the first iteration (orange line), the reconstruction is already very close to the exact solution. Bottom line: final reconstruction results for three different examples belonging to V * a,b .

  (a) Cut-off functions η i (ϕ(t, Lx)) plotted with respect to t < Lx

  Evolution of the corresponding solution a 0 i for 0 i 8.

Figure 6 :

 6 Figure 6: Illustration of the progressive process described in (4.15) on the first iteration of the recovery of the speed a * (x) = 6 + exp

Figure 8

 8 shows the results for the same a * (x) with different levels of noise in the measurements (θ = 5%, 10% and 20% in (4.2)).

  (a) without noise (b) with θ = 20% noise (c) after regularization

Figure 7 :Figure 8 :

 78 Figure 7: Examples of observations A * (t) at x = L x for the coefficient a * (x) = 6 + sin(2πx). We show the impact of adding noise following (4.2) and the result obtained after applying the low-pass filter.

  (a) T = Lx √ α 0 (b) T = 0.8 Lx √ α 0 (c) a * ∈ V * a,b

Figure 9 :

 9 Figure 9: Some examples of reconstructions using Algorithm 2 when one assumption of the convergence Theorem 1.6 is not satisfied. (a)(b) The observation time is smaller than the minimal time given in (1.13). (c) The exact coefficient a * ∈ V * a,b because it does not satisfy the condition (1.5) for any ρ > 0. The time T for this example is computed using (1.13) formula with ρ 0 = 0.01.

27

 27

Figure 10 : 4 ) 1

 1041 Figure 10: a * (x 1 , x 2 ) = 5 + exp(

Figure 11 :

 11 Figure 11: Reconstruction of a * (x 1 , x 2 ) = 6 + (cos(x 1 π) + 0.2 sin(5x 1 π)) cos( 2x2π 1.1 ). Top line: exact coefficient in 2d and 3d views. Bottom line: coefficient recovered numerically.

  -Rw) 2 dxdt a 2 , (A.2)

  1 -ρ) 2 + 12(4 -ρ) -(r + 1 -ρ) , (A.15)

(|∇w| 2 + s 2 w 2 )(∂ t w) 2 +

 222 (P w) 2 dxdt + C T -T Ω (|∇w| 2 + s 2 w 2 ) dxdt, using (A.2), we get s T -T Ω (∂ t w) 2 + |∇w| 2 + s 2 w 2 dxdt + T -T Ω (P 1 w) 2 + (P 2 w) |∇w| 2 + s 2 w 2 dxdt + Cs Ω (∂ t w(±T )) 2 + |∇w(±T )| 2 + s 2 w(±T ) 2 dx.

(∂ t w) 2 +

 2 |∇w| 2 + s 2 w 2 dxdt + T -T Ω (P 1 w) 2 + (P 2 w) 2 dxdt C T -T Ω (P w) 2 dxdt + Cs T -T Γ0 |∂ ν w| 2 dσdt + Cs {ϕ<0} (∂ t w) 2 + |∇w| 2 + s 2 w 2 dxdt +Cs Ω (∂ t w(±T )) 2 + |∇w(±T )| 2 + s 2 w(±T ) 2 dx.

s 2 Ωee

 2 2s φ|z| 2 dx C Ω e 2s φ|∇z| 2 dx.(A.17)Since a in Ω is larger than α 0 > 0 and since |∇ϕ(t)| = 2|x -x 0 | 2 √ ε > 0 and v(t) ∈ H 1 0 (Ω) for all t ∈ (-T, T ), we can deduce that Ω e 2sϕ(t) (∂ t v(t)) 2 + |∇v(t)| 2 + s 2 v 2 (t) dx CE s [v](t),∀t ∈ (-T, T ).On the other hand, defining a cut-off function χ ∈ C ∞ (R) such that 0 χ 1 andχ(τ ) = 1, if τ 0, χ(τ ) = 0, if τ -T, we can write E s [v](0) = χ(0)E s [v](0) -χ(-T )E s [v](-T ) = 0 -T d dt (χ(t)E s [v](t)) dt.Easy calculations (omitting the variable t in the integrals) gived dt E s [v](t) = s Ω e 2sϕ ∂ t ϕ (∂ t v) 2 + a|∇v| 2 dx + Ω e 2sϕ ∂ t v a v dx -2s Ω e 2sϕ a∂ t v∇v • ∇ϕ dx. Therefore, since χ ∈ C ∞ (R, [0, 1]), a < α 1 in Ω,and using Cauchy-Schwartz inequality, 2sϕ (∂ t v) 2 + |∇v| 2 dxdt + C T -T Ω e 2sϕ ( a v) 2 dxdt, the right hand side quantity being bounded by the Carleman estimate (A.1).Gathering these informations, one finally getsΩ e 2sϕ(0) (∂ t v(0)) 2 + |∇v(0)| 2 + s 2 v 2 (0) dx CE s [v](0) C T -T Ω e 2sϕ ( a v) 2 dxdt + Cs T -T Γ0 e 2sϕ |∂ ν v| 2 dσdt + Cs {ϕ<0} e 2sϕ (|∂ t v| 2 + |∇v| 2 + s 2 v 2 ) dxdt + Cs Ω e 2sϕ(±T ) (∂ t v(±T )) 2 + |∇v(±T )| 2 + s 2 v(±T ) 2 dx,that allows to conclude to (1.17) and ends the proof of Theorem 1.4.

Lemma A. 3 . 2 Ω|b| 2 e

 322 Let X ∈ W 1,∞ (Ω) be such that inf Ω |X • ∇ϕ 0 | > 0.Then there exist C > 0 and s 0 1 such that for any f ∈ L 2 (Ω), the solution b ofX • ∇b = f in Ω, b = 0, on Γ X , (A.18)(where Γ X is defined in (1.19)) satisfies, for all s s 0 ,s 2sϕ0 dx C Ω |f | 2 e 2sϕ0 dx.(A.[START_REF] Klibanov | Inverse problems and Carleman estimates[END_REF]

1 2 Ω 2 Ω 2 ∂Ω e 2sϕ0 b 2 XΩ e 2sϕ0 b 2 dx + 1 2 Ω e 2sϕ0 b 2 ∇ • X dx Ω e 2sϕ0 bf dx 1 2r 0 s Ω e 2sϕ0 f 2 dx + r 0 s 2 Ω e 2sϕ0 b 2

 122222222 e 2sϕ0 ∇(b 2 ) • X dx = Ω e 2sϕ0 bf dx.After an integration by parts we get-1 e 2sϕ0 b 2 ∇ • X dx -s Ω e 2sϕ0 b 2 ∇ϕ 0 • X dx + 1 • ν dx = Ω e 2sϕ0 bf dx.Since the boundary data in (A.18) yields that the boundary term is non-positive, and since X•∇ϕ 0 -r 0 in Ω, we get with Young's inequality r 0 s dx.

1 X 2 Ω |b| 2 e 2 Ω|b| 2 e

 12222 (Ω): ∇ • (bX) = g in Ω, b = 0, on Γ X . (A.21)Of course, this fits the setting of Lemma A.3 with f = g -b∇ • X, so that we have, for all s s 0 ,s 2sϕ0 dx C Ω |f | 2 e 2sϕ0 dx 2C Ω |g| 2 e 2sϕ0 dx + 2C Ω |b| 2 |∇ • X| 2 e 2sϕ0 dx.Choosing s 0 larger if needed, we easily get that for all s s 0 ,s 2sϕ0 dx C Ω |g| 2 e 2sϕ0 dx. (A.22)

s 2 ΩΩ|∇g| 2 e

 22 |∂ i b| 2 e 2sϕ0 dx C Ω |f i | 2 e 2sϕ0 dx 3C 2sϕ0 dx + 3C Ω |∇b| 2 |∇ • X| 2 e 2sϕ0 dx + 3C Ω |b| 2 |∇(∇ • X)| 2 e 2sϕ0 dx.

s 2 Ω 2 e 2 Ω |∇b| 2 e 2sϕ0 dx + s 4 Ω

 2224 |∇b| 2sϕ0 dx C Ω |∇g| 2 e 2sϕ0 dx + C Ω |b| 2 |∇(∇ • X)| 2 e 2sϕ0 dx. (A.23) Multiplying (A.22) by s 2 and summing to the last estimate, we obtain that for all s s 0 , s |b| 2 e 2sϕ0 dx C Ω |∇g| 2 e 2sϕ0 dx + Cs 2 Ω |g| 2 e 2sϕ0 dx + C Ω |b| 2 |∇(∇ • X)| 2 e 2sϕ0 dx. (A.24)

Table 1 :

 1 Numerical values of the variables used for all the numerical examples in the 1d setting according to equation (4.4).

	cos(πt)	see Remark 4.5

* ) max(a * )/0.95 given by (1.13) see

Remark 4.4

  • By Cauchy-Schwarz, all the terms in t = ±T can be easily bounded from below by w(±T)) 2 + |∇w(±T )| 2 + s 2 (w(±T )) 2 dx.• The boundary set Γ 0 is such that we can finally write(∂ t w) 2 + |∇w| 2 + s 2 w 2 dxdt -Cs {ϕ<0} (∂ t w) 2 + |∇w| 2 + s 2 w 2 dxdt w(±T )) 2 + |∇w(±T )| 2 + s 2 w(±T ) 2 dx

	T -T Ω	P 1 wP 2 w	dxdt a 2	sc *	T -T Ω
				-Cs
					T
				-s
					-T Γ0
		-Cs

Ω (∂ t Ω (∂ t

1-8|x-0.4| 2 1 1>8|x-0.4| 2 .

Acknowledgements.

Lucie Baudouin and Maya de Buhan have been partially supported by the Agence Nationale de la Recherche, Project ODISSE, grant ANR-19-CE48-0004-01. Sylvain Ervedoza has been partially supported by the Agence Nationale de la Recherche, Project IFSMACS, grant ANR-15-CE40-0010, and by the CIMI Labex, Toulouse, France, under grant ANR-11-LABX-0040-CIMI.