The main task is to estimate the Calderon-Zygmund type operator

where $(x-y)_{i}$ denotes the components of $x-y$ viewed as a vector in $\mathbb{R}^{4}$. We say $x\in\mathbb{R}^{4}$ belongs to the dyadic scale $|x|\sim 2^{n}$ where $n\in\mathbb{N}$, if either $n=0$ and $|x|\leq 1$, or $n>0$ and $2^{n}\leq|x|\leq 2^{n+1}$. To ensure the Green operator is well defined, we will temporarily assume $f$ to be compactly supported, with no quantitative restriction on the measure of its support.

Since $\delta>-3$ and $|f(y)|\lesssim\ell(y)^{\delta}|\vec{y}|_{a}^{\tau}$, we have $\left\lVert f\right\rVert_{L^{1}(|y|\sim 2^{m})}\lesssim 2^{m(\delta+\tau+4)}$. Thus if $|x|\sim 2^{n}$ does not belong to scale $m$, then the contribution of $|y|\sim 2^{m}$ to $G_{ij}(x)$ is bounded by $O(2^{m(\delta+\tau+4)}\min\{2^{-4n},2^{-4m}\})$. Adding up all contributions from $m\neq n$, we get

using $\delta+\tau<0$ for the summability of the series. Since the source is far from the observer, elliptic bootstrap implies higher order Hölder regularity.

Now that we are left with only one scale, it is clear that the claimed bound for second derivatives holds where $\ell$ is comparable to $|\vec{\mu}|_{a}$. We now focus on $x$ close to $\mathfrak{D}$. The contribution of $|y-x|_{a}\gtrsim\ell(x),|y|\sim|x|$ is estimated by

where we use $-3<\delta<0$ in the convergence of the integrals. Since the contribution comes from sources at distance at least $\ell(x)$ away from the observer, the higher Hölder norms are controlled. Finally, the estimates for the contribution from $|y-x|_{a}\lesssim\ell(x)$ follows simply from standard Schauder theory.

At this stage we have proved the second derivative bound

together with an implicit weighted $C^{k,\alpha}$-bound in the $g_{a}$-metric. Since $f$ is compactly supported by our temporary assumption, qualitatively $\Delta^{-1}_{a}f$ has quadratic decay at infinity. By integrating the second order derivatives from infinity, we can bound first order derivatives $d\Delta_{a}^{-1}f$:

using $\delta+\tau<-1$ and $\delta<-1$ in the integration. Alternatively the first order derivative bounds can be proved using the same singular integral operator method.

Now the Hessian $\nabla^{2}_{g^{(2)}}\Delta_{a}^{-1}f$ can be expanded as a linear combination of second derivatives $\frac{\partial^{2}}{\partial\mu_{i}\partial\mu_{j}}\Delta_{a}^{-1}f$ etc and first derivatives $\frac{\partial}{\partial\mu_{i}}\Delta_{a}^{-1}f$ etc. Hence

where the sum includes also $\eta$-derivatives. Higher order derivatives of the Hessian can be expanded by the Leibniz rule. Using $\left\lVert d\mu_{i}\right\rVert_{C^{k,\alpha}_{0,0}}\leq C$ and $\left\lVert d\eta\right\rVert_{C^{k,\alpha}_{0,0}}\leq C$, we obtain the Hessian bound $\left\lVert\nabla^{2}_{g^{(2)}}\Delta_{a}^{-1}f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C$ as claimed.

Finally, an approximation argument in the weak topology removes the compact support assumption on $f$, so we conclude that $\nabla^{2}_{g^{(2)}}\Delta_{a}^{-1}$ extends canonically to a bounded linear operator between the weighted Hölder spaces.

As a delicate side remark, to bound the integral operator $\Delta_{a}^{-1}$ itself we would need to impose further $\delta+\tau<-2$ and $\delta<-2$, which would not be adequate for our intended applications. It is crucial in the above argument that the integral kernel of $G_{ij}$ decays two orders faster than the Green kernel. ∎