For an equation that is actually hyperbolic, this is well-known. Here are some classical references:
Lars Garding, Cauchy's Problem for Hyperbolic Equations (1958)
Jean Leray, Hyperbolic Differential Equations (Institute of Advanced Study, 1953): see especially Chapter VI, section 4 as well as the extensions in Part III.
These are in addition to works of Friedrichs, Lewy, Petrowsky, and Sobolev; as well as treatments in Courant and Hilbert.
For general nonlinear second order hyperbolic systems derived from a Lagrangian formulation, there's also
- Demetri Christodoulou, The action principle and partial differential equations, (Princeton University Press, 2000): especially Section 5.3.
For a free resource, see
- Qian Wang's lecture notes; specifically section 2.2 which starts on slide 33. Note this statement applies to any quasilinear perturbation of the linear wave equation (and hence by a simple argument of Hormander applies to any fully nonlinear perturbation of the linear wave equation also). This can be generalized to the case where the background equation has variable coefficients too, and so by freezing coefficients you get also a uniqueness theorem.
Your equation, on the other hand, is not always hyperbolic. You can expand it to read
$$ u_{tt} - 3 u^2 \Delta u - 6 u |\nabla u|^2 = 0 $$
and when $u = 0$ the principal part of the equation is degenerate. In particular, it is not a perturbation of the linear wave equation.
So the standard energy methods cannot be applied directly.
In fact, for any constant $c \neq 0$, and any initial data $v_0, v_1\in C^\infty_c$, using the methods given in the linked resources above, you can show that if $u = c + v$ solves your nonlinear equation with initial data $u(0,x) = c + v_0(x)$ and $\partial_t u(0,x) = v_1(x)$, then $v$ has compact support.
Here's something interesting: consider the following three equations
\begin{gather}
u_{tt} = \Delta (u^3) \tag{NL1}\\
u_{tt} = \frac32 u \Delta (u^2) \tag{NL2} \\
u_{tt} = 3 u^2 \Delta u \tag{NL3}
\end{gather}
You see that all three equations can be written in the form
$$ u_{tt} = 3 u^2 \Delta u + \text{semilinear terms} $$
But for (NL1) and (NL3) you run into obvious difficulties when trying to prove energy estimates (like what you've seen already), but for (NL2) the following holds: multiplying the equation by $u_t$ you get
$$ u_{tt} u_t - \frac34 \partial_t (u^2) \Delta (u^2) = 0 $$
which implies
$$ \frac12 \partial_t (u_t)^2 + \frac38 \partial_t [\nabla (u^2)]^2 - \frac32 \nabla\cdot \left[ u u_t \nabla (u^2) \right] = 0 \tag{2} $$
Now, suppose $u$ is smooth and solves (NL2) with initial data $u(0,x)$ and $u_t(0,x)$ both smooth with compact support in the unit ball. Fix $\gamma > 0$. Consider the domain
$$ D_T := \{ t\in [0,T], |x| \geq \gamma t + 1 \} $$
and denote by
$$ \Sigma_T := \{T\} \times \{|x| \geq \gamma T + 1\} $$
and
$$ C_T := \{t\in [0,T]. |x| = \gamma t + 1\} $$
the boundaries. Integrating (2) on $D_T$ and applying the divergence theorem, and using the initial support hypothesis on $u$, gives us the energy estimate
$$ \int_{\Sigma_T} \frac12 (\partial_t u)^2 + \frac38 [\nabla (u^2)]^2 ~\mathrm{d}x - \int_{C_T} \gamma \left[ \frac12 (\partial_t u)^2 + \frac38 (\nabla u^2)^2 \right] + \frac32 u u_t \partial_r (u^2) ~\mathrm{d}\sigma
= 0 $$
where $~\mathrm{d}\sigma$ is some normalized surface measure on $C_T$. The key now is that by our assumption that $u$ is a smooth solution, since $C_T$ is compact, for sufficiently small $T$ it holds that $|u| \leq \gamma / \sqrt{3}$. (Here we use the initial support hypothesis.) And for such $T$ we have necessarily that the integrand in the $C_T$ integral is non-negative and coercive on $(\partial_t u, \nabla u)$. And so we conclude that $u \equiv 0$ on $D_T$.
A continuity argument then shows that provided $u$ remains a classical (in fact $C^1$ is enough) solution to (NL2) on the spacetime domain $[0,T]\times \mathbb{R}^d$, we must have that $u \equiv 0$ on the corresponding $D_T$. In other words, we can prove finite speed of propagation for (NL2) using a variation of the standard energy technique.
The same argument, however, fails for (NL1) and (NL3) due to a term that appears in the bulk integration. For those equations, the analogue of (2) has a nonzero term on the right hand side which we cannot Gronwall away.
This phenomenon is actually reminiscent of something that arises frequently in the study of linear hyperbolic systems. In Chapter 12 of Hormander's Analysis of Linear PDO you find the following statement:
Let $P_m$ be a homogeneous polynomial of degree $m$. Then the following are equivalent:
- $P_m$ is strictly hyperbolic
- For every polynomial $Q$ of degree at most $m-1$, the polynomial $P_m + Q$ is hyperbolic.
In the context of energy estimates, this roughly means that if the principal part of your equation is strictly hyperbolic, then you can perturb your equation by arbitrary lower order terms and still have energy estimates hold. But all bets are off when the principal part is not strictly hyperbolic. A simple example is the case of the biwave equation
$$ (- \partial_t^2 + c_1 \Delta)(-\partial_t^2 + c_2 \Delta) u = 0 $$
When $c_1$ and $c_2$ are distinct positive numbers, this equation can be perturbed by arbitrary third and lower order derivatives and one can prove energy estimates and finite speed of propagation more or less following the standard method. But when $c_1 = c_2$, you will find that there are certain lower order terms (even constant coefficient, linear) whose addition will break the energy estimates.