# Numerical solution of zero-potential time-dependent Schrödinger equation in 1D

I want to solve numerically the one-dimensional time-dependent Schrödinger equation $$i \psi_t(x,t)=-\frac{\hbar}{2m} \psi''(x,t)$$

My issue is that I don't have the physical background to understand what are the correct boundary conditions/initial state and I don't know how to know if the solution I get is the correct one. So I want to try to reproduce the same solution I found on wikipedia, shown below: What I've seen is that one usually discretize $$\psi''(x,t)$$ with the usual centeral finite difference scheme $$\psi''(x_i,t) = \frac{\psi_{i+1}(t) - 2\psi_i(t) + \psi_{i-1}(t)}{dx^2} + \mathcal{O}(h^2)$$

and hence the PDE becomes a system of ODEs that I can solve with an appropriate method.

Here are my questions:

• What boundary conditions do I have to impose to have a behaviour like in the figure? The solution does not appear to have a "fixed" value. How can I impose them (I'd need an answer in terms of what entries of the matrix I should change)
• What could be an initial condition do I have to impose to have a "wave" like the one in the picture?

Following the suggestion of @AloneProgrammer, I focus on the particle in a box case, where my domain now is $$[0,1]$$ and I have $$0$$ potential inside the domain, and $$V(x) = \infty$$ outside. In this configuration, the boundary conditions at $$0$$ and $$L$$ are Dirichlet homogeneous,i.e. $$\psi(0,t) = \psi(L,t) = 0$$

Hence the PDE becomes (I don't consider for the moment $$\hbar$$ and set $$m=1$$):

$$\begin{cases} \psi_t = \frac{\mathbf{i}}{2} \psi_{xx} \\ \psi(0,t) = \psi(1,t) = 0 \\ \psi(x,0) = \sin(2 \pi x) \end{cases}$$

where I choose as initial datum $$\psi(x,0)$$ a sinus. I integrate up to time $$T=1$$ using a suitable numerical method for the time integration and discretizing with finite difference in space as written above. I show in the following the plot of real and imaginary part at different times.

• Do you have any insist on reproducing that particular animation? I think it's better to solve a free particle in a box, where you have analytical solution and then compare your numerical results with that. If you agree I will elaborate it a bit more in an answer. Also, I don't think it's a good idea to convert Schrodinger equation, which is a PDE, to an ODE. I will elaborate a bit more about that as well in my answer if you want. – Alone Programmer Dec 19 '19 at 3:42
• Yes, if you could elaborate it as an answer it would be perfect ! Especially for the case where I have an analytical solution to compare. anyway, I prefer to keep my approach for the moment, since it's just a simple method of lines – VoB Dec 19 '19 at 6:46
• @AloneProgrammer I edited my answer considering a free particle in a box as you suggested. Is it okay in your opinion? I know that usually one end up with solving the time independent one, but I have to solve the time dependent. I don't know what could be the analytical solution, and moreover I don't know what am I supposed to plot once I found $\psi(x,t)$. I just plotted the real and imaginary part of the soluition at different times like $t=0.1,0.3$ – VoB Dec 19 '19 at 9:19
• Moreover, I noticed that the norm square of the solution at each time step is always equal to the norm of the initial data $\psi(x,0)= \sin(2 \pi x)$, which I know that it's a propery that has to be satisfied. – VoB Dec 19 '19 at 10:05
• It's a standard absorbing boundary condition. There are many references one can find for it regarding the Schroedinger equation. – Wolfgang Bangerth Dec 19 '19 at 20:06

## 2 Answers

Regarding the boundary conditions: Don't be fooled by Wikipedia. Yes, the scenario in the picture suggests an absorption at the boundaries, and yes, one could use absorbing boundary conditions in order to reproduce that numerically. In the simple case of a wavepacket these are readily available, because in the end, there exists an analytical solution for the wavefunction.

For more complex scenarios, however, absorbing boundary conditions are usually not that simple and the numerical schemes always entail approximations, i.e. they won't be completely reflectionless. In this answer, I tried to give a short overview on basically the simplest approach to absorbing boundaries (or better, absorbing boundary regions). In quantum mechanics, this often goes under the name complex absorbing potentials (CAPs).

So, if you simply want to reproduce this picture, just take a grid that is 10-times as large, and show only the inner region. The resulting picture will look exactly as the one from Wikipedia, but is obtained with a much smaller effort ... which is why I strongly believe the picture on Wikipedia has been produced in exactly this way.

Note further that the particle-in-a-box case is completely different to the wave equation, as this model necessarily implies reflections at the boundaries. It would absolutely make no sense to apply absorbing BCs here.

Initial conditions: The initial condition is a wavepacket. In general, this means you have an arbitrary mixture of plane waves $$e^{ikx}$$. A common example is to construct it by applying a another Gaussian to the plane wave on the basis of a plane-wave $$e^{i(kx)}$$ by applying a Gaussian with mean $$x_0$$ and width $$\lambda$$ to it, i.e. $$\Psi(x, t=0) = e^{i(kx)} e^{-\frac{(x-x_0)^2}{2\lambda^2}}$$ As usual, this should further be normalized so that $$||\Psi(x, t=0)||=1$$.

As David said, absorbing boundary conditions won't be completely reflectionless. That said, we can reduce relfections quite a bit, which helps to avoid influence from the boundaries while the particle still travelling inside.

Since this is a time dependent problem, one simple choice of boundary conditions will look like this.

• At the left boundary:

$$\frac{\partial \Psi}{\partial x}=-i k(E)\Psi =-i \sqrt{\frac{2mE}{\hbar^2}} \Psi$$

• At the right boundary:

$$\frac{\partial \Psi}{\partial x}=i k(E)\Psi =i \sqrt{\frac{2mE}{\hbar^2}} \Psi$$

What we are doing here, is telling our solution that it could pass the left boundary if it wants to travel in the negative x direction and it could pass the right boundary if it wants to travel in the positive direction.

As for the energy, we should clearly set:

$$E=\frac{i}{\hbar} \frac{\partial }{\partial t}$$

Now you may wonder how to extract the square root from a derivative, and one good choice is to approximate it:

$$\sqrt{\frac{2mE}{\hbar^2}} \approx \sqrt{\frac{2mE_0}{\hbar^2}} \frac{3E+E_0}{E+3E_0}$$

It's a good and well known rational approximation to square root.

In case of wavepacket, we can take:

$$E_0=\frac{\hbar^2 k_0^2}{2m}$$

Where $$k_0$$ is the initial momentum.

Now that we have a rational function, we can write for the left boundary.

$$\left(\frac{i}{\hbar} \frac{\partial }{\partial t}+3E_0 \right) \frac{\partial \Psi}{\partial x}=-i \sqrt{\frac{2mE_0}{\hbar^2}} \left(3\frac{i}{\hbar} \frac{\partial }{\partial t}+E_0 \right) \Psi$$

Now any proper finite difference scheme can be used with these boundary conditions in a very direct way.

Now, since we have taken an approximation for the square root, with an arbitrary parameter $$E_0$$, there's still going to be reflections. There are ways to improve this approximation, for example use rational functions of higher order.

Here's an illustrations on how these boundary conditions work. You can see a small part of the wavepacket being reflected at the right boundary, but there seems to be no reflection on the left. Note that it depends on both the group velocity and the spread of the packet. Here's a more impressive example. Note how we can clearly observe the two parts of the initial wavepacket - reflected and transmitted. If we change the boundaries to hard walls here, they would quickly reflect back at each other and make a mess. 