Generally you would Fourier transform the differential equation. The elements of the D'Alembertian transform to $\vec{k}\cdot\vec{k}$ (the Laplacian) and $-\omega^2$ (the 2nd time derivative), while the delta function transforms to $1$.
Then you get something like
$$ (2\pi)^4 \left(\left\lVert \vec{k}\right\lVert^2 - \frac{\omega^2}{c^2}\right) g= 4\pi $$
Solving for $g$ gives you the answer.
This is obviously greatly simplified and makes many assumptions about certain integrals converging, but it's the general approach.
I'm not sure in what generality and rigor you are expecting an answer; you may want to look at the theory of pseudodifferential operators on vector bundles. But perhaps the following will be more relevant to you:
It helps to think of $ P $ as an operator with matrix coefficients or equally as a matrix of differential operators.
The identity operator $ \mathcal{I} $ defined by $ \mathcal{I}[v] = v $ has the kernel $ I \; \delta(r) $ where $ I: \mathbb{R}^3 \rightarrow \mathbb{R}^3 $ is the identity matrix,
$$
\mathcal{I}[v](r) = \int \delta(r-r') \; I \cdot v(r') \; d^3r' = \int \delta(r-r') \; I \cdot v(r') \; d^3r' = \int \delta(r-r') \; v(r') \; d^3r' = v(r).
$$
Since $ I $ is a constant, it just factors out of the fourier transform, so that transform of the kernel of $\mathcal{I}$ is just $ I \; e^{i k \cdot r } $.
Let's look at those examples you had. The vector laplacian is
$$
\begin{pmatrix} \Delta & 0 & 0 \\ 0 & \Delta & 0 \\ 0 & 0 & \Delta \end{pmatrix},
$$ and in the fourier basis $$ ( k_x^2 + k_y^2 + k_z^2 ) \; I, $$ so it is easy to invert, and after inverse fourier transform will give you the usual Green's function.
Curl is
$$
\begin{pmatrix} 0 & -\partial_z & \partial_y \\ \partial_z & 0 & - \partial_x \\ -\partial_y & \partial_x & 0\end{pmatrix}.
$$
In fourier basis
$$
C = i \begin{pmatrix} 0 & -k_z & k_y \\ k_z & 0 & -k_x \\ - k_y & k_x & 0 \end{pmatrix}
$$
This one is tricky to invert because it has a kernel (gradients) and cokernel (div of curl is zero). You can check though that
$$
C \cdot \left( \frac{1}{k_x^2 + k_y^2 + k_z^2} C \right)
$$
is the identity when restricted to subspace $ \left\{ v \; \big| \; v \cdot k = 0 \right\} $, which is precisely the subspace of divergence free fields. Inverse fourier transform to obtain
$$
\frac{1}{(x^2+ y^2 + z^2)^{3/2}} \begin{pmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{pmatrix},
$$
which is a weird way of writing the kernel
$$
\frac{(r - r')}{|r-r'|^3} \times \cdot
$$
that you will recognize from the Biot-Savart law.
(Sorry for sloppy notation, missing $ -1 $s and $ \pi $s, etc.)
Best Answer
Well, $$\Delta G=\delta$$ by definition, so taking the Fourier transform of this gives that $$\hat{G}(\xi)=-(2\pi)^{-3/2}|\xi|^{-2}\in\mathcal{S}'.$$
In fact, $$\hat{G}(\xi)=-(2\pi)^{-n/2}|\xi|^{-2}$$ on $\mathbb{R}^n,$ and this is in $\mathcal{S}'$ for all $n\geq 3$, as $|\xi|^{-2}\in L^1_{loc}$ for such $n$.
Here, I used the convention that $$\hat{u}(\xi)=(2\pi)^{-n/2}\int\limits_{\mathbb{R}^n} u(x)e^{-ix\xi}\, dx.$$