Fortunately your boundary conditions are very simple, such that you can express the analytical solution of your problem explicitly. In fact, by using the ansatz of separation of variables one derives that for dimension $d=1$ and sufficiently smooth initial function $u_0$ the unique solution of your problem can be derived by Fourier expansion and it is given by
\begin{align}\tag{1}\label{1}
u(x,t)=\sum_{n=1}^\infty B_n\sin(n\pi x)e^{-n^2\pi^2t},
\end{align}
where
\begin{align*}
B_n:=2\int_0^1\sin(n\pi x)u_0(x)\,dx
\end{align*}
for $n\in\mathbb{N}.$
In order to verify $\partial_tu(x,t)=\partial_x^2u(x,t)$ for all $t>0$ and $x\in[0,1],$ one has to justify the identities
\begin{align}\tag{2}\label{2}
\partial_tu(x,t)=\sum_{n=1}^\infty\partial_t B_n\sin(n\pi x)e^{-n^2\pi^2t}
\end{align}
and
\begin{align}\tag{3}\label{3}
\partial_x^2u(x,t)=\sum_{n=1}^\infty\partial_x^2 B_n\sin(n\pi x)e^{-n^2\pi^2t}
\end{align}
for all $t>0$ and $x\in[0,1].$ In oder to prove e.g. identity \eqref{2}, it suffices to show that \eqref{1} converges uniformly, the series of derivatives $$\sum_{n=1}^\infty\partial_t B_n\sin(n\pi x)e^{-n^2\pi^2t}$$
converges uniformly and the derivatives $\partial_t B_n\sin(n\pi x)e^{-n^2\pi^2t}$ are continuous for all $n\in\mathbb{N}.$
Those uniform convergence results can easily be verified by using the Weierstraß M-test. Now, identity \eqref{3} can be shown in the exact same way.
Furthermore, it is straightforward to see that $u(0,t)=u(1,t)=0$ for all $t>0.$ Finally, it holds
\begin{align*}
u(x,0)&=\sum_{n=1}^\infty B_n\sin(n\pi x)\\
&=2\sum_{n=1}^\infty\sin(n\pi x)\int_0^1\sin(n\pi x)u_0(x)\,dx\\
&=u_0(x)
\end{align*}
for sufficiently smooth initial conditions $u_0(x)$ satisfying $u_0(0)=u_0(1)=0.$ Minimal regularity assumptions for the justification of this Fourier Sine expansion can be found in the literature.
Exactly as above one now derives
\begin{align*}
\partial_xu(x,t)&=\sum_{n=1}^\infty\partial_x B_n\sin(n\pi x)e^{-n^2\pi^2t}\\
&=\sum_{n=1}^\infty n\pi B_n\cos(n\pi x)e^{-n^2\pi^2t},
\end{align*}
such that elementary estimates show
\begin{align*}
|\partial_x u(x,t)|\leq 2\pi\int_0^1|u_0(\tilde x)|\,d\tilde x\sum_{n=1}^\infty ne^{-n^2\pi^2t},
\end{align*}
which vanishes uniformly in $x\in[0,1]$, as $t\to\infty.$ Furthermore, if $\int_0^1|u_0(\tilde x)|\,d\tilde x$ is chosen sufficiently small, it also holds $|\partial_xu(x,t)|\leq \varepsilon$ for $t>0,\,x\in[0,1]$ and a given $\varepsilon>0.$
The same arguments should also hold true for dimensions $d>1.$ In fact, two dimensional separation of variables yields the candidate
\begin{align*}
u(x,y,t)=\sum_{m,n=1}^\infty B_{m,n}\sin(m\pi x)\sin(n\pi y)e^{-(m^2+n^2)\pi^2t},
\end{align*}
where
\begin{align*}
B_{m,n}:=4\int_0^1\int_0^1\sin(m\pi x)\sin(n\pi y)u_0(x,y)\,dxdy
\end{align*}
for $m,n\in\mathbb{N}.$
By repeating the above arguments you should therefore also be able to obtain your desired decay and smallness properties for $d=2.$ As far as i can judge, the extension to $d\in\mathbb{N}$ is now more or less straightforward.
