Positivity of a convolution integral Let $f\in\mathcal{S}(\mathbb{R}^3)$ a real function. Consider the following integral
$$I_f:=\int_{\mathbb{R}^3}\int_{\mathbb{R}^3}\frac{f(x)f(y)}{|x-y|^2}dydx$$
Observe that,  either if $f$ is positive or negative, the integral becomes positive.
Is it true that $I_f\geq 0$ for a generic $f$?
 A: In the following, we will construct a suitable real-valued Schwartz
function $f\in\mathcal{S}\left(\mathbb{R}^{3}\right)$ satisfying
$$
A_{f}:=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{f\left(x\right)\cdot f\left(y\right)}{\left|x-y\right|^{2}}dy\,dx<0.
$$
For the idea, note that (as we will see below) $A_{f}$ is given by
a certain convolution, i.e., $A_{f}=\left\langle \phi\ast f,\,f\right\rangle _{\mathcal{S}',\mathcal{S}}$.
Using the Fourier transform, we see that $A_{f}=\left\langle \widehat{f}\cdot\widehat{\phi},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}}$.
We will show below that $\widehat{\phi}$ is given by integration
against a locally integrable function; in fact by integration against
$\widehat{\phi}\left(\xi\right)=\frac{\pi}{\left|\xi\right|}$.
Now, note for $U:=B_{1}\left(\left(\begin{smallmatrix}10\\
0\\0
\end{smallmatrix}\right)\right)$ (the open ball with radius $1$ around $(10,0,0)$) that $U\cap\left(-U\right)=\emptyset$. Indeed, for $\xi\in U\cap\left(-U\right)$
we would have
$$
18=2\left[\left|\left(\begin{matrix}10\\
0\\0
\end{matrix}\right)\right|-1\right]\leq2\left|\xi\right|=\left|\xi-\left(-\xi\right)\right|\leq\left|\xi-\left(\begin{matrix}10\\
0\\0
\end{matrix}\right)\right|+\left|-\left(\begin{matrix}10\\
0\\0
\end{matrix}\right)-\left(-\xi\right)\right|<2,
$$
a contradiction. Now, choose an arbitrary nonnegative function $h\in C_{c}^{\infty}\left(U\right)$
and define
$$
g\left(\xi\right):=i\cdot h\left(\xi\right)-i\cdot h\left(-\xi\right)\qquad\forall x\in\mathbb{R}^{3}.
$$
Note $\overline{g\left(-\xi\right)}=-i\cdot h\left(-\xi\right)+i\cdot h\left(\xi\right)=g\left(\xi\right)$,
so that $f:=\mathcal{F}^{-1}g\in\mathcal{S}\left(\mathbb{R}^{3}\right)$
is real-valued.
But (assuming $\widehat{\phi}\left(\xi\right)=\frac{\pi}{\left|\xi\right|}$,
which we will show below),
\begin{align*}
A_{f} & =\left\langle \widehat{f}\cdot\widehat{\phi},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}}\\
 & =\int_{\mathbb{R}^{3}}\frac{\pi}{\left|\xi\right|}\cdot\left[\widehat{f}\left(\xi\right)\right]^{2}\,d\xi\\
 & =\int_{\mathbb{R}^{3}}\frac{\pi}{\left|\xi\right|}\cdot\left[i\cdot h\left(\xi\right)-i\cdot h\left(-\xi\right)\right]^{2}\,d\xi\\
