The Green's function is known for a single d'Alembertian in the presence of a delta-function source, i.e.
$\Box_x G(x^\sigma-y^\sigma)=\delta^{(4)}(x^\sigma - y^\sigma)$
and is given by
$G(x^\sigma - y^\sigma) = -\frac{1}{4\pi|\textbf{x}-\textbf{y}|}\delta\left[|\textbf{x}-\textbf{y}|-\left(x^0-y^0\right)\right]\theta\left(x^0-y^0\right)d^3 y$
but is it known for two or more d'Alembertians?
i.e. $\Box^2_x G(x^\sigma-y^\sigma)=\delta^{(4)}(x^\sigma - y^\sigma)$
or
$\Box^n_x G(x^\sigma-y^\sigma)=\delta^{(4)}(x^\sigma - y^\sigma)$
You have $$\Box^{n}_{x} G_{n}(x-y)=\delta^{(4)}(x - y)$$ By writing $$G_{n}(x-y)=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}\tilde{G}_{n}(q)e^{iq(x-y)}d^{4}q$$ and $$\delta^{(4)}(x - y)=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}e^{iq(x-y)}d^{4}q$$ we obtain $$(|k|^{2}-\omega^{2})^{n}\tilde{G}_{n}(q)=1$$ So $$\hat{G}_{n}(t-\tau, k)=\int_{\mathbb{R}}\frac{d\omega}{2\pi}\frac{e^{i\omega(t-\tau)}}{(k^{2}-\omega^{2})^{n}}$$ We let $k^{2}=\alpha$, so that $$\hat{G}_{n}(t-\tau, k)=\mathcal{P}\int_{\mathbb{R}}\frac{d\omega}{2\pi}\frac{e^{i\omega(t-\tau)}}{(\alpha-\omega^{2})^{n}}=\frac{1}{(-1)^{n}\Gamma(n)}\frac{d^{n-1}}{d\alpha^{n-1}}\mathcal{P}\int_{\mathbb{R}}\frac{d\omega}{2\pi}\frac{e^{i\omega(t-\tau)}}{\omega^{2}-\alpha}=$$ $$=-\frac{1}{2(-1)^{n}\Gamma(n)}\frac{d^{n-1}}{d\alpha^{n-1}}\frac{\sin(\sqrt{\alpha}(t-\tau))}{\sqrt{\alpha}}=-\frac{(t-\tau)^{2n-1}}{2^{n}(-1)^{n}\Gamma(n)}\Big(\frac{1}{k(t-\tau)}\frac{d}{d(k[t-\tau])}\Big)^{n-1}\frac{\sin(k(t-\tau))}{k(t-\tau)}=$$ $$=\frac{(t-\tau)^{n}}{2^{n-1}\Gamma(n)k^{n-1}}j_{n-1}(k(t-\tau))$$ Where $j_{n}(z)$ is the spherical Bessel function of the first kind and $n^{th}$ order. Now, $$G_{n}(t-\tau, r-r')=\int_{\mathbb{R}^{3}}\frac{d^{3}k}{(2\pi)^{d}}\hat{G}(t, k)e^{ik(r-r')}=$$ $$=\frac{1}{(2\pi)^{3}}\int_{0}^{\infty}k^{2}dk\int_{0}^{2\pi}d\varphi\int_{-1}^{1}d[\sin{\vartheta}]e^{i{k}|r-r'|\sin{\vartheta}}\hat{G}(t-\tau, k)=$$ $$=\frac{2}{(2\pi)^{2}}\int_{0}^{\infty}dkk^{2}j_{0}(k|r-r'|)\hat{G}(t-\tau, k)=$$ $$\frac{(t-\tau)^{n}}{2^{n-2}\Gamma(n)(2\pi)^{2}}\int_{0}^{\infty}dk\frac{1}{k^{n-3}}j_{0}(k|r-r'|)j_{n-1}(k(t-\tau))$$ With $$\int_{0}^{\infty}dk\frac{1}{k^{n-3}}j_{0}(k|r-r'|)j_{n-1}(k(t-\tau))=I(n, |r-r'|, t-\tau)$$ We have $$I(1, |r-r'|, t-\tau)=\int_{0}^{\infty}dkk^{2}j_{0}(k|r-r'|)j_{0}(k(t-\tau))=$$ $$=\frac{\pi}{2|r-r'|^{2}}\delta((t-\tau)-|r-r'|)$$ So $$G_{1}(t-\tau, r-r')=\frac{(t-\tau)}{4\pi|r-r'|^{2}}\delta(|r-r'|-(t-\tau))$$ Which is indeed a correct retarded Green's function (which is equivalent to yours in distribution sense apart from the sign convention (i've the basis of the oposite handidness)), i.e. consider the solution of $\Box_{x}F(t, r)=K(t, r)$, then $$F(t, r)=(G_{1}*K)(t, r)=\int_{\mathbb{R}^{3}}d^{3}r'\frac{K(r', t-|r-r'|)}{4\pi|r-r'|}$$ Now consider $$\frac{\partial{I(n, r, t)}}{\partial{t}}=\int_{0}^{\infty}dk\frac{1}{k^{n-3}}j_{0}(kr)\frac{\partial}{\partial{t}}j_{n-1}(kt)=$$ $$=\int_{0}^{\infty}dk\frac{1}{k^{n-4}}j_{0}(kr)j_{n-2}(kt)-n\int_{0}^{\infty}dk\frac{1}{k^{n-3}}j_{0}(kr)j_{n-1}(kt)=I(n-1, r, t)-nI(n, r, t)$$ Hence $$\frac{\partial{I(n, r, t)}}{\partial{t}}=I(n-1, r, t)-nI(n, r, t)$$ We let $$I(n, r, t)=\int_{\mathbb{R}}\frac{d\omega}{2\pi}\tilde{I}(n, r, \omega)e^{i\omega{t}}$$ Which gives $$\tilde{I}(n, r, \omega)=\frac{1}{i\omega+n}\tilde{I}(n-1, r, \omega)=$$ $$=\Big(\prod_{k=1}^{n-1}\frac{1}{i\omega+n-(k-1)}\Big)\tilde{I}(1, r, \omega)$$ With $$\tilde{I}(1, r, \omega)=\int_{\mathbb{R}}\frac{\pi}{2r^{2}}\delta(t-r)e^{-i\omega{t}}dt=\frac{\pi}{2r^{2}}e^{-i\omega{r}}$$ We have $$I(n, r, t)=\frac{1}{4r^{2}}\int_{\mathbb{R}}d\omega\Big(\prod_{k=1}^{n-1}\frac{1}{i\omega+n-(k-1)}\Big)e^{i\omega(t-r)}$$ Consider a new complex variable $z=i\omega$, the integrand has $n-1$ simple poles on the negative real axis. So, we consider an semi-circular contour with radius R $\gamma_{R}$ on the left half of the complex plane, as the integrand behaves like $\sim\frac{e^{i\Im{z}(t-r)}e^{-|\Re{z}|(t-r)}}{z^{n-1}}$ as $|z|\rightarrow\infty, \ \Re{z}<0$ we send $R$ to infinity and the integral becomes $$I(n, r, t)=\frac{1}{4ir^{2}}\oint_{\gamma_{\infty}}dz\Big(\prod_{k=1}^{n-1}\frac{1}{z+n-(k-1)}\Big)e^{z(t-r)}=\frac{\pi}{2r^{2}}\sum_{l=1}^{n-1}e^{(l-(n+1))(t-r)}\prod_{k=1, k\neq{l}}^{n-1}\frac{1}{l-k}=$$ $$=\frac{\pi{(-1)^{n+1}e^{-(n+1)(t-r)}}}{2r^{2}}\sum_{l=1}^{n-1}\frac{(-e^{(t-r)})^{l}}{\Gamma(l)\Gamma(n-l)}=\frac{\pi{(-1)^{n+2}e^{-n(t-r)}}(1-e^{(t-r)})^{n-2}}{2r^{2}(n-2)!}$$ hence, for $n\geq{2}$ we have $$G_{n}(t-\tau, r-r')=\frac{{(\tau-t)^{n}e^{-n((t-\tau)-|r-r'|)}}(1-e^{((t-\tau)-|r-r'|)})^{n-2}}{2^{n+1}\pi|r-r'|^{2}\Gamma(n-1)\Gamma(n)}$$