I am looking for a proof for the following fact: for $U>0,\beta>0,$ there exists $C>0,\epsilon>0$ such that $$\forall u\in [0,U],\alpha\in\left]0,1\right],\int\_0^u\int\_{[-1,1]^2}\int\_{[-1,1]^2} \frac{1}{r}e^{-\alpha^2|x-y|^2/r} \, dx\,dy\,dr\leq Cu^{\epsilon}\alpha^{-2\beta},$$ where $|\cdot|$ is the euclidean norm.



$\int_0^u\int_{[-1,1]^2}\int_{[-1,1]^2}\frac{1}{r}e^{-\alpha^2|x-y|^2/r} \, dx\,dy\,dr\leq Cu^{\epsilon}\alpha^{-2\beta}$
$\newcommand{\al}{\alpha}\newcommand{\be}{\beta}\newcommand\ep\epsilon$The integral in question is
\begin{equation\*}
 I(u,\al):=\int\_0^u \frac{dr}r \int\_{[-1,1]^2} dy \int\_{[-1,1]^2} dx\,e^{-\al^2|x-y|^2/r}. 
\end{equation\*}
Using the substitution $r=\al^2 s$ (suggested by Giorgio Metafune), we get
\begin{equation\*}
 I(u,\al)=J(u/\al^2),
\end{equation\*}
where
\begin{equation\*}
 J(v):=\int\_0^v \frac{ds}s \int\_{[-1,1]^2} dy \int\_{[-1,1]^2} dx\,e^{-|x-y|^2/s}
 =\int\_0^v \frac{ds}s\, K(s)^2,
\end{equation\*}
where
\begin{equation\*}
 K(s):=\int\_{-1}^1 da\int\_{-1}^1 db\,e^{-(a-b)^2/s}
 \le\int\_{-1}^1 da\int\_{-\infty}^\infty db\,e^{-(a-b)^2/s} 
 =\int\_{-1}^1 da\, \frac c2\,\sqrt s=c\sqrt s;
\end{equation\*}
here and in what follows $c$ denotes universal positive real constants (possibly different even within the same formula), and $s$ is any positive real number. Also,
\begin{equation\*}
 K(s)\le\int\_{-1}^1 da\int\_{-1}^1 db=4.
\end{equation\*}
So,
\begin{equation\*}
 K(s)^2\le c\min(1,s).
\end{equation\*}
So,
\begin{equation\*}
 J(v)\le c\int\_0^v \frac{ds}s\,\min(1,s)=c(\min(1,v)+\ln\max(1,v))\le c\ln(1+ v); 
\end{equation\*}
the latter inequality follows because (say) $\min(1,v)+\ln\max(1,v)$ and $\ln(1+ v)$ are (i) positive and continuous in $v\in(0,\infty)$ and (ii) asymptotically equivalent to each other as $v\downarrow0$ and as $v\to\infty$.


So,
\begin{equation\*}
 I(u,\al)\le c\ln\Big(1+\frac u{\al^2}\Big). 
\end{equation\*}


So, upon the substitution $u=\al^2 v$, the inequality in question reduces to the inequality  

\begin{equation\*}
 \ln(1+v)\le C v^\ep\al^{2\ep-2\be} \tag{1}\label{1}
\end{equation\*}
for some real $C>0$, some real $\ep>0$, all real $v>0$, and all $\al\in(0,1]$.
If we now take $\ep=\min(1,\be)$, then \eqref{1} will hold for some real $C>0$ depending only on $\be$. $\quad\Box$


