Higher order energy decay rates for damped wave equations with variable coefficients

Under appropriate assumptions the energy of wave equations with damping and variable coefficients $c(x)u_{tt}-\hbox{div}(b(x)\nabla u)+a(x)u_t =h(x)$ has been shown to decay. Determining the rate of decay for the higher order energies involving the $k$th order spatial and time derivatives has been an open problem with the exception of some sparse results obtained for $k=1,2,3$. We establish estimates that optimally relate the higher order energies with the first order energy by carefully analyzing the effects of linear damping. The results concern weighted (in time) and also pointwise (in time) energy decay estimates. We also obtain $L^\infty$ estimates for the solution $u$. As an application we compute explicit decay rates for all energies which involve the dimension $n$ and the bounds for the coefficients $a(x)$ and $b(x)$ in the case $c (x)=1$ and $h(x)=0.$


Introduction
We will study hyperbolic equations of the form (1.1) c(x)u tt − div(b(x)∇u) + a(x)u t = h(x, t), x ∈ R n , t > 0, where a, b, c, and h satisfy the assumptions listed below. This system is generally accepted as a model of wave propagation in a heterogeneous medium with friction (given by a(x)u t ) and source terms h(x, t). The coefficient c(x) accounts for variable mass density, while b(x) is responsible for temperature changes as the wave travels in space (see the derivation in [6]). More surprisingly, (1.1) has been considered as a model for heat propagation by Cattaneo [2], and independently by Vernotte [17]. In this setting the constitutive equation for the flux q as given by replaces the classical Fourier's law q = b(x)∇u. The Cattaneo-Vernotte equation (1.2) is usually considered as appropriate to describe unsteady heat conduction. For more on the applicability of this model see [3,8,9]. The finite speed of propagation feature that the solutions enjoy under this formulation renders this system as a more accurate way to study the real effects of heat diffusion. We will show later how these aspects regarding the physical interpretation of (1.1) tie in with the mathematical analysis of finding the decay rates, but we believe that this paper offers compelling evidence as to why the diffusion phenomenon is an essential feature for the damped wave equations.
The literature is replete with results on energy estimates for hyperbolic systems with linear damping, but they mainly concern the case when b(x) = c(x) = 1 and h(x, t) = 0. In [13] we considered wave equations with variable coefficients in the case c(x) = 1, h(x, t) = 0 for which we obtained explicit decay rates that depend on the dimension and the growth of the coefficients. The main goal of this work is to show how the energies are related to each other; more precisely, we deduce decay rates for the higher order energies in terms of the decay rate for the first energy E(t). This issue was addressed for example in [11] where the author shows that for c(x) = b(x) = 1, and h(x, t) = 0 the energy has the following decay rate while for the second order energy the decay rate improves so one has E(t; u t ) = 1 2 (||u tt (t)|| 2 + ||∇u t (t)|| 2 ) (1 + t) −2 .
The idea of using the diffusion effect in order to show a faster decay for the higher order energies is also present in [5] where the authors show that for the elasticity system with linear damping u t one has One expects, however, that the gain in the decay rate when moving from the first energy to the second should be t −2 . This is motivated by the fact that the long time behavior of the linearly damped wave equation with constant coefficients (1.3) u tt − ∆u + au t = 0 resembles the behavior of the corresponding parabolic equation as it was suggested in [10], [12], [16]. By estimating the higher order energies for the Gaussian as a solution of (1.4), one notices that E(t; ∂ k+1 t u) decays faster than E(t; ∂ k t u) by a factor of t −2 . This suggests that a similar phenomena should occur for (1.3). Indeed, our paper shows that this conjecture is correct, so in the absence of other effects caused by inhomogeneities, one basically has .. In this spirit, by matching the gain in decay t −2 of E(t; ∂ k+1 t u) versus E(t; ∂ k t u) obtained for diffusion problems we show the optimality of our decay results for higher order energies. One of the most important contributions of this paper is developing a methodology for determining the rates of decay for energies of all orders for a variety of hyperbolic systems with linear damping.
The paper is organized as follows. In the next section we present the assumptions for our results, after which we present the main results regarding the decay of the higher order energies involving time derivatives. The weighted decay of the spatial derivatives is obtained for orders k = 2 and k = 4 from which we deduce L ∞ estimates of the solution in dimension n = 3. We apply these results in the case c = 1, h = 0 in section 6. In the Appendix one can find the theorem and proof regarding the finite speed of propagation of solutions of (1.1) for variable b and c.

