I would like to show the following function is $\gamma$-Hölder continuous. Said function $F:\mathbb{R}^n \rightarrow \mathbb{R}$ is defined by a singular integral operator of convolution type as follows (for some fixed $i$ and $j$)
$$ F(x)= \lim_{\epsilon \rightarrow 0} \int_{|x-y| \geq \epsilon} \partial_k \partial_j k_i(x-y) \left(Y_k(x)- Y_k(y) \right)f(y)dy $$
where $f \in C_c^{\gamma}(\mathbb{R}^n)$ (ie a compactly supported $\gamma$-Hölder continuous function) and the integral kernel $k:\mathbb{R}^n\backslash \lbrace \vec{0} \rbrace \rightarrow \mathbb{R}^n$ is defined as
\begin{equation} k_i(x) =\frac{x_i}{|x|^n} \end{equation}
Note $k$ is the derivative of the Newtonian kernel up to a constant factor. In addition, we have $Y: \mathbb{R}^n \rightarrow \mathbb{R}^n$ such that
$$ \| \nabla Y \|_{1, \gamma} = |Y(0)|+\|\nabla Y\|_0 + |\nabla Y|_{1,\gamma} < \infty $$ where $$ |W|_{\gamma} = \sup_{x \neq y} \frac{|W(x)-W(y)|}{|x-y|^{\gamma}} $$ is the $\gamma$-Hölder seminorm. Note $Y$ has bounded $\gamma$-Hölder continuous derivatives.
-----------------------------Attempt at Solution -------------------------------------
In attempting to prove the $\gamma$-Holder continuity of $F$, my strategy thus far has involved expanding $Y$ as a first order Taylor polynomial as follows $$ Y_k(y)=Y_k(x) + \partial_l Y_k(x)(y_l-x_l) +(y_l-x_l)\int_0^1 \left( \partial_l Y_k(x+t(y-x)) -\partial_l Y_k(x) \right)dt $$ I'll denote the remainder term as $R_x(y)$. Note $|R_{x}(y)| \leq C|y-x|^{1+\gamma}$.
So I write \begin{eqnarray} && F(x+h)-F(x)=\lim_{\epsilon \rightarrow 0} \bigg( \int_{|x+h-y|\geq \epsilon} \partial_k \partial_j k_i(x+h-y)\partial_lY(x+h)(x_l+h_l-y_l)f(y)dy \\ &&- \ \int_{|x-y|\geq \epsilon} \partial_k \partial_j k_i(x-y)\partial_lY(x)(x_l+h_l-y_l)f(y)dy -\int_{|x+h-y|\geq \epsilon} \partial_k \partial_j k_i(x+h-y)\partial_lR_{x+h}(y)f(y)dy \\ && + \int_{|x-y|\geq \epsilon} \partial_k \partial_j k_i(x-y)\partial_lR_{x}(y)f(y)dy \bigg) \end{eqnarray}
I can show the absolute value of the difference of the first two terms is bounded by $C|h|^{\gamma}$. This is done by first observing that $\partial_j k_i$ is the integral kernel of a singular integral operator (ie it is smooth everywhere but $\vec{0}$, is homogeneous of degree $-n$, and has mean value on the unit sphere). In fact, this mean-value property extends as follows:
\begin{equation} \int_{|x|=1}x^{\beta}D^{\alpha}k_i(x)dS=0 \quad \Rightarrow \int_{r \leq |x|\leq R} x^{\beta}D^{\alpha}k_i(x) dx = 0 \end{equation}
for any multi-indices $\alpha$ and $\beta$ such that $|\alpha|=|\beta|+1$. In particular terms of the form
\begin{equation} C\int_{r \leq |x|\leq R} x_l \partial_k\partial_jk_i(x) dx = 0 \end{equation} for C constant can be added freely.
The first two terms can be dealt with (roughly, this is jut a sketch) in two steps. The first is
\begin{eqnarray} && \int_{ |h|>|x+h-y|\geq \epsilon} \partial_k \partial_j k_i(x+h-y)\partial_lY(x+h)(x_l+h_l-y_l)f(y)dy \\ && = \int_{ |h|>|x+h-y|\geq \epsilon} \partial_k \partial_j k_i(x+h-y)\partial_lY(x+h)(x_l+h_l-y_l)\left(f(y)-f(x) \right) dy \leq c|f|_{\gamma}|h|^{\gamma} \end{eqnarray}
where we made use of the $\gamma$-Holder continuity of $f$ and the mean-value property of the kernel by adding a zero term. The second is
\begin{eqnarray} && \int_{ |x+h-y|\geq|h|,|x-y|\geq 2|h|} \Big( \partial_k \partial_j k_i(x+h-y)(x_l+h_l-y_l) - \partial_k \partial_j k_i(x-y)(x_l-y_l) \Big)\left(f(y)-f(x+h) \right) dy \\ && \leq c|f|_{\gamma}|h|^{\gamma} \end{eqnarray}
through the use of the mean-value theorem on the integral kernel $\partial_k \partial_j k_i$ and the further addition of zero terms.
--------------- Remaining Problems -----------------
I have trouble bounding the difference of the remainder terms by $C|h|^{\gamma}$. Following a similar method as above I can bound them by $C|h|^{\xi}$ for $\xi < \gamma$. But perhaps there is a more general theorem I could use? I'd also be interested in textbook recommendations covering similar subject matter.