Stability for solutions of wave equations with C^{1,1} coefficients

We consider the stable dependence of solutions to wave equations on metrics in C^{1,1} class. The main result states that solutions depend uniformly continuously on the metric, when the Cauchy data is given in a range of Sobolev spaces. The proof is constructive and uses the wave packet approach to hyperbolic equations.


Introduction
We consider the wave equation Here A(x, D x ) = a ij (x)D x i D x j is a uniformly elliptic operator, satisfying a ij = a ji , and a ij (x)ξ i ξ j ≥ c|ξ| 2 for ξ ∈ R n . We assume that the functions a ij are in the space C 1,1 (R n ), with norm a C k−1,1 = |α|≤k ∂ α a L ∞ . The question investigated in this article is the stable dependence of the solution u on the metric (a ij ). Intuition for the problem can be obtained from the simplest possible case, namely the one-dimensional wave equation with constant sound speed c. Given f ∈ L 2 (R), the equation (∂ 2 t − c 2 ∂ 2 x )u(t, x) = 0, u(0) = f, ∂ t u(0) = 0 has the solution u(t, x) = 1 2 [f (x − ct) + f (x + ct)]. Thus, the solution at time t = 1 is obtained by translating f by c units in the positive and negative directions. Since translation of L 2 functions is a uniformly continuous operation, we see that u(1, · ) depends uniformly continuously in L 2 norm on the sound speed c. A stronger result may be obtained if the initial data is smoother: if f ∈ H 1 , then f ( · − c) − f ( · − c ′ ) L 2 ≤ ∇f L 2 |c − c ′ |, and u(1, · ) depends Lipschitz continuously in L 2 on the sound speed.
Our main result is the following theorem, which shows that uniformly continuous or Lipschitz dependence are valid also for general wave equations with C 1,1 metrics. In this introduction, we state the theorem only in the case where an initial velocity g is present. In Section 7 we will give the straightforward extensions to cases where driving terms F and initial positions f are present, and also where A(x, D x ) is replaced by a divergence form or Laplace-Beltrami operator.
If −1 ≤ α < 2, and if A = (a ij ), B = (b ij ) satisfy (1) and u A , u B are the corresponding solutions, then for any ε > 0 there is δ > 0 such that Further, if −1 ≤ α ≤ 1 and g ∈ H α+1 , then where C depends only on M and n.
Here Existence and uniqueness of solutions to linear wave equations is of course classical and can be established via energy estimates under quite general conditions. We refer to [8] for a comprehensive account and further references. A constructive method, valid in the setting of our theorem, for proving existence and uniqueness of weak solutions was introduced in [6] using a wave packet approach.
Stability estimates such as the ones in Theorem 1.1 are also classical, see [8] and references given there. In Section 7 we give an easy argument which uses just the existence and uniqueness of solutions, without any precise knowledge about the solutions.
The novelty is that our proof of Theorem 1.1 is constructive: an explicit expression for the solution is given, and the stability properties are deduced from that. The constructive method gives the same intuition to stability as in the n = 1 case, namely that stability for solutions should be the same as for translating functions. We use the wave packet approach introduced in [6], where the initial data is decomposed into wave packets, and an approximate solution to the equation is obtained by translating wave packets along the Hamilton flow. The stability is exactly governed by this translation.
The main motivation for this study comes from inverse problems in seismic imaging, where the analysis of solutions of wave equations has a key role. Many existing results, see [3] for a survey, assume a linearization about a smooth sound speed c 0 and use microlocal analysis and calculus of Fourier integral operators (FIOs).
There has been recent interest, see [4], in the practically more realistic case where c 0 is not smooth, and in this case few results are known. One reason is that there is no calculus of nonsmooth FIOs. However, solution operators for nonsmooth wave equations are understood quite well due to the wave packet approach of [6], and a more precise analysis of these operators is expected to lead to new results. This work is an attempt in this direction. The related work [4] discusses propagation of singularities (in terms of concentration of wave packets) and application of wave packets in numerical computations. We remark that wave packets are the same as curvelets [1], [2], which have been introduced in image processing as an efficient way of representing functions with singularities on smooth curves.
The proof of Theorem 1.1 is based on constructing an explicit solution operator for the nonsmooth wave equation, following the method of [6]. The idea is to localize the initial data to dyadic frequency shells, and write each localized piece as a superposition of wave packets at the given frequency. The action of the wave group on a wave packet is well approximated by translating the packet along the Hamilton flow, and this gives an approximate solution operator for frequency localized initial data. These are added up to obtain a parametrix for the full equation.
Actually, to handle the nonsmooth coefficients, at frequency level 2 k one truncates the coefficients to frequencies less than 2 k/2 and uses the Hamilton flow for the truncated metric. The error terms resulting from this will add up to a bounded operator, and this error can be iterated away by solving a Volterra equation.
The precise construction of solution operator will be a combination of methods in [6] and [7] and it does not involve new ideas. However, in view of the stability result, we need to give the construction in detail to see how the operator depends on the metric. The main outline is the same as in [6], which used a discrete wave packet frame. Some computations are simplified if one uses instead a continuous wave packet representation (i.e. the FBI transform, see [5]). This was used in [9], [10], [11] for wave packets based on the Gaussian. We will follow [7] which used instead wave packets compactly supported in frequency, a property which will keep the dyadic frequency annuli separated.
The crux of the stability proof is Lemma 6.4, which considers the stability of translating along Hamilton flow. Lipschitz stability for translation involves a loss of one derivative, and the main point in the proof is that there is a smooth deformation of two Hamilton flows obtained by deforming the corresponding metrics. For the full result, one needs to check that the several corrections required to obtain an exact solution operator, most importantly the Volterra iteration, do not affect the stability given by translation.
We first prove Lipschitz stability with a loss of one derivative, and the uniform continuity is an immediate consequence. Since there is a loss of one derivative arising from the flow, we can afford to lose one derivative in certain other estimates as well. In this way, we get stability in terms of the C 0,1 norm of the metric instead of C 1,1 norm.
The plan of the paper is as follows. Section 2 contains some basic facts about Hamilton flows, and Section 3 introduces the FBI transform. In Section 4 we outline the construction of the solution operator, and Section 5 contains the details. The stability result, Theorem 1.1, is proved in Section 6. Section 7 discusses variations of the stability result.
The gradient with respect to x or ξ is denoted by d x or d ξ , and D = 1 i d. Throughout the paper, M will be a large constant so that (1) and (2) are satisfied. We write a b if a ≤ Cb where C > 0 depends only on M and the dimension n. Also, a ∼ b means that a b and b a. We write L p = L p (R n ), similarly for the L 2 Sobolev spaces H α , the spaces C 0,1 and C 1,1 , and the Schwartz space S . The mixed norm spaces are denoted by Acknowledgements. Research partly supported by the Academy of Finland. Part of this research was carried out during visits at MSRI and at the University of Washington, and I wish to express my gratitude to these institutions for their hospitality and support. I would like to thank Hart Smith for generous advice, and Gunther Uhlmann for suggesting the problem.

