In general relativity one wants to find "harmonic" functions $u$ on (a neighborhood $U$ in) a Lorentzian manifold $(M, g)$. In arbitrary coordinates $\{x^\mu\}$ the equation $\Delta_g u = 0$ reads $$ g^{\mu\nu}(\partial_{\mu\nu}u - \Gamma^\alpha_{\mu\nu}\partial_\alpha u) = 0. \tag{$1$} $$ My question is how exactly to show existence/uniqueness of solutions via results for second order linear equations (equation (1) is such an equation). If $g$ were Riemannian, we could apply linear elliptic theory to this equation. However, even then, after specifying boundary conditions, I am not sure if we are guaranteed existence, since Lax-Milgram and the Fredholm alternative only give existence under certain conditions on the coercivity of the associated bilinear form.
However, for the Lorentzian metric, we want to get this into a hyperbolic form $\partial_{tt} u + L'u = 0$, where $L'$ is an elliptic operator. However, a priori (1) involves all sorts of mixed partials, and so I am unsure how to get (1) into this form, or if this is even the right approach.
For context: this is on the road to showing local existence of wave coordinates on a Lorentzian manifold. I would be interested in any resources on this topic. Many thanks!