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

$$G_{ij}f(x)=\int_{\mathbb{R}^{4}}\frac{(x-y)_{i}(x-y)_{j}}{|x-y|_{a}^{6}}f(y)d\text{Vol}_{a}(y).$$

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

$$|G_{ij}f(x)-\int_{|y|\sim 2^{n}}\frac{(x-y)_{i}(x-y)_{j}}{|x-y|_{a}^{6}}f(y)d\text{Vol}_{a}(y)|\lesssim 2^{n(\delta+\tau)}\lesssim(1+|x|_{a})^{\delta+\tau},$$

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

$$\begin{split}&C\int_{|y-x|_{a}\gtrsim\ell(x),|y|\sim|x|}\frac{1}{|x-y|_{a}^{4}}\ell(y)^{\delta}|y|^{\tau}d\text{Vol}_{a}(y)\\ \leq&C(1+|x|_{a})^{\tau}\int_{|y-x|_{a}\gtrsim\ell(x)}\frac{1}{|x-y|^{4-\delta}}(\frac{\ell(y-x)}{|x-y|_{a}})^{\delta}d\text{Vol}_{a}(y)\\ \leq&C(1+|x|_{a})^{\tau}\int_{r>\ell(x)}r^{\delta-1}dr\int_{S^{3}}(\frac{\ell(y^{\prime})}{|y^{\prime}|})^{\delta}d\text{Area}_{S^{3}}(y^{\prime})\\ \leq&C(1+|x|_{a})^{\tau}\ell(x)^{\delta}\end{split}$$

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

$$|\nabla^{2}_{g_{a}}\Delta_{a}^{-1}f|_{g_{a}}\leq C\ell^{\delta}(|\vec{\mu}|_{a}+1)^{\tau}$$

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$:

$$|d\Delta_{a}^{-1}f|_{g_{a}}\leq C\ell^{\delta+1}(|\vec{\mu}|_{a}+1)^{\tau},$$

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

$$|\nabla^{2}_{g^{(2)}}\Delta_{a}^{-1}f|\leq\sum|\frac{\partial^{2}}{\partial\mu_{i}\partial\mu_{j}}\Delta_{a}^{-1}f||\nabla\mu_{i}||\nabla\mu_{j}|+\sum|\frac{\partial}{\partial\mu_{i}}\Delta_{a}^{-1}f||\nabla_{g^{(2)}}^{2}\mu_{i}|,$$

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. ∎

This follows from Lemma 2.6, Remark 2.7, and the fact that for the flat model $g_{\text{flat}}$ the Laplacian on $T^{2}$-invariant functions coincides with the base Laplacian $\Delta_{a}$. ∎

As in the proof of Lemma 2.18 we may assume $f$ has compact support to ensure a priori the well definition of $G_{\text{Taub}}f$. We use cutoff functions to decompose $f$ into a sum of functions $f_{n}$ supported on $\{n\lesssim\mu_{2}\lesssim n+1,\text{dist}_{g_{a}}(\cdot,\mathfrak{D}_{1})<C_{2}\}$ centred around points $x_{n}\in\mathfrak{D}$, with Hölder bound $\left\lVert f_{n}\right\rVert_{C^{k,\alpha}(B(x_{n},C_{2}))}\lesssim n^{\tau}$. At a fixed point $x$ bounded away from $\text{supp}(f_{n})$, the contribution $G_{\text{Taub}}f_{n}$ is estimated by $|G_{\text{Taub}}f_{n}|\lesssim n^{\tau}(|x-x_{n}|_{a}+1)^{\epsilon-2}$, where $\epsilon>0$ is any given small number (cf. Corollary 2.17 and notice the translational symmetry of $g_{\text{Taub}}$ along $\mathfrak{D}_{1}$). Elliptic bootstrap gives

$$\left\lVert G_{\text{Taub}}f_{n}\right\rVert_{C^{k+2,\alpha}_{0,0}(B(x,\ell(x)/10))}\lesssim n^{\tau}(|x-x_{n}|_{a}+1)^{\epsilon-2}.$$

Summing over all $n\in\mathbb{N}$,

