Analytical solution for a Gaussian wave packet / free particle Is there any closed formula for the time-dependent gaussian wave packet?
$$\Psi(x,0) = \frac{1}{ \sqrt[4]{ 2 \pi {\sigma_x}^2 }} \exp \left( i k_0 x - \frac{{\left( x - x_0 \right)}^2}{4 {\left(\sigma_x \right)}^2} \right)$$
$$\Psi(x,t) = ?$$
I found some formulas for the case when $x_0=0$


*

*https://en.wikipedia.org/wiki/Wave_packet#Dispersive

*https://www.colorado.edu/physics/phys2170/phys2170_fa06/downloads/Gaussian.pdf
There is also these:


*

*https://ocw.mit.edu/courses/nuclear-engineering/22-02-introduction-to-applied-nuclear-physics-spring-2012/lecture-notes/MIT22_02S12_lec_ch6.pdf
But I am struggling in the algebra for including the $x_0$ in them.
Maybe I should go through all those Fourier integrals, but it would be a nice shortcut if someone could give a final formula.
Edit
My goal is to use it for validating numerical results.
 A: Given a wavefunction, $\psi(x,0)$, the corresponding quantity, $\psi(x,t)$, is obtained from the integral:
$$
\boxed{\psi(x,t) = \int_{-\infty}^{\infty}dy\,Z[x,t;y,0]\psi(y,0)\,}
$$
where $Z[x,t;y,0]$ is a path integral, $\int_{x,0}^{y,t} \mathcal{D}x\,e^{\frac{i}{\hbar}I[x]}$, and $I[x]$ the single particle action of interest. For free non-relativistic particles, $I[x]=\int_0^td\tau\frac{m}{2}(\frac{dx}{d\tau})^2$, and the correctly normalised path integral reads:
$$
Z[x,t;y,0] = \sqrt{\frac{m}{2\pi i\hbar t}}\,\exp\Big(i\frac{m(x-y)^2}{2\hbar t}\Big).
$$
The initial wave packet is of the form:
$$
\psi(y,0)=Ne^{\frac{i}{\hbar}ky}e^{-\frac{1}{2}b(y-x_0)^2}, 
$$
where you can read off the $b,N$ of interest from your formula, e.g. $b=1/(2\sigma_x^2)$. 
The objective is to compute $\psi(x,t)$, which amounts to carrying out the $y$ integral in the above boxed formula. This is a Gaussian and so can easily be carried out. 
You'll primarily need the result:
$$
\int_{-\infty}^{\infty} da\,e^{iza^2} = \sqrt{\frac{\pi i}{z}},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0,
$$
which follows from Cauchy's theorem applied to $\oint_C da\,e^{iza^2}=0$, where the contour $C$ is that in the figure:

It is convenient to write $z=\rho e^{i\phi}$ (with $0\leq \phi<\pi$), and then given any allowed $\phi$ you choose $\theta$ such that $\phi+2\theta=\pi/2$. This enables you to take $R\rightarrow \infty$ while ensuring the angular contour contributions vanish, and you are left with $\int_{-\infty}^{\infty}da\,e^{iza^2}=\sqrt{i}e^{-i\phi/2}\int_{-\infty}^{\infty}da\,e^{-\rho a^2}$. Evaluating the remaining Gaussian integral and rewriting the result in terms of $z=\rho e^{i\theta}$ leads to the displayed result. 
The next step is to redefine your integration variable, $a\rightarrow a'=a-b$ (with $b\in \mathbb{R}$), 
$$
\int_{-\infty}^{\infty} da\,e^{iz(a+b)^2} = \sqrt{\frac{\pi i}{z}},
$$
which implies:
$$
\int_{-\infty}^{\infty} da\,e^{iza^2+2izba} = \sqrt{\frac{\pi i}{z}}\,e^{-izb^2},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0.
$$
We are almost there. The only thing remaining to do is to notice that both the left- and right-hand sides of this last expression are analytic in $b$, so we can extend $b$ into the complex plane, $b\rightarrow w=b+ib'$. Finally, we redefine $w\rightarrow zw$. Therefore, for generic complex numbers $z$ (with ${\rm Im}\,z\geq0$, $z\neq0$) and $w$,
$$
\boxed{\,\int_{-\infty}^{\infty} da\,e^{iza^2+2iwa} = \sqrt{\frac{\pi i}{z}}\,e^{-iw^2/z},\quad {\rm if }\quad z,w\in \mathbb{C},\quad {\rm and}\quad{\rm Im}\,z\geq0,\quad z\neq0\,\,\,}
$$
This is all you need in order to evaluate the above integral that gives you $\psi(x,t)$. You just read off $z,w$ and $a$ as appropriate for your case of interest and apply this result. (Of course the same formula applies for general $x_0$.) The result is:
$$
\boxed{\psi(x,t) = \frac{1}{\sqrt{1+i\frac{\hbar b }{m}t}}\exp\Bigg\{\!\!-\frac{i}{\hbar}\frac{k^2t}{2m}\frac{\big[1+i\frac{\hbar b}{k}(x-x_0)\big]^2}{1+i\frac{\hbar b}{m}\,t}\Bigg\}\psi(x,0)\,\,}
$$
Recall that $
\psi(x,0)=Ne^{\frac{i}{\hbar}kx}e^{-\frac{1}{2}b(x-x_0)^2}.
$
