Your textbook says the fundamental solution is $\Phi(x,t) = \frac{1}{(4\pi t)^{n/2}}e^{-\frac{|x|^2}{2t}}$ then we plug into the definition of the energy ball:
\begin{eqnarray} E(x,t,r) &=& \left\{ (y,s): s \leq t \text{ and } \frac{1}{(4\pi (t-s))^{n/2}}e^{-\frac{1}{2}\frac{|x-y|^2}{t-s}} \geq \frac{1}{r^2}\right\} \\
&=& \bigcup_{s \leq t} \left\{ (y,s): |x-y|^2 \leq -2(t-s)\log \frac{(4\pi (t-s))^{n/2}}{r^2}\right\} \end{eqnarray}
So the energy ball lives in (plain-old Galilean) space-time $\mathbb{R}^n \times \mathbb{R}$ and it's fibered by Euclidean balls. This is something like a Euclidean light-cone in physics.
Stochastic Point of View
One can show the heat equation can be solved by random walk. I believe the formula is $f(x) = \mathbb{E}[f(B_\tau)]$ where $\tau$ is the hitting time of a Brownian motion with $B_0 = x$ hitting the boundary $B_\tau \in \partial V$.
One can imagine heat diffusing by way of a Brownian motion. See Greg Lawler Heat Equation and Random Walk.
Actually come to think of it, you solve the heat equation by convolving the initial solution $u(t=0, x)$ with the heat kernel $\frac{1}{\sqrt{t}}e^{-x^2/t}$. Convolving with the heat kernel is as if diffusing the original solution via Brownian motion.
- for $t \ll 1$ this is a point distribution $\frac{1}{\sqrt{2\pi / t}}e^{-x^2/t} \to \delta(t)$
- for $t \gg 1$ this is uniform across space $\frac{1}{\sqrt{2\pi/ t}}e^{-x^2/t} \to \frac{1}{\sqrt{t}}$.
For linear PDE, that convolution can be thought of as just the Minkowski sum, or the theory of "wavefronts"as mentioned in the book of Hormander: Analysis of Linear Partial Differential Operators I or Tao's blog Computing Convolutions of Measures.
The mean value property is rather intuitive in the stochastic point of view:
$$\Delta^2 f \approx \frac{f(x+h,y)+ f(x,y+h)+f(x-h,y)+f(x,y-h) }{4} = \mathbf{E}f\big((x,y) + (\Delta x, \Delta y)\big)$$
where $(\Delta x, \Delta y)= (\pm 1, \pm 1)$ each with probability $\mathbb{P}=\frac{1}{4}$. Or in a very geometric way the Laplacian is just the average of the values of $f$ over a circle:
$$ \Delta^2f = \frac{1}{2\pi} \oint f\bigg((x,y) + \epsilon(\cos \theta, \sin \theta)\bigg)d\theta = f(x)$$