Hamilton flow
We record for later use some elementary facts related to Hamilton flows. Let χ(ξ) be a smooth cutoff supported in the unit ball with χ = 1 for |ξ| ≤ 1/2. We define smooth approximations of the coefficients a ij by . Then a ij k is supported in frequency in {|ξ| ≤ 2 k/2 }, and satisfies derivative bounds Consider the Hamilton flow related to It will be useful to define this in terms of the functions . The Hamilton equations are given bẏ , so given an initial condition (x(0), ξ(0)) = (y, η) ∈ T * R n , the Hamilton equations have a solution (x(t), ξ(t)) depending smoothly on t, y, η at least for small time.
It is well known that the solution (x(t), ξ(t)) exists for all time. For, if it exists for t in an interval I = (−t 0 , t 0 ), then ξ(t) = 0 for t ∈ I, and h(t) = |ξ(t)| 2 satisfies by the homogeneity of ṗ , ω(t)) and ω(t) = ξ(t)/|ξ(t)|. Solving the equation gives h(t) = e R t and this argument shows that Thus (x(t), ξ(t)) stays in a compact subset of T * R n when t ∈ I, and we may extend the solution past the endpoints of I and to all time.

Wave packet representation
Let g be a real, even Schwarz function in R n with g L 2 = (2π) − n 2 , and assumeĝ is supported in the unit ball. For λ ≥ 1 and y, x, ξ ∈ R n define This is a wave packet at frequency level λ, centered in space at x and in frequency at ξ. The Fourier transform is given bŷ The FBI transform of a function f ∈ S (R n ) is given by Suppose λ ≥ 2 6 . Then iff is supported in 1 4 λ < |ξ| < λ, then T λ f vanishes unless 1 8 x,ξ ), the adjoint has the form It follows that T * λ T λ = I, and T λ f L 2 (R 2n x,ξ ) = f L 2 (R n ) . The following result, stating the L 2 boundedness of FBI transform type operators, is from [7, Lemma 3.1].
Lemma 3.1. Let g x,ξ be a family of Schwartz functions whose Schwartz seminorms are bounded uniformly in x and ξ. Then the operator T , defined for Schwartz functions by , and the adjoint T * given by The norms of T and T * are bounded by a sum of finitely many Schwartz seminorms of the g x,ξ .

