Let $\mathcal{H}$ be a Hilbert space and let $(A(t))_{0\leqslant t\leqslant 1}$ be a family of operators in $\mathcal{H}$. Consider the differential equation
$$
x'(t)=A(t)x(t). \tag{1}
$$
If, say $A:[0,1]\rightarrow B(\mathcal{H})$ is a continuous map into the bounded operators on $\mathcal{H}$, then there is a unique solution $x:[0,1]\rightarrow B(\mathcal{H})$ to equation (1), and that unique solution is given by the Peano-Baker series
$$
x(t) = \sum_{n=0}^\infty \int_{t>s_1>\ldots >s_n>0} A(s_1)\cdots A(s_n) ds_n\cdots ds_1.
$$
Furthermore, if $A(t)=A$ for some fixed bounded operator, then the Peano-Baker series reduces to the exponential function series
$$
x(t) = \sum_{n=0}^\infty \frac{t^nA^n}{n!}=e^{tA} \tag{2}.
$$
Now, consider the heat equation, which is equation (1) with $A(t)=\Delta$, where $\Delta$ denotes the Laplacian. Here, $\Delta$ is no longer bounded, but rather an unbounded, self-adjoint, non-positive operator in $\mathcal H=L^2(\mathbb{R}^d)$. Then it is still correct that $x(t)=e^{t\Delta}$ solves equation (1), but the exponential is now understood in the sense of functional calculus of a self-adjoint operator. Since $\Delta\leqslant 0$, one finds $\lVert x(t)\rVert \leqslant 1$, i.e. $x(t)$ is still bounded for all $t$. It is then natural to ask if we also have the series expansion (2), corresponding to the Peano-Baker series.
However, it is clear that we cannot have convergence of the series with respect to the operator norm, since $\Delta$ is unbounded. Nonetheless, it is reasonable to expect
$$
x(t)\psi=\sum_{n=0}^\infty \frac{t^n\Delta^n\psi}{n!} \tag{3}
$$
for some suitable subspace of functions $\psi\in L^2(\mathbb{R}^d)$. For instance, any function $\psi$ which satisfies $\hat\psi(k)=0$ if $\lvert k\rvert > M$ (where $\hat\psi$ denotes the Fourier transform of $\psi$) will also satisfy
$$
\lVert \Delta^n \psi\rVert^2=\int \lvert k\rvert ^{4n}\lvert \hat\psi(k)\rvert^2 \, dk \leqslant M^{4n} \int \lvert \hat\psi(k)\rvert^2 \, dk = M^{4n} \lVert \psi\rVert^2,
$$
and therefore
$$
\sum_{n=0}^\infty \frac{t^n \lVert \Delta^n \psi\rVert}{n!} \leqslant \sum_{n=0}^\infty \frac{t^n M^{2n}\lVert \psi\rVert}{n!} = e^{t M^2}\lVert \psi\rVert < \infty.
$$
It follows that the series $\sum_{n=0}^\infty \frac{t^n\Delta^n\psi}{n!}$ converges! With some additional work, one can also show that the identity (3) in fact holds true for such $\psi$.
On the other hand, consider the Gaussian function $\psi(x)=(2\pi)^{-1/4}e^{-\lvert x \rvert^2/4}$. Recalling that the moments of the normal distribution are well known, we find
$$
\lVert \Delta^n \psi\rVert^2 = \int \lvert k\rvert^{4n}\lvert\psi(k)\rvert^2\, dk = (4n-1)!!= \frac{(4n)!}{2^{2n}(2n)!}.
$$
In particular, by Stirling's approximation,
$$
\frac{t^n \lVert \Delta^n \psi\rVert }{n!} = \frac{t^n}{2^n}\sqrt{\frac{(4n)!}{(2n)!(n!)^2}} \approx \frac{2^{1/4}}{\sqrt{2\pi}} \frac{t^n 4^{n}}{\sqrt{n}}.
$$
Thus, when $t>1/4$, the series (3) does not converge with $\psi(x)=(2\pi)^{-1/4}e^{-\lvert x \rvert^2/4}$.
The conclusion is that one should be careful when considering the Peano-Baker series as the solution to equation (1). Even the Gaussian function, which is well-behaved by most standards, turns out to be problematic if one tries to plug it directly into the series expansion in the case of the heat equation.
One resource for material on these topics is the book 'Methods of Modern Mathematical Physics Volume 2' by Reed and Simon.