Testing the time dependent Schrodinger Equation with an analytical solution?

Question

I am numerically solving the Schrodinger Equation in 1D first and in higher dimension later, but I want to know the convergence rate of my numerical solver in different grid size and numerical methods.

To calculate the error (or convergence rate), I will need an analytical solution to compare with, but I don't seem to find one for TDSE, even in 1D when potential is 0.

For example, numerically solved Poisson equation can be tested with an analytical solution sin(x). Is there a simple analytical solution for 1D TDSE, or in higher dimension?

Thank you

Answer 1

It depends what you consider as an analytical solution. If you consider a Fourier series as an analytical solution, you have this:

$$i \hbar \frac{\partial \Psi (\mathbf{r},t)}{\partial t} = -\frac{\hbar^{2}}{2m} \nabla^{2} \Psi(\mathbf{r},t)$$

Let's say you can use separation of variables and assume this for your wave function as:

$$\Psi(\mathbf{r},t) = F(\mathbf{r}) T(t)$$

So:

$$i \hbar F T^{'} = -\frac{\hbar^{2}}{2m} T \nabla^{2} F$$

or:

$$i \frac{2 m}{\hbar} \frac{T^{'}}{T} = -\frac{1}{F} \nabla^{2} F = \lambda^{2}$$

You have this for temporal part:

$$T(t) = \exp(-i \frac{\hbar}{2m} \lambda^{2} t)$$

For spatial part:

$$\nabla^{2}F + \lambda^{2} F = 0$$

The spatial part equation is just a Helmholtz equation. Now, it depends what geometry you have here. If it is just a cube in three-dimensional space:

$$F(\mathbf{r}) = X(x)Y(y)Z(z)$$

You have:

$$X^{''}YZ + XY^{''}Z + XYZ^{''} + \lambda^{2} XYZ = 0$$

or:

$$\frac{X^{''}}{X} + \frac{Y^{''}}{Y} + \frac{Z^{''}}{Z} + \lambda^{2} = 0$$

or:

$$\frac{X^{''}}{X} = - \frac{Y^{''}}{Y} - \frac{Z^{''}}{Z} - \lambda^{2} = -k_{x}^{2}$$

So:

$$X^{''} + k_{x}^{2} X = 0$$

$$Y^{''} + k_{y}^{2} Y = 0$$

$$Z^{''} + k_{z}^{2} Z = 0$$

and:

$$k_{x}^{2} + k_{y}^{2} + k_{z}^{2} = \lambda^{2}$$

$$X(x) = A_{x} \cos{(k_{x}x)} + B_{x} \sin{(k_{x}x)}$$

$$Y(y) = A_{y} \cos{(k_{y}y)} + B_{y} \sin{(k_{y}y)}$$

$$Z(z) = A_{z} \cos{(k_{z}z)} + B_{z} \sin{(k_{z}z)}$$

If you have Dirichlet boundary condition of zero at all boundaries of cube, you have:

$$X(x) = B_{x} \sin{(\frac{n\pi}{L_{x}}x)}$$

$$Y(y) = B_{y} \sin{(\frac{m\pi}{L_{y}}y)}$$

$$Z(z) = B_{z} \sin{(\frac{l\pi}{L_{z}}z)}$$

Where $L_{x}$, $L_{y}$, and $L_{z}$ are the length of the cube, $n$, $m$, and $l$ are integers, so:

$$k_{x} = \frac{n \pi}{L_{x}}$$

$$k_{y} = \frac{m \pi}{L_{y}}$$

$$k_{z} = \frac{l \pi}{L_{z}}$$

$$\lambda^{2}_{n,m,l} = \Bigg(\frac{n \pi}{L_{x}}\Bigg)^{2} + \Bigg(\frac{m \pi}{L_{y}}\Bigg)^{2} + \Bigg(\frac{l \pi}{L_{z}}\Bigg)^{2}$$

Finally:

$$\Psi_{n,m,l}(x,y,z,t) = B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i \frac{\hbar}{2m} \lambda^{2}_{n,m,l} t)$$

and

$$\Psi(x,y,z,t) = \sum_{n}\sum_{m}\sum_{l} B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i \frac{\hbar}{2m} \lambda^{2}_{n,m,l} t)$$

If for initial condition: $$\Psi(x,y,z,0) = \Psi_{0}(x,y,z)$$

$$B_{n,m,l} = \frac{1}{L_{x}L_{y}L_{z}}\int_{0}^{L_{x}} \int_{0}^{L_{y}} \int_{0}^{L_{z}} \Psi_{0}(x,y,z) \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} dx dy dz$$

Note that the frequency of your system is defined as:

$$\omega = \frac{\hbar}{2m} \mathbf{k} \cdot \mathbf{k}$$

Where:

$$\mathbf{k} = (k_{x},k_{y},k_{z})$$

The energy of the system is defined as:

$$E = \hbar \omega = \frac{\hbar^{2}}{2m} \mathbf{k} \cdot \mathbf{k}$$

So rewriting the solution:

$$\Psi(x,y,z,t) = \sum_{n}\sum_{m}\sum_{l} B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i\omega_{n,m,l} t)$$

Where:

$$\omega_{n,m,l} = \frac{\hbar}{2m} \Bigg( \Big ( \frac{n \pi}{L_{x}} \Big )^{2} + \Big ( \frac{m \pi}{L_{y}} \Big )^{2} + \Big ( \frac{l \pi}{L_{z}} \Big )^{2} \Bigg)$$

So, now you have an analytical solution by means of Fourier series to examine the error and convergence rate your numerical implementation.

Update: For an initial Gaussian wave packet of the form:

$$\Psi(x,y,z,0) = \Psi_{0}(x,y,z) = \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{\sigma_{0}^{2} \pi}} \exp{(i \mathbf{k} \cdot \mathbf{r} - \frac{|\mathbf{r} - \mathbf{r}_{0}|^{2}}{2 \sigma_{0}^{2}})}$$

Where: $\mathbf{k} = (k_{x}, k_{y}, k_{z})$, $\mathbf{r} = (x, y, z)$, $\mathbf{r}_{0} = (\frac{L_{x}}{2}, \frac{L_{y}}{2}, \frac{L_{z}}{2})$.

So:

$$B_{n,m,l} = \frac{1}{L_{x}L_{y}L_{z}}\int_{0}^{L_{x}} \int_{0}^{L_{y}} \int_{0}^{L_{z}} \Psi_{0}(x,y,z) \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} dx dy dz$$

But this 3D integral could be broken down to three 1D integrals of:

$$\mathcal{I}_{x} = \frac{1}{L_{x}} \int_{0}^{L_{x}} \exp{(i k_{x}x - \frac{(x - \frac{L_{x}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{x} x)} dx$$

$$\mathcal{I}_{y} = \frac{1}{L_{y}} \int_{0}^{L_{y}} \exp{(i k_{y}y - \frac{(y - \frac{L_{y}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{y} y)} dy$$

$$\mathcal{I}_{z} = \frac{1}{L_{z}} \int_{0}^{L_{z}} \exp{(i k_{z}z - \frac{(z - \frac{L_{z}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{z} z)} dz$$

But:

$$\mathcal{I}_{x} = i\frac{\sigma_{0}}{2L_{x}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{x}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{x}^{2} + i k_{x} L_{x})} (\mathrm{erf} (\frac{- \frac{L_{x}}{2} + i 2 \sigma_{0}^{2} k_{x}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{x}}{2} + i 2 \sigma_{0}^{2} k_{x}}{\sqrt{2} \sigma_{0}})))$$

$$\mathcal{I}_{y} = i\frac{\sigma_{0}}{2L_{y}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{y}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{y}^{2} + i k_{y} L_{y})} (\mathrm{erf} (\frac{- \frac{L_{y}}{2} + i 2 \sigma_{0}^{2} k_{y}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{y}}{2} + i 2 \sigma_{0}^{2} k_{y}}{\sqrt{2} \sigma_{0}})))$$

$$\mathcal{I}_{z} = i\frac{\sigma_{0}}{2L_{z}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{z}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{z}^{2} + i k_{z} L_{z})} (\mathrm{erf} (\frac{- \frac{L_{z}}{2} + i 2 \sigma_{0}^{2} k_{z}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{z}}{2} + i 2 \sigma_{0}^{2} k_{z}}{\sqrt{2} \sigma_{0}})))$$

So:

$$B(k_{x},k_{y},k_{z}) = \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{\sigma_{0}^{2} \pi}} \mathcal{I}_{x} \mathcal{I}_{y} \mathcal{I}_{z}$$

Answer 2

If you have a solution $\psi$ to the stationary Schroedinger equation $$ H\psi(x) = E \psi(x) $$ then the time dependent Schroedinger equation $$ i\hbar \frac{\partial}{\partial t}\Psi(x,t) = H\Psi(x,t) $$ has the solution $\Psi(x,t) = e^{-iE t/\hbar}\psi(x)$.

(In other words, every pure eigenstate remains a pure eigenstate -- only the phase rotates with a frequency proportional to $E$.)

Answer 3

Take the one-dimensional harmonic oscillator and its (analytical) ground state function $\Psi_0(x)$ (which is just a Gaussian) and displace it by a constant length $a$. Use the resulting state $\Psi(x) = \Psi_0(x-a)$ as initial state to the time-dependent Schrödinger equation. One knows that it should oscillate with a frequency of $1/2\pi$ (in atomic units). Evolve the wavepacket for a chosen number of periods, and compare the final state with the initial state (which for an exact propagation method should be identical).

Here is a picture I produced more than a decade ago ($M$ in this case is the order of the generalized Crank-Nicolson method). I hope this gets the idea.

Answer 4

AloneProgrammer and Wolfgang provided a general answer, which are very useful for a general initial condition. Here I want to provide a more explicit answer, which is simply the analytical solution of a 1D box.

Let's assume $V=0$, and the TDSE has the following form $$i\dfrac{\partial \psi}{\partial t} = -\dfrac{1}{2}\dfrac{\partial^2\psi}{\partial x^2}$$

If the boundary condition is $\psi(x=0) = \psi(x=1) = 0$. The time evolution of the real part follows

$$Re(\psi(x,t)) = \sqrt{2}\cos(\frac{n^2 \pi^2}{2}t)\sin(n\pi x)$$

where n is the number of node that you want. You can now compare it with the numerical solution of the TDSE!