Outline of construction of solution operator
We will outline the construction of a solution operator S(t) : g → u(t, · ) for the problem (3). More details are given in Section 5.
Let us start with the standard Littlewood-Paley frequency localization. Let χ(ξ) be a smooth cutoff supported in the unit ball with χ = 1 for |ξ| ≤ 1/2, and take β k (D) to be a Littlewood-Paley partition of unity with We let g k = β k (D)g, and for each k we want to find an approximate solution of ( This will be done by reducing matters to first order hyperbolic equations. Consider the pseudodifferential symbol k is positively homogeneous of degree 1 in ξ and smooth when ξ = 0, and satisfies for |ξ| = 1 . we see that p ± k (x, D)f will be supported in |η| ∼ 2 k iff is supported in |η| ∼ 2 k . We will also use the symbols We are ready to write down the first approximate solution operator for (3). It will be given byS where k 0 is a constant depending on M which will be chosen later, and where This implies that when one applies the wave operator D 2 t −A(x, D x ) toS(t), the resulting operator will be of order 0, which is one order lower than expected. Therefore we really have an approximate solution operator.
However, since ∂ tS (t)g| t=0 is only close to g but not equal to g in general, we need to correct the initial values of the operator. Thus instead ofS(t) we use . If k 0 is chosen large enough then K will be small on H α (R n ), and S(t) will be an operator of order −1 with S(t)g| t=0 = 0 and ∂ t S(t)g| t=0 = g. Thus S(t)g will be an approximate solution of (3) with correct Cauchy data.
It remains to show that one may iterate away the error and obtain an exact solution. To do this we seek a solution of (3) of the form Here S(t, s) = S(t − s) is the operator corresponding to S(t) but where the initial surface is {t = s} instead of {t = 0}. Then S(t, s)g| t=s = 0 and ∂ t S(t, s)g| t=s = g, and one has We obtain where T (t, s) = (D 2 t − A(x, D x )) S(t, s). As remarked above, we will show that T (t, s) is an operator of order 0, which in the present setting with a C 1,1 wave operator will mean that it is bounded on H α (R n ) for −1 ≤ α ≤ 2. Then u will be a solution provided that G(t, x) = V (T (t)g(x)), where G = V F solves the Volterra equation Thus the full solution operator for (3) will be We will now give the details.

Details of construction of solution operator
Let (a ij ) be a symmetric matrix of C 1,1 functions satisfying (1). We take M large enough so that (1) holds also for the truncated metrics (a ij k ), and we also assume (2). We will use all the notations in Section 4. Also, k 0 will be a sufficiently large integer, depending on M and n.
We start by noting that if the Cauchy data is localized near frequency 2 k , then the operators constructed in the preceding section preserve this localization.
Using this lemma, it will be enough to consider a fixed frequency and the result will follow by summing over dyadic annuli. Thus, let λ = 2 k , and write p = p ± k , U (t) = U ± k (t), etc. The main idea is that the wave evolution of a wave packet at frequency λ can be well approximated by transport along the Hamilton flow. The transport operator will be L x,ξ = L ± k , given by It is easy to see that U (t)F will solve the corresponding transport equation.
Let now f be a function localized near frequency λ. We write f as a superposition of wave packets at frequency λ, To get an approximate solution u to the half-wave equation (D t +P y )u(t, y) = 0 with u(0, y) = f (y), where P y = p(y, D y ), we transport the wave packets for time t along the Hamilton flow of p. Then u is given by u(t, y) = T λ f (x, ξ)g λ (y; χ t,0 (x, ξ)) dx dξ.
Using that χ t,0 is a symplectic map with inverse χ 0,t , u will be equal to To measure how far u is from an exact solution, we compute the last equality by integration by parts. Thus The following lemma, corresponding to [7, Lemma 3.2], will be crucial.
where g x,ξ is a family of functions whose Fourier transform is supported in a ball of radius 2 and the Schwartz seminorms of g x,ξ are uniformly bounded in x, ξ. In fact g Here m x,ξ satisfies symbol estimates uniform in x and ξ, Proof. We define g x,ξ by the relation (g x,ξ ) λ = (P y + L x,ξ )g λ . Recalling the formula (7) for g λ , it follows that since ξ · d ξ p = p by homogeneity. By (8), we have which shows that (g x,ξ )ˆ= 0 outside a ball of radius 2. From Taylor's formula we see that the term in brackets is equal to m x,ξ (z, ζ) given by (11).
Next we apply a second order wave operator D 2 t − P 2 y to the approximate solution operator, and show that there is a loss of one derivative.
Since u and (D t + P y )u are supported in |ξ| ∼ λ and p is of order 1, we have P y (D t + P y )u L 2 λ f L 2 by Lemma 5.4.
Using (10) and writing (P y + L x,ξ )g λ = (g x,ξ ) λ , we get Since the argument in Lemma 5.3 gives The Schwartz seminorms of g x,ξ are uniformly bounded in x, ξ, and (12) shows that the same applies to the Schwartz seminorms of m x,ξ (z, D z )g x,ξ . Lemma 3.1 gives It remains to study the symbol seminorms ofm x,ξ (z, ζ) = L x,ξ m x,ξ (z, ζ).
Since |ξ| ∼ λ a computation as in (13) shows and Lemma 3.1 implies This concludes the proof.
We move on to the full parametrix, where all the frequencies are added up. The errors resulting from truncation are handled as in [6,Theorem 4.5].
Proof. Since q ± k is order −1 one has E ± k (t) : H α → H α+1 , and since E ± k (t)g is localized near |ξ| ∼ 2 k the sum converges in H α+1 andS(t) : and write the last expression asT 1 (t)g +T 2 (t)g +T 3 (t)g. It is clear that which shows that (P + k ) 2 − A k (x, D x ) is of order 1 and Similar results hold for u − k . ForT 3 (t) we need to show that Hereβ k (ξ) are cutoffs to |ξ| ∼ 2 k which are 1 on the frequency support of u + k . We let l 0 = l 0 (M ) be an integer so that supp(β k (ξ)) ⊆ {2 −l 0 +k ≤ |ξ| ≤ 2 l 0 +k }.
By looking at the supports on the Fourier side, we get Using that a − a k L ∞ 2 −k we obtain It follows that Γ is bounded on H α with |α| ≤ 1. We compute and note that this is bounded H α+1 → H α if |α| ≤ 1, due to the C 1,1 regularity of a and the argument above. This shows boundedness of Γ for −1 ≤ α ≤ 2.
In the next lemma we correct the value of ∂ tS (t) at t = 0.

Proof.
One has R ± k (0)g L 2 2 −k g L 2 by Lemma 5.4. Also, R ± k β k (D) is of order −1 by looking at the symbol expansion of (P ± k Q ± k −I)β k (D) and using (16). The norms do not depend on any previous value of k 0 , and we may choose k 0 so large that K will have norm ≤ 1/2 on H α for α ∈ [−α 0 , α 0 ], for any α 0 > 0. Thus I + K will be invertible on these spaces and the norm of the inverse will be ≤ 2.
It now remains to show that one obtains a full solution operator from S(t) by introducing a correction by solving a Volterra equation. The solution operator for the Volterra equation is the following. with norm bounded by e CM , such that for any F (t, x) ∈ L ∞ t H α x , G = V F solves the equation It is easy to check that the series converges in L ∞ t H α x to a solution which satisfies the desired norm estimate.
As in the outline, we write S(t, s) = S(t−s) as the operator corresponding to S(t) but with initial surface {t = s}. If s, t ∈ [−M, M ] then S(t, s) has the same properties as S(t), except that S(t, s)g| t=s = 0 and ∂ t S(t, s)g| t=s = g. We let so T (t, s) and T (t) are bounded on H α when −1 ≤ α ≤ 2 by Lemma 5.6. We let V be the solution operator for the Volterra equation corresponding to T (t, s).
The argument in the outline shows that if then (D 2 t − A(x, D x ))S(t)g = 0 and S(t)g| t=0 = 0, ∂ t S(t)g| t=0 = g. One only needs to check that the time derivatives are justified, but this may be done as in [6] and will not be needed for stability considerations. This ends the construction of the solution operator.

Stability
We now proceed to prove the stability part of the result. Let A(x) = (a ij (x)) and B(x) = (b ij (x)) be two symmetric matrices satisfying (1), and take M so large that the truncated metrics (a ij k ) and (b ij k ) also satisfy (1). Also assume that t satisfies (2).
We write P A = P ± k for the operator at frequency λ = 2 k defined in terms of the metric A, and similarly Q A , R A ,R A (t) etc. The following operators depend norm continuously on the metric.
Proof. Consider h(x, ξ) = (λ n/2 χ(λ 1/2 · ) * [F (A k ( · , ξ))− F (B k ( · , ξ))])β λ (ξ), whereβ λ is a cutoff to |ξ| ∼ λ and F (t) = t 1/2 . We wish to show that This will show the estimate for P A − P B , and the estimate for Q A − Q B follows from a similar result with F (t) = t −1/2 . In ∂ α x ∂ β ξ h we let the x-derivatives hit the mollifier, which gives the desired growth. We may thus assume that α = 0. Each ξ-derivative hittingβ λ (ξ) gives λ −1 , so we only need to look at the case when the ξ-derivatives hit If P, Q are pseudodifferential operators the symbol of P Q is where the last terms are oscillatory integrals, and Note that |α + β − γ| ≥ N . We have for a suitableβ λ . Suppressing the cutoffs, this has the symbol , and the arguments above show that the corresponding operator has norm λ −1 A − B L ∞ on L 2 . It is easy to see that the terms with 0 < |α| < N have the same bound. Finally, if |γ| + |δ| ≤ 2n + 1 then Taking N large enough and using standard estimates for oscillatory integrals gives the L 2 bound for R A − R B .
Let g A x,ξ be the Schwartz functions in Lemma 5.3, defined in terms of the metric A. Proof. It is enough to show that This follows from the expression (11) for m x,ξ , the computation (13), and the estimate (17).
Remark. It is also true that the Schwartz seminorms of g A x,ξ − g B x,ξ are A − B C 1,1 , which follows by replacing (17) with the alternate estimate However, due to translation along Hamilton flow, there is a loss of one derivative in the stability estimate in any case. Therefore we can afford to lose one derivative in other estimates as well. This results in stability in terms of the C 0,1 norm of the metric instead of C 1,1 .
Proof. One has The result follows from Lemmas 3.1, 5.3, 6.1 and 6.2.
From the preceding results, and from the factorization we see that (I +K) −1 , the operator which corrects the initial values, depends norm continuously on the metric: Next we consider the stability of the flow operator. Here we will lose one derivative to get Lipschitz stability. Lemma 6.4. Iff vanishes unless |ξ| ∼ λ, then Proof. Fix t and write We let C r (x, ξ) = rA λ (x, ξ) + (1 − r)B λ (x, ξ), and let Φ r = χ Cr = (x r , ξ r ) be the flow corresponding to the metric C r . Then (Φ r ) r∈[0,1] is a smooth family of symplectic diffeomorphisms of T * R n , and we have Let h(s) = h(s, r, x, ξ) = (x r (s), ξ r (s)/λ) where (x r (0), ξ r (0)/λ) = (x, ξ/λ). Differentiating the Hamilton equations for (x r , ξ r /λ) with respect to r, and using Lemmas 3.1 and 5.3 imply the desired estimate.
We proceed to prove stability results for the operators where all the frequencies are summed up.
by (22) and Lemma 6.1. Writing C ± k,A − C ± k,B = (P ± k,A ) 2 − (P ± k,B ) 2 − (A k (x, D x ) − B k (x, D x )), and using the argument in the end of Lemma 6.1, we see that whenβ k (ξ) is a cutoff to |ξ| ∼ 2 k . Lemma 6.5 shows For the last term we write (T A,3 (t)−T B,3 (t))g = The discussion in Lemma 5.6 and (23) give the bound for the first term.
The second term is handled as in Lemma 5.6, except that we use One first gets the bound for −1 ≤ α ≤ 0 and then for |α| ≤ 1 by computing the derivative. The result follows.
We note that the last estimate holds for the operators T (t, s) with uniform constants when t, s ∈ [−M, M ]. The final estimate we need is for the Volterra solution operator. T A (s l+1 , s l+2 ) · · · T A (s j−1 , s j )F (s j , x) ds j · · · ds 1 Choose C = C(M ) such that for t, s ∈ [−M, M ], Using the H α estimate for each T B , the H α+1 → H α estimate for T A − T B , and then the H α+1 estimate for each T A gives There are j terms of the form I(t, x) at level j. It follows that x .