\left({\scriptstyle \text{support of }h\text{ and }h\left(-\bullet\right)\text{ is disjoint}}\right) & =-2\cdot\int_{\mathbb{R}^{3}}\frac{\pi}{\left|\xi\right|}\cdot\left[h\left(\xi\right)\right]^{2}\,d\xi\\
\left({\scriptstyle \text{integrand }\geq0\text{ and does not vanish a.e.}}\right) & <0.
\end{align*}
All that remains to show is that indeed $A_{f}=\left\langle \phi\ast f,\,f\right\rangle _{\mathcal{S}',\mathcal{S}}$
with $\widehat{\phi}\left(\xi\right)=\frac{\pi}{\left|\xi\right|}$.
To this end, define 
$$
\phi:\mathbb{R}^{3}\to\left[0,\infty\right),x\mapsto\begin{cases}
\left|x\right|^{-2}, & \text{if }x\neq0,\\
0, & \text{if }x=0
\end{cases}
$$
and note $\phi\in L_{{\rm loc}}^{1}\left(\mathbb{R}^{3}\right)$;
in fact, we have (using polar coordinates):
\begin{align*}
\int_{\left|x\right|<1}\left|\phi\left(x\right)\right|\,dx & =\int_{0}^{1}r^{3-1}\int_{S^{2}}\left|\phi\left(r\xi\right)\right|\,d\mathcal{H}^{2}\left(\xi\right)\,dr\\
 & =\mathcal{H}^{2}\left(S^{2}\right)\cdot\int_{0}^{1}\,dr<\infty,
\end{align*}
where $\mathcal{H}^{2}$ denotes the $2$-dimensional Hausdorff measure
and $S^{2}\subset\mathbb{R}^{3}$ is the unit sphere. Likewise, a
similar calculation shows
\begin{align*}
\int_{\left|x\right|>1}\left|\phi\left(x\right)\right|^{2}\,dx & =\mathcal{H}^{2}\left(S^{2}\right)\cdot\int_{1}^{\infty}r^{3-1}\cdot r^{-4}\,dr\\
 & =\mathcal{H}^{2}\left(S^{2}\right)\cdot\int_{1}^{\infty}r^{-2}\,dr\\
 & =\mathcal{H}^{2}\left(S^{2}\right)\cdot\left(-r^{-1}\right)\bigg|_{1}^{\infty}=\mathcal{H}^{2}\left(S^{2}\right)<\infty.
\end{align*}
Hence, we have shown $\phi\in L^{1}\left(\mathbb{R}^{3}\right)+L^{2}\left(\mathbb{R}^{3}\right)\subset\mathcal{S}'\left(\mathbb{R}^{3}\right)$.
Now, with this notation, it is not hard to see
\begin{align*}
\left\langle \phi\ast\overline{f},\,f\right\rangle _{\mathcal{S}',\mathcal{S}} & =\int_{\mathbb{R}^{3}}\left(\phi\ast\overline{f}\right)\left(x\right)\cdot f\left(x\right)\,dx\\
 & =\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\overline{f\left(x\right)}}{\left|x-y\right|^{2}}\,dy\,f\left(x\right)\,dx=A_{f}.
\end{align*}
Here, we take the pairing $\left\langle \cdot,\cdot\right\rangle _{\mathcal{S}',\mathcal{S}}$
to be bilinear (not sesquilinear). Furthermore, I use the normalization
$\widehat{f}\left(\xi\right)=\int f\left(x\right)\cdot e^{-2\pi i\left\langle x,\xi\right\rangle }\,dx$
of the Fourier transform. Now, the left-hand side of the above equation
can be rewritten as
$$
A_{f}=\left\langle \phi\ast\overline{f},\,f\right\rangle _{\mathcal{S}',\mathcal{S}}=\left\langle \widehat{\phi\ast\overline{f}},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}}=\left\langle \widehat{\overline{f}}\cdot\widehat{\phi},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}}=\left\langle \widehat{\overline{f}}\cdot\widehat{\phi},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}}=\left\langle \widehat{f}\cdot\widehat{\phi},\,\widehat{f}\right\rangle _{\mathcal{S}',\mathcal{S}},
$$
where the last step is only valid if $f$ is \emph{real-valued}. This
identity shows that the Fourier transform 
$$
\psi:=\widehat{\phi}\in\mathcal{F}\left(L^{1}+L^{2}\right)\subset C_{0}\left(\mathbb{R}^{3}\right)+L^{2}\left(\mathbb{R}^{3}\right)
$$
is important. In fact, if we can show that $\psi\geq0$ almost everywhere,
we are in business.
Now, recall that taking the Laplacian in Fourier domain corresponds
(up to a constant) to multiplication with $\left|x\right|^{2}$ in
space domain, which is nice for us since this will cancel out $\phi$.
More precisely, we have 
$$
\mathcal{F}^{-1}\left(\Delta\psi\right)=\sum_{j=1}^{3}\mathcal{F}^{-1}\left(\partial_{j}^{2}\psi\right)=\sum_{j}\left(-2\pi i\cdot x_{j}\right)^{2}\cdot\mathcal{F}^{-1}\psi=-4\pi^{2}\cdot\left|x\right|^{2}\cdot\mathcal{F}^{-1}\psi=-4\pi^{2}.
$$
Taking the Fourier transform of this identity yields
$$
\Delta\left(\frac{-\psi}{4\pi^{2}}\right)=\delta_{0},
$$
with $\left\langle \delta_{0},g\right\rangle _{\mathcal{S}',\mathcal{S}}=g\left(0\right)$
for $g\in\mathcal{S}\left(\mathbb{R}^{3}\right)$. But another solution
to this equation is given by the fundamental solution (cf. https://en.wikipedia.org/wiki/Fundamental_solution#Laplace_equation
and also https://en.wikipedia.org/wiki/Green\%27s_function_for_the_three-variable_Laplace_equation)
of the Laplace equation,
$$
\Gamma:\mathbb{R}^{3}\to\mathbb{R},x\mapsto\begin{cases}
\frac{-1}{4\pi\cdot\left|x\right|}, & \text{if }x\neq0,\\
0, & \text{if }x=0.
\end{cases}
$$
Hence, $\Delta\left(\Gamma-\frac{\psi}{4\pi^{2}}\right)=0$. Taking
again the Fourier transform of this identity, we arrive at
$$
-4\pi^{2}\cdot\left|\xi\right|^{2}\cdot\left[\widehat{\Gamma}+\frac{\widehat{\psi}}{4\pi^{2}}\right]=0.
$$
Since $\left|\xi\right|^{2}$ only vanishes at $0$, it is not too
hard to see that this implies 
$$
\mathrm{supp}\left[\widehat{\Gamma}+\frac{\widehat{\psi}}{4\pi^{2}}\right]\subset\left\{ 0\right\} ,
$$
which has to be understood in the sense of tempered distributions.
But the only distributions with this support condition are (finite(!))
linear combinations of $\partial^{\alpha}\delta_{0}$ with $\alpha\in\mathbb{N}_{0}^{3}$.
Hence, $\mathcal{F}\left[\Gamma+\frac{\psi}{4\pi^{2}}\right]\in\left\langle \partial^{\alpha}\delta_{0}\,:\,\alpha\in\mathbb{N}_{0}^{3}\right\rangle $.
By again taking the inverse Fourier transform, we conclude
$$
\Gamma+\frac{\psi}{4\pi^{2}}\in\left\langle X^{\alpha}\,:\,\alpha\in\mathbb{N}_{0}^{3}\right\rangle \qquad\text{ i.e., }\qquad\Gamma+\frac{\psi}{4\pi^{2}}=p
$$
for some polynomial $p$. We want to show $p\equiv0$. To this end,
recall $\psi\in C_{0}\left(\mathbb{R}^{3}\right)+L^{2}\left(\mathbb{R}^{3}\right)$.
Furthermore, $\Gamma\left(x\right)\to0$ as $\left|x\right|\to\infty$.
But for $g\in L^{2}\left(\mathbb{R}^{3}\right)$, it is well-known
that $g\left(nx\right)\xrightarrow[n\to\infty]{}0$ for almost all
$x\in\mathbb{R}^{3}$. Hence, we get $p\left(nx\right)\xrightarrow[n\to\infty]{}0$
for almost all $x\in\mathbb{R}^{3}$. This implies that $p\equiv0$
(I am quite sure that this is true, but at the moment I am lacking
a convincing proof. I will post one if I find one). Hence, we have
shown (up to this last slight problem) that
$$
\psi=-4\pi^{2}\Gamma=\left(x\mapsto\frac{\pi}{\left|x\right|}\right).
$$
Remark: There is probably some easier way of computing the
Fourier transform of $\psi=\mathcal{F}\phi$.