Main Assumptions
Let a, b, and c be smooth coefficients satisfying the conditions: where a i , b i , c i , i = 0, 1 are positive constants and the exponents α, β, γ satisfy: We consider initial data and a source h ∈ L ∞ ([0, ∞), L 2 (R n )). The three functions have compact support: where q which measures the increase in the support of u is given in the Appendix.

Weighted L 2 estimates for time derivatives
The strategy behind determining the rate of decay for E(t; ∂ k t u) in terms of the decay for E(t; u) is based on some simple facts. First, since ∂ k t u solves the equation (2.3) with the different inhomogeneity ∂ k t h, the higher order energy E(t; ∂ k t u) satisfies (almost) the same decay estimate as E(t; u) with possibly different constants (the actual decay will depend on ∂ k t h). The question is whether these decay rates can be improved similarly to the case a = b = c = 1 and h = 0 where the Fourier representation of u shows faster decay rates of norms involving higher derivatives.
Higher order spatial norms are estimated from since taking directly x-derivatives produces new terms and changes the dissipative nature of this equation. In other words, we will express the x-derivatives in terms of the t-derivatives. Let us recall the definition and energy identity A simple consequence is the following. ( Proof. The energy identity (3.1) implies Notice that |hu t | ≤ au 2 t /2 + a −1 h 2 /2. The result follows from integration on [T 0 , T ].
We can now proceed with higher-order norms. Let v = u t to simplify notations. Then In order to obtain integrals of E(t; v), we multiply the equation for v with 1 2 W (t)v where W is to be chosen later. Integration on R n yields We will choose W satisfying The support of the solution u at time t is given by supp u(·, t) = {x ∈ R n : |x| ≤ R + Ct Assume that W also satisfies where C 0 is a constant; later we give explicit conditions in terms of W, a and c.
Next multiply identity (3.5) by (1 + t) ν and integrate on [T 0 , T ]: Hence, Our assumptions on the quadratic form yield a simpler estimate: Choose T 0 sufficiently large, such that so by combining similar terms we obtain We can estimate the last integral as follows: for some 0 < ω < 1 and w 1 , w 2 > 0. This implies that which together with (3.9) simplifies (3.8) to the following: The growth assumptions from (3.10) used in this last inequality yield: By using the estimate from Proposition 3.1 with µ = ν − ω − 1 and (3.4) we derive the following estimate: The proof is complete provided we show: Existence of the weight function W . Let (3.14) where : (i) w 1 ≤ w 0 ≤ w 2 with w 1 , w 2 given by (3.10).
(ii) the exponent ω is chosen such that This choice for the weight W satisfies all the constraints as we show below: • The inequality (3.4) 1 . By (2.1) it suffices to show: but for sufficiently large times t this holds since ω > 2(α−γ) 2−β−γ by (3.15). • (3.4) 2 holds since for W given by (3.14) we have W tt ≥ 0, W t < 0 and the coefficients a and c are positive. • For the left inequality in (3.6) we first complete the square in v and v t , so we are left to show By using (2.1) and (3.14) we reduce the problem of proving the above inequality to showing: The assumption (3.15) on ω gives us that the above inequality is true for large times. We follow a similar approach for the right inequality in (3.6); we complete the square so we have: We need to show The constant C 0 can be chosen sufficiently large so that we have 4C 0 −2 > 0. By (2.1) and (2.2) in order to prove (3.18) it is enough to show , which is done exactly as above when proving (3.17).
Denote by θ = ν − ω. From (3.13) with a simple induction argument we obtain Proposition 3.2. Let a, b, c, and h be smooth coefficients which satisfy (2.1) and (2.2), and let θ > 0. Then the solution u of equation (1.1) with initial conditions (2.3), which satisfy (2.4), satisfies the following weighted energy estimate: The above arguments allow us to also obtain pointwise decay rates for the energy as stated in the following: Proof. In (3.1) written for ∂ k t u we apply Young's inequality with appropriate coefficients such that the damping term disappears. This yields Simplifying the above inequality and using (3.19) we obtain the desired estimate.
Let us go back to lower order derivatives and estimate which by the Cauchy inequality gives We obtain the following The explicit decay can be obtained by using (3.20) with k = 0 and k = 1; the decay rate of energy associated with the damping term is approximately −1 over the decay rate of total energy, in agreement with Proposition 3.1 on the "average" in time decay rates. 4. Weighted L 2 estimates for spatial derivatives Decay estimates of higher order spatial derivatives are important for studying nonlinear perturbations of equation (1.1). Such problems require more regular solutions and naturally lead to W k,p norms with k ≥ 0 and 2 < p ≤ ∞. Although the latter norms are expected to decay faster than H k norms, this is not easy to verify. The first difficulty is that L p norms are not convenient to estimate by the multiplier method when p = 2. Here the standard approach is to use interpolation or embedding, i.e., the Gagliardo-Nirenberg or Sobolev inequalities. We thus reduce the question to multiplier estimates with losses of derivatives. The variable coefficients present an additional difficulty as x-differentiation produces terms that change the dissipative form of (1.1). A simple solution is to express x-derivatives in terms of t-derivatives from the equation: where M u = a −1 div(b∇u). The diffusion phenomenon means M u ≈ u t , so secondorder spatial derivatives are related with the first-order time derivative. Similarly we can derive an identity for M 2 u. Below we give two sufficient conditions on a, b, and c to guarantee decay estimates of M u and M 2 u, respectively. Given λ 1 , λ 2 ∈ [0, 1], the coefficients satisfy Proposition 4.1. Let M u = a −1 div(b∇u) and assume that u is a solution of (1.1).

(ii) If (4.1) and (4.2) hold, then
Proof. (i) The first claim follows from (4.1) and (ii) To verify the second claim, we apply M (uv) = uM v + vM u + 2(b/a)∇u · ∇v. We have the chain of identities Using the expression for M u, we further obtain The square of M 2 u is bounded by the sum of squares times 7: Finally, we rewrite last estimate to match conditions (4.1) and (4.2): Thus, we have It is possible to derive expressions for M k u, k ≥ 3, and find conditions on a, b, and c which yield decay estimates of higher order spatial derivatives. However, considering k ≤ 2 is sufficient for most applications in R n when n ≤ 3.

Then the following weighted estimates hold for the solution u of equation (1.1) with initial conditions (2.3) which satisfy (2.4):
Proof. (i) Recall the estimate of u in Proposition 3.1 where µ is replaced by θ + 1: (1 + t) θ E(t; u)dt We also need Proposition 3.2 with k = 1: Since λ 1 ≤ 1, the claim follows from Proposition 4.1(i).
(ii) A simple adaptation of Proposition 3.1 to u t yields Applying Proposition 3.2 to the integral of E(t; u t ), we obtain It is clear from Proposition 4.1(ii) and λ 1 , λ 2 ≤ 1 that au 2 tt is the main term in the upper bound of M 2 u. We readily estimate the remaining terms by Proposition 3.2 with k = 2, 3.
Remark 4.3. It is straightforward to establish pointwise estimates of a 1/2 M u L 2 and a 1/2 M 2 u L 2 for large t. Such results will be presented elsewhere.

L ∞ estimates in dimension n = 3
The most important applications of Propositions 3.4 and 4.2 concern L ∞ decay estimates. We give an example in n = 3, although we can treat all n ≤ 7 due to the embedding H 4 (R n ) ⊂ L ∞ (R n ) for 4 > n/2. It is important to mention that the standard Sobolev estimate u|| L ∞ ( u|| L 2 + ∆u|| L 2 ) is not suitable, since the L ∞ norm is expected to decay faster than the L 2 norm.
Proposition 5.1. Assume that n = 3 and a, b are C 1 -functions. If for every sufficiently regular u, then Hence u L ∞ is bounded in terms of E(t; u) and E(t; u t ) whenever u is a solution of problem (1.1), (2.3).
Proof. We set a −1/2 div(b∇u) = f and multiply with u to obtain The well-known formula which is valid for compactly supported u, yields the estimate u 2 (x) |a 1/2 (y)f (y)u(y)| b(y)|x − y| dy + |u(y)||∇ ln b(y) · ∇u(y)| |x − y| dy. (5.1) It follows from the Cauchy inequality that where the second factor comes from the Hardy inequality The assumptions on a and b give |∇(a 1/2 (y)b −1 (y)u(y))| 2 dy b(y)|∇u(y)| 2 dy, so the final estimate becomes The second integral in (5.1) admits similar estimates. From the Cauchy and Hardy inequalities, we have It is easy to find more explicit conditions on a and b which will allow us to use Proposition 5.1. This part is about integration by parts and Hardy's inequality.
Corollary 5.2. Assume that a ∈ C 1 and b ∈ C 2 are such that and the matrix with entries is non-negative definite. Then the conditions on a and b in Proposition 5.1 hold; hence the L ∞ estimate also holds.
Proof. Conditions (i), (ii), and Hardy's inequality readily show that for sufficiently regular u. Thus, it remains to verify the condition about ∆u L 2 . Let a −1/2 div(b∇u) = f and multiply this equation with ∆u. Then b(∆u) 2 dx = a 1/2 f ∆u dx − (∇b · ∇u)∆u dx.
Integrating by parts in the second term on the right side, we obtain (∇b · ∇u)∆u dx = 1 2 We have from assumption (iii) that b(∆u) 2 dx ≤ a 1/2 f ∆u dx.
Now assumption (i) and Cauchy's inequality imply 6. Application in the homogeneous case c = 1, h = 0 In this section we will derive explicit decay estimates for (1.1) when c(x) = 1 and h(x, t) = 0, i.e. for as an application of the results proven here and in [13]. To state these results we introduce the following: Hypothesis A. Let a and b satisfy the growth conditions listed in (2.1). Then there exists a subsolution A(x) which satisfies and has the following properties: As we will see below in Theorem 6.1, the quantity µ gives the decay of the weighted L 2 norm of u, while µ + 1 gives the decay of the energy of u. In several cases one can construct explicit subsolutions A which satisfy (a1)-(a3), so in these cases one can compute the value of µ as given by the formula These cases are: (1) a, b are radial functions.
(3) b is radial and a is arbitrary.
(4) a, b are arbitrary (no radial symmetry), α + β < 2 and For a full discussion of these special situations see Section 7 in [13].
Theorem 6.1. Assume that a and b satisfy (2.1) with α, β such that Also, assume that Hypothesis A holds. Then for every δ > 0 the solution of (6.1) satisfies for all t ≥ 1. The constant C δ depends also on R, a, b, and n.
An immediate consequence of the above theorem which follows by (a1) is Under the assumptions of Theorem 6.1, we have These results in conjunction with the decay estimates presented in Sections 3-5 yield the following: Theorem 6.3. Assume that a, b satisfy (2.1) and (6.6). Then for every δ > 0 and T 0 > 0 sufficiently large the following inequalities hold: (i) Weighted (in time) L 2 energy decay: (ii) Pointwise (in time) energy decay: (6.10) (iii) Weighted kinetic energy decay (the energy associated with the damping): (iv) Decay of spatial derivatives (recall that M u = a −1 [div(b∇u)]): (1 + t) θ+1 a(M u) 2 dxdt (6.12) Remark 6.4. Note from the above decay estimate (6.10) that the k-th order energy has a polynomial decay of order µ + 1 + 2k − δ, exactly 2k units below the decay of the first order energy as given by (6.7). The damping term has the decay rate µ + 2 − δ, as it follows from (6.11), which is 2 units below the decay of the L 2 norm of u given by (6.8). Finally, the average decay rates of second and fourth order spatial derivatives resemble those of u t and u tt , respectively. Here (1 + t) λ2 accounts for possibly fast oscillations in b(x)/a(x) as |x| → ∞.
Appendix A. Bounds on the support of solutions Recall the non-homogeneous hyperbolic equation with damping where the coefficients a, b, and c satisfy conditions (2.1) in the introduction. The propagation speed is determined by the ratio c(x)/b(x). We assume the following: A better bound on the propagation speed is given in terms of q satisfying see Evans [4], but this equation is globally solvable only in special cases. Fortunately the corresponding inequality is sufficient for propagation speed estimates and solvable in most cases. It is also convenient to have bounds depending on |x|, so we look for radial solutions q. We consider h ∈ C((0, ∞), L 2 (R n )) and There exist a unique solution u ∈ C((0, ∞), H 1 (R n ))∩C 1 ((0, ∞), L 2 (R n )) whose support is described below.
Proof of Proposition A.1. We can assume that (u 0 , u 1 ) and h are more regular, so u ∈ C 1 ((0, ∞), H 1 (R n )) ∩ C 2 ((0, ∞), L 2 (R n )); the general case follows from a standard approximation argument. Multiplying (A.1) with u t , we obtain This identity will be integrated over the exterior of a suitable truncated cone. The outer normal to M (T ) at (x, t) is given by the (n + 1)-vector We integrate the energy identity over C(T ) and apply the divergence theorem: for all T > 0. The integral over M (T ) is non-negative, since for all u t and ∇u whenever (bq ′ ) 2 ≤ cb, or 0 < q ′ ≤ c/b. Thus, Notice that |hu t | ≤ au 2 t /2 + h 2 /(2a). Hence, 1 2 q(|x|)>T +q (R) c(x)u 2 T + b(x)|∇u| 2 dx The integrands on the right side are 0 due to our assumptions on u 0 , u 1 , and h. This shows that u(x, T ) = 0 if q(|x|) > T + q(R).
Proof of Corollary A.2. We already know that u(x, t) = 0 if q(|x|) > t + q(R), where q ′ (|x|) ≤ c(x)/b(x). It is sufficient to find q, such that where the constants are given in (2.1). Solving for q, we obtain for some q 0 > 0.