  We study the large-time behavior in all $L^p$ norms and in different
space-time scales of solutions to a heat equation with a Caputo $\alpha$-time
derivative posed in $\mathbb{R}^N$. The initial data are assumed to be
integrable, and, when required, to be also in $L^p$. A main difficulty in the
analysis comes from the singularity in space at the origin of the fundamental
solution of the equation when~$N>1$.
  The rate of decay in $L^p$ norm in the characteristic scale, $|x|\asymp
t^{\alpha/2}$, dictated by the scaling invariance of the equation, is
$t^{-\frac{\alpha N}{2}(1-\frac1p)}$. In compact sets it is $t^{-\alpha/2}$ for
$N=1$, $t^{-\alpha}$ for $N\ge 3$, and $t^{-\alpha}\log t$ in the critical
dimension $N=2$. In intermediate scales, going to infinity but more slowly than
$t^{\alpha/2}$, we have an intermediate decay rate. In fast scales, going to
infinity faster than $t^{\alpha/2}$, there is no universal rate, valid for all
solutions, as we will show by means of some examples. Anyway, in such scales
solutions decay faster than in the characteristic one.
  When divided by the decay rate, solutions behave for large times in the
characteristic scale like $M$ times the fundamental solution, where $M$ is the
integral of the initial datum. The situation is very different in compact sets,
where they converge to the Newtonian potential of the initial datum if $N\ge
3$, one of the main novelties of the paper, and to a constant if $N=1,2$. In
intermediate scales they approach a multiple of the fundamental solution of the
Laplacian if $N\ge 3$, and a constant in low dimensions. The asymptotic
behavior in scales that go to infinity faster than the characteristic one
depends strongly on the behavior of the initial datum at infinity. We give
results for certain initial data with specific decays.