$$\begin{split}&\left\lVert G_{\text{Taub}}f\right\rVert_{C^{k+2,\alpha}_{0,0}(B(x,\ell(x)/10))}\lesssim\sum n^{\tau}(|x-x_{n}|_{a}+1)^{\epsilon-2}\\ \lesssim&\int_{1}^{\infty}y^{\tau}(\ell(x)^{2}+|\mu_{2}(x)-y|^{2})^{\epsilon/2-1}dy\\ \lesssim&\begin{cases}(|x|_{a}+1)^{\tau}\ell(x)^{\epsilon-1}\quad-1<\tau<1-\epsilon,\\ (|x|_{a}+1)^{\epsilon-2},\quad\tau\leq-1.\end{cases}\end{split}$$

Thus

$$\begin{cases}\left\lVert G_{\text{Taub}}f\right\rVert_{C^{k+2,\alpha}_{-1+\epsilon,\tau}}\leq C,\quad-1<\tau<1-\epsilon,\\ \left\lVert G_{\text{Taub}}f\right\rVert_{C^{k+2,\alpha}_{0,-2+\epsilon}}\leq C,\quad\tau\leq-1.\end{cases}$$

which controls $\left\lVert\nabla_{\text{Taub}}^{2}G_{\text{Taub}}f\right\rVert_{C^{k,\alpha}_{\delta,\tau}}$ under the numerical conditions on weight exponents. ∎

The error $\Delta_{g^{(2)}}(\chi_{1}G_{\text{Taub}}f)-f$ comes from two sources: the deviation of the metric $g_{\text{Taub}}$ from $g^{(2)}$, and the cutoff error. The metric deviation error is estimated in Lemma 2.6 which we recall as $\left\lVert g_{\text{Taub}}-g^{(2)}\right\rVert_{C^{k,\alpha}_{0,-1}}\leq C.$ In particular for $|\vec{\mu}|_{a}>C_{3}$ and on the support of $\chi_{1}$, we have $\left\lVert g_{\text{Taub}}-g^{(2)}\right\rVert_{C^{k,\alpha}_{0,0}}\leq CC_{3}^{-1},$ so the metric deviation error is $O(C_{3}^{-1})$.

We turn to the cutoff error. By Lemma 2.20

$$\begin{cases}\left\lVert\chi_{1}G_{\text{Taub}}f\right\rVert_{C^{k+2,\alpha}_{-1+\epsilon,\tau}}\leq C,\quad-1<\tau<1-\epsilon,\\ \left\lVert\chi_{1}G_{\text{Taub}}f\right\rVert_{C^{k+2,\alpha}_{0,-2+\epsilon}}\leq C,\quad\tau\leq-1.\end{cases}$$

which implies

$$\begin{cases}\left\lVert\nabla^{2}_{\text{Taub}}(\chi_{1}G_{\text{Taub}}f)\right\rVert_{C^{k+2,\alpha}_{-3+\epsilon,\tau}}\leq C,\quad-1<\tau<1-\epsilon,\\ \left\lVert\nabla^{2}_{\text{Taub}}(\chi_{1}G_{\text{Taub}}f)\right\rVert_{C^{k+2,\alpha}_{-2,-2+\epsilon}}\leq C,\quad\tau\leq-1.\end{cases}$$

so in particular on $\text{supp}(d\chi_{1})\cap\{|\vec{\mu}|_{a}>2C_{3}\}$ where $\ell\sim C_{3}C_{1}^{-1}$, we have

$$\left\lVert\nabla^{2}_{\text{Taub}}(\chi_{1}G_{\text{Taub}}f)\right\rVert_{C^{k,\alpha}_{\delta,\tau}}\leq CC_{3}^{-\epsilon}.$$

By the support assumptions $f=0$ on $\text{supp}(d\chi_{1})\cap\{|\vec{\mu}|_{a}>2C_{3}\}$, hence the cutoff error is $O(C_{3}^{-\epsilon})$. Combining the two errors give the claim. ∎

The absolute value is estimated by Corollary 2.16:

$$|G_{g^{(2)}}f|\leq C(1+|\vec{\mu}|_{a})^{-2+\epsilon}.$$

The higher order estimate $\left\lVert\nabla^{2}_{g^{(2)}}G_{g^{(2)}}f\right\rVert_{C^{k,\alpha}_{-2,-2+\epsilon}(\mathbb{C}^{3})}\leq C$ follows by bootstrapping, which controls $C^{k,\alpha}_{\delta,\tau}$ norm for the given range of weight exponents $\delta,\tau$. ∎

