As we all learned in a first course in PDEs, for "reasonable" initial data $u_0$, the solution to the linear heat equation \begin{equation} \partial_t u = \triangle u \end{equation} exists classically and converge to zero uniformly.
The proof can be given using the explicit Gaussian kernel \begin{equation} K_t(x) = \frac{1}{(4\pi t)^{d/2}} \exp( -|x|^2 / 4t ) \end{equation} for which \[ u(t,x) = K_t \star u_0 \] solves the heat equation. The fact that $K_t$ is a Schwartz function for every $t > 0$, and that it depends smoothly on $t$, means that as long as $u_0$ is locally integrable and has no more than polynomial growth (for example), then $u(t,x)$ is well-defined everywhere and smooth.
Furthermore, by H\"older's inequality we see \begin{equation} |u(t,x)| \leq |K_t|_{L^p} |u_0|_{L^q} \end{equation} where $p$ and $q$ are conjugate exponents. Observe that by explicit computation, $K_t$ always have $L^1$ norm $1$, but has decreasing $L^p$ norm for every $p > 1$. This shows that, in particular, if $u_0 \in L^q$ for any $q\in [1,\infty)$, then the classical solution converges uniformly to zero.
Notice that the theory leaves out $L^\infty$. In particular, we see easily that if $u_0 \equiv C$, then the convolution defines the classical solution $u(t,x) \equiv C$ as the solution to the heat equation, and this does not converge to $0$.
An obvious question is: is it the case that for bounded, locally integrable functions $u_0$, that convergence to a constant is always the case for solutions to the linear heat equation?
On compact manifolds
Our case is bolstered by the situation on compact manifolds. When solving the heat equation on compact manifolds, if the initial data is bounded, it is also in $L^2$. So we can decompose the initial data using the eigen-expansion relative to the Laplace-Beltrami operator $\triangle$. The system then diagonalizes into decoupled ordinary differential equations which, except for the kernel of $\triangle$,have solutions that are exponentially decaying.
In fact, this is sufficient to tell us that on compact manifolds, solutions to the heat equation converge exponentially to its mean.
We immediately however see a difference with the case on $\mathbb{R}^d$. On Euclidean space, the solution frequently only decay polynomially to zero, when the data is in $L^p$. One can check this by letting $u_0$ be a Gaussian function, for which the heat equation solution decay at rate $t^{-d/2}$.
Solutions are asymptotically locally constant
On the positive side, we do have that, for every $u_0$ that is bounded, the solutions are asymptotically locally constant at every time. This is because if $u_0$ is bounded (say by the number $b$), we have that \[ u(t,x) - u(t,y) = \int (K_t(x-z) - K_t(y-z))u_0(z) ~\mathrm{d}z \] so \begin{equation} |u(t,x) - u(t,y)| \leq b \int |K_t(x-z) - K_t(y-z) | ~\mathrm{d}z \end{equation} Performing an explicit change of variables we get \begin{equation} |u(t,x) - u(t,y)| \leq \frac{b}{\pi^{d/2}} \int | \exp(- |z|^2) - \exp(- |z - (y-x)/ (2\sqrt{t})|^2) | ~\mathrm{d}z. \end{equation} For $(y-x)/\sqrt{t}$ sufficiently small, the integral can be very roughly bounded by $O(|y-x| / \sqrt{t})$. This shows that on every compact set $K$, \begin{equation} \sup_K u(t,\cdot) - \inf_K u(t,\cdot) \lesssim t^{-1/2} \end{equation} giving that the solutions are asymptotically locally constant.
This argument, however, says nothing about whether the "constant" value is the same between different times.
Oscillatory solutions
In fact, it is fairly easy to construct a bounded initial data for which the solution, evaluated at $x = 0$, fails to converge as $t \to \infty$. Performing the explicit change of variable we get \begin{equation} u(t,0) = \int \frac{1}{(4\pi)^{d/2}} \exp( - |z|^2/4) u_0( - \sqrt{t} z) \mathrm{d}z. \end{equation} Denote by $\mu_G$ the Gaussian measure $(4\pi)^{-d/2} \exp( - |z|^2/4) \mathrm{d}z$ (which we observe has total mass 1) we see that we can write $u(t,0)$ as the integral of a rescaling of $u_0$ against $\mu_G$.
First choose $\epsilon \ll 1$. We can choose $0 < \lambda < \Lambda < \infty$ such that \[ \int_{ |z| \leq \lambda} \mu_G < \epsilon, \quad \int_{|z| \geq \Lambda} \mu_G < \epsilon. \] Let $\kappa = \Lambda / \lambda$. Let $u_0$ be defined by \begin{equation} u_0(x) = \begin{cases} 0 & |x| = \frac12 \newline 1 & |x| \in [\lambda \kappa^n, \lambda \kappa^{n+1}), n \text{ is even}\newline 0 & |x| \in [\lambda \kappa^m, \lambda \kappa^{m+1}), m \text{ is odd} \end{cases} \end{equation} Observe that this function satisfies \[ u_0(\kappa^{2n} x) = u_0(x) \] and \[ u_0(x) + u_0(\kappa x) \equiv 1.\] By construction, there exists a value $\gamma \in [1 - 2\epsilon, 1]$ such that \[ \int u_0(x) \mu_G = \gamma, \quad \int u_0(\kappa x) \mu_G = 1 - \gamma.\] Returning to the formula for $u(t,0)$, we see that \begin{equation} u(\kappa^{2n},0) = \begin{cases} \gamma & n \text{ is even} \newline 1 - \gamma & n \text{ is odd} \end{cases} \end{equation} This shows an oscillatory behavior at $x = 0$, showing that the limit as $t\to \infty$ of $u(t,0)$ to not exist.
On the other hand, this particular solution has the nice property that, on any compact set $K$, and for any fixed $t_0 > 0$, the sequence of functions $v_n(x) = u(t_0 \kappa^{4n}, x)$ converges uniformly to $u(t_0,0)$ as $n \to \infty$.
