Calderon-Zygmund type estimates for nonlocal PDE with H\"older continuous kernel

We study interior $L^p$-regularity theory, also known as Calderon-Zygmund theory, of the equation $$\int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{K(x,y)\ (u(x)-u(y))\, (\varphi(x)-\varphi(y))}{|x-y|^{n+2s}}\, dx\, dy = \langle f, \varphi \rangle \quad \varphi \in C_c^\infty(\mathbb{R}^n).$$ For $s \in (0,1)$, $t \in [s,2s]$, $p \in [2,\infty)$, $K$ an elliptic, symmetric, H\"older continuous kernel, if $f \in \left (H^{t,p'}_{00}(\Omega)\right )^\ast$, then the solution $u$ belongs to $H^{2s-t,p}_{loc}(\Omega)$ as long as $2s-t < 1$. The increase in differentiability is independent of the H\"older coefficient of $K$. For example, our result shows that if $f\in L^{p}_{loc}$ then $u\in H^{2s-\delta,p}_{loc}$ for any $\delta\in (0, s]$ as long as $2s-\delta < 1$. This is different than the classical analogue of divergence-form equations ${\rm div}(\bar{K} \nabla u) = f$ (i.e. $s=1$) where a $C^\gamma$-H\"older continuous coefficient $\bar{K}$ only allows for estimates of order $H^{1+\gamma}$. In fact, it is another appearance of the differential stability effect observed in many forms by many authors for this kind of nonlocal equations -- only that in our case we do not get a "small" differentiability improvement, but all the way up to $\min\{2s-t,1\}$. The proof argues by comparison with the (much simpler) equation $$\int_{\mathbb{R}^n} K(z,z) (-\Delta)^{\frac{t}{2}} u(z) \, (-\Delta)^{\frac{2s-t}{2}} \varphi(z)\, dz = \langle g,\varphi\rangle \quad \varphi \in C_c^\infty(\mathbb{R}^n).$$ and showing that as long as $K$ is H\"older continuous and $s,t, 2s-t \in (0,1)$ then the "commutator" $$\int_{\mathbb{R}^n} K(z,z) (-\Delta)^{\frac{t}{2}} u(z) \, (-\Delta)^{\frac{2s-t}{2}} \varphi(z)\, dz - c\int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{K(x,y)\ (u(x)-u(y))\, (\varphi(x)-\varphi(y))}{|x-y|^{n+2s}}\, dx\, dy$$ behaves like a lower order operator.