The key technique is to construct a parametrix $P_{g^{(2)}}$ for the Green operator semi-explicitly, with precise control on its mapping properties.

Given a function $f$ with $\left\lVert f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}=1$, temporarily assumed to have sufficient decay at infinity, we will construct an approximate solution $u=P_{g^{(2)}}f$ to the Poisson equation as follows. Take a smooth cutoff function $\chi^{\prime}$

$$\chi^{\prime}=\begin{cases}1\quad&\text{dist}(\cdot,\mathfrak{D})\geq 2,\\ 0\quad&\text{dist}(\cdot,\mathfrak{D})\leq 1,\end{cases}$$

then $\chi^{\prime}f$ has norm $\left\lVert\chi^{\prime}f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\lesssim 1$ and is supported in $\{\text{dist}(\cdot,\mathfrak{D})\geq 1\}$. Applying Lemma 2.18, the function $u_{0}=\Delta_{a}^{-1}(\chi^{\prime}f)$ satisfies $\left\lVert\nabla^{2}_{g^{(2)}}u_{0}\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C$. By Lemma 2.19 we can choose $C_{2}\gg 1$ large enough independent of $f$ to ensure

$$\left\lVert\Delta_{g^{(2)}}u_{0}-f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3}\cap\{\text{dist}(\cdot,\mathfrak{D})>C_{2}/2\})}\ll 1.$$

Next we take smooth cutoff functions $\chi_{1}^{\prime},\chi_{2}^{\prime},\chi_{3}^{\prime}$ near $\mathfrak{D}_{1},\mathfrak{D}_{2},\mathfrak{D}_{3}$, such that

$$\chi_{1}^{\prime}=\begin{cases}1\quad&\text{dist}_{g_{a}}(\cdot,\mathfrak{D}_{1})\leq C_{2}/2\text{ and }|\vec{\mu}|_{a}>2C_{2},\\ 0\quad&\text{dist}_{g_{a}}(\cdot,\mathfrak{D}_{1})\geq C_{2}\text{ or }|\vec{\mu}|_{a}<C_{2}\end{cases}$$

and similarly with $\chi_{2}^{\prime},\chi_{3}^{\prime}$. The function

$$f_{1}=\chi_{1}^{\prime}(f-\Delta_{g^{(2)}}u_{0})=\chi_{1}^{\prime}(f-\Tr_{g^{(2)}}\nabla^{2}_{g^{(2)}}u_{0})$$

is supported in $\{\text{dist}_{g_{a}}(\cdot,\mathfrak{D}_{1})\leq C_{2},|\vec{\mu}|_{a}\geq C_{2}\}$ with bound $\left\lVert f_{1}\right\rVert_{C^{k,\alpha}_{\delta,\tau}}\leq C$. So we can apply Lemma 2.20 and Lemma 2.21 to find $u_{1}=\chi_{1}G_{\text{Taub}}f_{1}$ with bounds

$$\left\lVert\nabla^{2}_{g^{(2)}}u_{1}\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C,\quad\left\lVert\Delta_{{g^{(2)}}}u_{1}-f_{1}\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3}\cap\{|\vec{\mu}|_{a}>2C_{3}\})}\leq CC_{3}^{-\epsilon}\ll 1.$$

Completely analogous constructions are made near $\mathfrak{D}_{2}$ and $\mathfrak{D}_{3}$, where we obtain $u_{2},u_{3}$ with similar bounds.

Let $\chi_{4}^{\prime}$ be a smooth cutoff function

$$\chi_{4}^{\prime}=\begin{cases}1\quad&|\vec{\mu}|_{a}\leq 2C_{3},\\ 0\quad&|\vec{\mu}|_{a}\geq 4C_{3},\end{cases}$$

and define $f_{4}=\chi_{4}^{\prime}(f-\Delta_{g^{(2)}}(u_{0}+u_{1}+u_{2}+u_{3}))$, which is supported in the ball $\{|\vec{\mu}|_{a}\leq 4C_{3}\}$ and admits the bound $\left\lVert f_{4}\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C$. Then we can apply Lemma 2.22 to obtain $u_{4}=G_{g^{(2)}}f_{4}$ with bounds

$$\left\lVert\nabla^{2}_{g^{(2)}}u_{4}\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C,\quad\Delta_{{g^{(2)}}}u_{4}=f_{4}.$$

