I'm trying to understand characteristics of second order hyperbolic equations. When the number of dimensions is 2, say $$ au_{xx} + 2b u_{xy} + cu_{yy} = 0, \tag{1} $$ we have a change of variables $(x,y) \mapsto (\xi, \eta)$ that satisfies very particular conditions. In particular we reduce (1) to $$ w_{xy} = 0 \tag{2} $$ which can be solved. In doing so we learn a wealth of information about the characteristics of this equation, and the solution's behavior along these characteristics. My question is: is there a generalization of this technique to higher dimensions? For higher dimensions there are more than one mixed partial term, which is where I'm getting stuck.
The motivation for this is understanding wave propagation in 3 spatial (and 1 time) dimensions.
Edit: For an equation of the form $g^{ij}(x)\partial_{ij}u = 0$, with $g_{ij}$ a Lorentzian metric, I think this is a bit easier (theoretically, at least). Since it's a Lorentzian metric, we can pick "null coordinates" $\xi$ so that, say, in $n$-dimensions, $\xi^1$ and $\xi^2$ are null coordinates. Then along the null hypersurfaces $\{\xi^1 = a\}$ and $\{\xi^2 = b\}$ the solution's "outgoing" derivative (i.e. the outgoing derivative to $\{\xi^1 = a\}$ is $\partial_{\xi^1} u$) must satisfy a first order PDE. In a manner similar to the transport equation, I think this gives "propagation equations" for the outgoing derivative along the null (characteristic) hypersurfaces. See Rendall's paper (1990) on the characteristic initial value problem.