Are P and S waves well-defined if we account for heterogeneity of Lamé parameters? Consider the isotropic elastic wave equation.
\begin{align}
\rho \ddot{\boldsymbol{u}} = \nabla \cdot \boldsymbol{\sigma} + \boldsymbol{f}
\end{align}
where
\begin{align}
\boldsymbol{\sigma}(\boldsymbol{x}) = \lambda(\boldsymbol{x}) (\nabla\cdot \boldsymbol{u}(\boldsymbol{x})) \boldsymbol{I} + \mu(\boldsymbol{x})\left( \nabla \boldsymbol{u}(\boldsymbol{x}) + (\nabla \boldsymbol{u}(\boldsymbol{x}))^{T}\right)
\end{align}
is the Cauchy stress tensor and $\boldsymbol{f}$ are internal body forces.
Expanding this equation out with the product rule gives
\begin{align}
\rho \ddot{\boldsymbol{u}} = \color{blue}{(\lambda + 2\mu) \nabla(\nabla \cdot \boldsymbol{u}) - \mu \nabla \times \nabla \times \boldsymbol{u}} + \color{red}{\nabla \lambda (\nabla \cdot \boldsymbol{u}) + \nabla \mu \cdot \left( \nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^{T}\right)}.
\end{align}
The standard derivation for P waves and S waves involves dropping the term in red, valid when the Lamé parameters are constant. This leads to a simplified elastic wave equation of
\begin{align}
\rho \ddot{\boldsymbol{u}} = (\lambda + 2\mu) \nabla(\nabla \cdot \boldsymbol{u}) - \mu \nabla \times \nabla \times \boldsymbol{u}
\end{align}
From there, use the Helmholtz decomposition to write
\begin{align}
\boldsymbol{u} = -\nabla (\nabla \cdot \boldsymbol{u}_{p}) + \nabla \times \boldsymbol{u}
_{s}.\end{align}
From there, you can show that $\nabla \cdot \boldsymbol{u}_{p}$ and $\nabla \times \boldsymbol{u}_{s}$ are solutions to an acoustic wave equation with velocities
\begin{align}
c_p &= \sqrt{\frac{\lambda + 2 \mu}{\rho}} \\
c_s &= \sqrt{\frac{\mu}{\rho}}
\end{align}
respectively. This derivation is done in two steps. First, take the divergence of both of sides of the simplified wave equation and note that curl of divergence is zero and divergence of gradient is the Laplacian. Similarly, take the curl of both sides to get a form for $\nabla \times \boldsymbol{u}_{s}$. If we include the red term in the full isotropic elastic wave equation (blue and red colored version), then this trick no longer works! Therefore, are P and S waves still defined for the case with fully heterogeneous  Lamé parameters?
Here is my approach. Simply rewrite the Lamé parameters in terms of what the P and S-wave velocities should be while assuming that the Helmholtz decomposition trick works to get orthogonal components. That is, we rewrite
\begin{align}
\lambda &= \rho (c_p^2 - 2 c_s^2)\\
\mu &= \rho c_{s}^2
\end{align}
to get a new elastic wave equation written as
\begin{align}
\rho \ddot{\boldsymbol{u}} &= \rho c_p^2 \nabla (\nabla \cdot \boldsymbol{u}) - \rho c_s^2 \nabla \times \nabla \times \boldsymbol{u} + \nabla(\rho(c_p^2-2c_s^2)) (\nabla \cdot \boldsymbol{u}) + \nabla(\rho c_p^2) \cdot \left(\nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^{T}\right).
\end{align}
Then just stop there, knowing that the terms "P wave" and "S wave" do not actually refer to orthogonal components to build the full solution $\boldsymbol{u}$. This at least gives us a way to qualitatively look at our solution in terms of what the "P velocity" and "S velocity" are.
Is this the most complete approach or is there a way to get something as nice as the orthogonal P-wave and S-wave decomposition for when we drop the heterogeneous Lamé terms?
 A: For plane wave fields ${\bf u} = (u,v,w)^T$ propagating along $x$ (i.e., we therefore set ${\bf u} = {\bf u}(x,t)$ without loss of generality), doing so is fairly easy:
$$
\begin{aligned}
\rho u_{tt} &= \left[(\lambda+2\mu) u_x\right]_x \\
\rho v_{tt} &= [\mu v_x]_x \\
\rho w_{tt} &= [\mu w_x]_x
\end{aligned}
$$
The first component describes a wave polarized along the propagation direction, i.e. a compression waves. The last two components describe waves polarized in the transverse direction, i.e. shear waves.
Next, let us write the Helmholtz decomposition of the displacement field ${\bf u} = {\bf u}_h + {\bf u}_d$ where ${\bf u}_h = \nabla\varphi$ is curl-free and ${\bf u}_d = \nabla\times {\bf a}$ is divergence-free. It follows that
$$
\rho (\ddot {\bf u}_h + \ddot {\bf u}_d) = \nabla (\lambda \nabla\cdot {\bf u}_h ) + \nabla\cdot(2\mu \nabla{\bf u}_h ) + \nabla\cdot(\mu (\nabla+\nabla^T){\bf u}_d )
$$
Independently on the fact that $\rho$, $\lambda$, $\mu$ are constants or not, we can still consider the special cases ${\bf u} = {\bf u}_h$ and ${\bf u} = {\bf u}_d$. However, I doubt that doing so would be of any use, in general. In fact, it seems difficult to decouple ${\bf u}_h$, ${\bf u}_d$ for an arbitrary displacement field $\bf u$ with both components.