We set $P_{g^{(2)}}f=u=u_{0}+u_{1}+u_{2}+u_{3}+u_{4}$. The key point is that by construction

$$\left\lVert\nabla^{2}_{g^{(2)}}u\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C,\quad\left\lVert\Delta_{{g^{(2)}}}u-f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\ll 1,$$

namely $u$ is an approximate solution to the Poisson equation with bounds. A subtlety is that $\nabla^{2}_{g^{(2)}}u$ is fully controlled while $u$ is only controlled up to an additive constant. In any event, after extending the definition of the operator by removing the fast decay hypotheses on $f$, we have defined a bounded linear operator between weighted Hölder spaces of $T^{2}$-invariant functions and symmetric 2-tensors on $\mathbb{C}^{3}$

$$\nabla^{2}_{g^{(2)}}P_{g^{(2)}}:C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})\to C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3},\text{Sym}^{2}),$$

such that the operator $\Tr_{g^{(2)}}\nabla^{2}_{g^{(2)}}P_{g^{(2)}}$ is an approximation to the identity. Thus

$$\mathcal{R}^{\prime}=\nabla^{2}_{g^{(2)}}P_{g^{(2)}}(\Tr_{g^{(2)}}\nabla^{2}P_{g^{(2)}})^{-1}$$

is a bounded right inverse to $\Tr_{g^{(2)}}$. Composing with the projection to the type (1,1)-forms defines the operator

$$\mathcal{R}=\sqrt{-1}\partial\bar{\partial}P_{g^{(2)}}(\Tr_{g^{(2)}}\nabla^{2}P_{g^{(2)}})^{-1},$$

which takes value in closed (1,1)-forms and is a bounded inverse to $\Tr_{g^{(2)}}$. It is worth commenting that the same operators work for different exponents $\delta,\tau$.

It remains to relate $\mathcal{R}$ and $\mathcal{R}^{\prime}$ to the Green operator $G_{g^{(2)}}$ when $f$ has sufficient decay at infinity. The point is that for fast decay weights $\delta+\tau<-2$ and $\delta<-2$, the Hessian control $\left\lVert\nabla^{2}_{g^{(2)}}u\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}\leq C$ together with the a priori qualitative decay $u\to 0$ at infinity, imply the quantitative bound $\left\lVert u\right\rVert_{C^{k+2,\alpha}_{\delta+2,\tau}(\mathbb{C}^{3})}\leq C.$ This enables us to extend ${P}_{g^{(2)}}$ to a bounded linear operator

$${P}_{g^{(2)}}:C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})\to C^{k+2,\alpha}_{\delta+2,\tau}(\mathbb{C}^{3})$$

and the operator

$$P_{g^{(2)}}(\Tr_{g^{(2)}}\nabla^{2}P_{g^{(2)}})^{-1}:C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})\to C^{k+2,\alpha}_{\delta+2,\tau}(\mathbb{C}^{3})$$

defines an inverse to the Laplacian $\Delta_{g^{(2)}}$. By the uniqueness of decaying solution to the Poisson equation $G_{g^{(2)}}=P_{g^{(2)}}(\Tr_{g^{(2)}}\nabla^{2}P_{g^{(2)}})^{-1}$. Hence

$$\mathcal{R}^{\prime}=\nabla^{2}_{g^{(2)}}G_{g^{(2)}},\quad\mathcal{R}=\sqrt{-1}\partial\bar{\partial}G_{g^{(2)}}$$

as required. ∎

For $f$ with sufficient decay at infinity, we can find $u=P_{g^{(2)}}f$ with estimate $\left\lVert\nabla^{2}u\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3},\text{ Sym}^{2})}\leq C\left\lVert f\right\rVert_{C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})}.$ Using $\delta+\tau<-1$ and $\delta<-1$, we can integrate from spatial infinity to obtain the required gradient bound.

For a general $f$ without fast decay assumption, take a weakly convergent sequence of fast decaying functions $f_{k}\to f$ bounded in $C^{k,\alpha}_{\delta,\tau}(\mathbb{C}^{3})$, and find $u_{k}=P_{g^{(2)}}f_{k}$ with gradient bounds. After adjusting $u_{k}$ by additive constants to make $u_{k}(0)=0$, we can extract the subsequential limit $u$ of $u_{k}$, which solves $\Delta_{g^{(2)}}u=f$ with the gradient bound. ∎