Solution (source code)

= Solution

Take the spatial <Fourier transform> with the same convention. The <partial differential equation> becomes $i\partial_t\widehat f-k^2\widehat f=0$, so $\widehat f(k,t)=e^{-ik^2t}\widehat f_0(k)$. Inverting and substituting the transform of the initial <function> yields the <convolution>
$$
f(x,t)=\int_{\mathbb R}f_0(\xi)K_t(x-\xi)d\xi,\qquad K_t(s)=\frac1{2\pi}\int_{\mathbb R}e^{iks-itk^2}dk.
$$
Complete the square, $iks-itk^2=is^2/(4t)-it[k-s/(2t)]^2$, and put $\ell=\sqrt t[k-s/(2t)]$. This gives
$$
\boxed{f(x,t)=\frac{c}{2\pi\sqrt t}\int_{-\infty}^{\infty}f_0(\xi)e^{i(x-\xi)^2/(4t)}d\xi,\qquad c=\int_{\mathbb R}e^{-i\ell^2}d\ell=\sqrt\pi e^{-i\pi/4}.}
$$
The oscillatory Gaussian is a <Fresnel integral>, understood as the limit with a positive Gaussian damping factor. That prescription also justifies the transform manipulations, first for rapidly decreasing smooth data and then by <L2 space> extension. The resulting kernel is the <free-particle propagator> in units $\hbar=1$, $m=1/2$; its phase and normalization agree with the multiplier $e^{-ik^2t}$.