It's not true that all numerical relativity papers show $\Psi_4$; there are some groups directly producing $h = h_+ - i h_\times$ itself. But it is true that most numerical relativity groups only produce $\Psi_4$. Basically, there are two reasons:
- $\Psi_4$ is easier to measure than $h$ in practice.
- As far as GW detections are concerned, it doesn't really matter.
It's also not true that $\Psi_4$ is "more important" than $h$. While it's possible to compute radiated power from $\Psi_4$, that requires integrating $\Psi_4$, which means somehow coming up with a (complex) constant of integration. It's easier to compute power from $h$, because that just requires differentiation, and no constants of integration. IMO, it's better practice to measure $h$ from a simulation precisely because it contains all the information that $\Psi_4$, as well as two (complex) constants of integration.
The reason $h$ is not as easy to measure as $\Psi_4$ is because it is more sensitive to gauge effects (coordinate transformations). It may seem like you just have to take components of the metric tensor to get $h$, but if you look at discussions of gravitational waves in linearized gravity, you see that you have to jump through lots of hoops like Lorenz gauge and transverse-traceless gauge before you get to anything that's less sensitive to gauge effects. So you have to jump through similar hoops with numerical data, and it's not trivial to do that. On the other hand, $\Psi_4$ is constructed by just contracting the Riemann tensor with some "vectors", which is pretty quick to program.(*)
Now, as for why the difference doesn't really matter for GW detections, it helps to look in the frequency domain. If you know anything about LIGO, Virgo, or GW detectors in general, you probably know that they have a particular range of frequencies that they're sensitive to. In particular, pretty much no Earth-bound GW detector will be sensitive to frequencies much below 10Hz (or much above a few thousand Hz). So look at the frequency-domain representation of $h$. If we write the Fourier transform as
\begin{equation}
\tilde{h}(f) = \int_{-\infty}^{\infty}h(t)\, e^{-2\pi i f t}\, dt,
\end{equation}
it's not hard to see(**) that
\begin{equation}
\tilde{\ddot{h}}(f) = -4\pi^2 f^2 \tilde{h}(f).
\end{equation}
And the frequency domain is actually what gravitational-wave astronomers use to do their analyses, so this is what really matters. So if they want to use some numerical-relativity $\Psi_4$ "waveform", all they have to do to get whatever they would have gotten from $h$ is to divide:
\begin{equation}
\tilde{h}(f) = \frac{\tilde{\Psi}_4} {4\pi^2 f^2}.
\end{equation}
This works for every frequency except $f=0$ (because you can't divide by 0), but that doesn't matter anyway, because GW detectors are not sensitive to $0$Hz! To put this all another way, $\Psi_4$ contains all the information that $h$ does, except for two integration constants in the form of $a + bt$. But $a + bt$ is an extremely low-frequency signal, so GW detectors can't measure it.
Finally, as for why $\Psi_4$ is called a waveform, well that's just a very general term from signal processing that applies to basically any function of time.
(*)
Of course, $\Psi_4$ depends on how you pick those "vectors". In practice, those "vectors" depend on the coordinates you happen to be using, so they don't actually transform like vectors, which is why I used scare quotes. And in any case, you need to represent any waveform as a function of coordinates. So in the end, $\Psi_4$ still depends on basically arbitrarily coordinates, just as $h$ does. But it is easier to program.
(**)
As usual, "not hard to see" really means "a little subtle and not entirely obvious". The most obvious way to see that this equation is true is to integrate $\int_{-\infty}^{\infty} \ddot{h}(t)\, e^{-2\pi i f t}\, dt$ by parts twice, and throw in an assumption along the lines of $h(\pm \infty) \approx 0$ or something.
