Is there a proof of the fact that the operator $\mathfrak{D}$ defined by $\mathfrak{D}[f](x,y):=\dfrac{f(x)-f(y)}{x-y}$ is continuous from $H^n(\mathbb{R})$ to $H^m(\mathbb{R}^2)$ or even $L^2(\mathbb{R}^2)$ assuming that $n$ is sufficiently big (i guess $n\geq m+1$ would suffice) ?
Some insights why there exists $C>0$ such that for all $f\in H^n(\mathbb{R})$, $\|\mathfrak{D}[f]\|_{L^2(\mathbb{R}^2)}\leq C \|f\|_{H^1(\mathbb{R})}$ might be very useful to guess how to show the result for derivatives of $\mathfrak{D}[f]$.
Also some formula for the derivatives of $\mathfrak{D}[f]$ hold, it basically has the form of a sum of Taylor sums divided by $(x-y)^{\alpha-j}$ where $\alpha=m_1+m_2+1$ where $m_1$ (resp. $m_2$) is the order of the derivative with respect to $x$ (resp. $y$).
We start by the fact that $\mathfrak{D}[f](x,y)$ verifies the following identity:
$$\mathfrak{D}[f](x,y)=\int_0^1dtf'\left(x+t(y-x)\right)$$
Hence assuming $f$ is smooth enough (in $C^\infty_c(\mathbb{R})$ for example):
$$\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)=\int_0^1dt(1-t)^\alpha t^\beta f^{(\alpha+\beta+1)}\left(x+t(y-x)\right)$$
Now we want to show the $L^2$-ness of $\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)$ close to the singularity namely on $\{|x-y|<1\}$ and uniformly far from it.
$\textbf{Close to singularity:}$
The computation uses Fubini, change of variables on $y$ and Jensen with convex function is $x\mapsto x^2$ goes like
$$\iint_{|x-y|<1}|\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)|^2dxdy\leq\int_0^1dt(1-t)^{2\alpha} t^{2\beta}\int_{-1}^1dy'\int_\mathbb{R}dxf^{(\alpha+\beta+1)}\left(x+ty'\right)^2$$
The last integral is precisely the $L^2$ norm squared of $f^{(\alpha+\beta+1)}$ and the two first integrals are finite and independent of $f$.
$\textbf{Far from the singularity:}$ It can be shown, using Leibniz formula that $$\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)=\sum_{j=0}^{\beta}C_{\alpha,\beta,j}\dfrac{\left(f^{(j)}(y)-\displaystyle\sum_{k=0}^{\alpha-j}\dfrac{f^{(k+j)}(x)}{k!}(y-x)^{k}\right)}{(y-x)^{\alpha+\beta+1-j}}.$$ where $C_{\alpha,\beta,j}:=\binom{\beta}{j}(\alpha+\beta-j)!(-1)^{m\beta-j}$.
Hence by Jensen: $$\iint_{|x-y|>1}|\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)|^2dxdy\leq C\sum_{j}\int_\mathbb{R}dxf^{(j)}(y)^2\int_{|x-y|>1}\dfrac{dy}{(x-y)^{2\alpha+2\beta+2-2j}}+...$$
We have then decomposed $\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f](x,y)|^2$ into a sum of integrals on $\mathbb{R}$ of some derivatives $f$ squared times a integral that is finite. Hence it is smaller than the $L^2$ norm of $f^{(max(\alpha,\beta)}$ (the order of the higher derivatives appearing in the sum.)
Hence there exists $C>0$ such that $\|\partial^{\alpha}_x\partial^{\beta}_y\mathfrak{D}[f]\|_{L^2(\mathbb{R}))}\leq C\|f^{(\alpha+\beta+1)}\|_{L^2(\mathbb{R})}$