Why this abuse of notation correctly solves the heat equation

Here's a simpler example. Let $v(t)$ be a vector-valued function of a real variable, and consider the differential equation

$$v'(t) = A v(t)$$

where $A$ is a square matrix. You can again "pretend that $A$ is a constant" to arrive at the formal solution

$$v(t) = e^{At} v(0)$$

and if you now interpret this exponential as a matrix exponential this is in fact the correct solution. The reason is exactly the same as in the constant case: it's because the derivative of $e^{At}$ is $A e^{At}$.

A similar but more complicated (with more analytic details) phenomenon occurs when $A$ is a linear operator on an infinite-dimensional vector space, with the right hypotheses.

Abstractly, time evolution is an exponential when the system does not depend on time. And that's what you have here. If you know the state $f$ of the system (function of x) at some time $t$, then the state $E(t)f$ observed $t$ seconds later is a function only of the system (which does not depend on $t$ in this discussion,) and depends on time $t$. Therefore, if you start with state $f$, then evolve through $t$ seconds, and then evolve this new state $E(t)f$ through $s$ seconds, you must obtain the same as if you were to start with $f$ and evolve through $t+s$ seconds. In other words, $$ E(s)(E(t)f) = E(t+s)f \mbox{ or } E(s)E(t)=E(s+t). $$ This applies to any Physical system whose characteristics do not vary with time, assuming, of course, that the rules of the Universe don't change along the way. And $E(0)=I$ because you're evolving through no time. This is a very abstract principle first described in operator form by Oliver Heaviside in the late 1800's. He was the really the first to describe things in terms of operators, and to play with them in an algebraic way. He discovered this exponential property, which eventually morphed into what we now call the Laplace transform. Heaviside's operator methods included the $D$ operator method in ODE's, with annihilator, partial fractions, and the cover-up method for partial fractions. These were devised by Heaviside to solve electric circuit equations. The exponential solution operator was a brilliant observation by Heaviside, but it definitely confused people at the time.

When you add stability in some norm sense, which is $\lim_{t\downarrow 0}E(t)f=f$, then there is a way to interpret this in terms of an exponential $e^{tA}$ where $A$ is the "generator" of the semigroup. The generator is computed through a derivative: $$ Af = \lim_{t\downarrow 0}\frac{1}{t}\{S(t)-I\}f. $$ In your case, you already know the generator $$ \lim_{t\downarrow 0}\frac{u(x,t)-u(x,0)}{t}=u_{t}(x,0)=\Delta u $$ So the generator is $A=\Delta$. You found a way to represent this exponential operator $E$ with the Fourier transform. Adding the stability condition at $t=0$ makes the semigroup unique. You found it, and it is unique, and it is stable at $t=0$ in every $L^p$ for $1 \le p < \infty$.