Test: Free Particle

In order to test our code, we first propagate a free particle, since we know the analytic result from the quantum dynamics lecture. As an initial form of the wavepacket we choose a Gaussian shape with some (arbitrary) width and initial momentum \(p_{0}\):

$$ \psi(x, 0) = \left[ 2\pi(\Delta x)_ {0}^{2} \right]^{\frac{1}{4}} \cdot \exp\left( \frac{-x^{2}}{4(\Delta x)_ {0}^{2}} + \frac{\mathrm{i}}{\hbar}p_{0}x\right ) $$

Analytically, the time propagation will be given by:

$$ \psi(x, t) = \left[ 2\pi(\Delta x)_ {0}^{2} \right]^{\frac{1}{4}} \cdot \left[1 + \frac{i\hbar t}{2m(\Delta x)_ {0}^{2}}\right]^{-\frac{1}{2}} \exp\left( \frac{\frac{-x^{2}}{4(\Delta x)_ {0}^{2}} + \frac{\mathrm{i}}{\hbar}p_{0}x

• \mathrm{i}\frac{p_{0}^{2}t}{2m\hbar}} {1 + \frac{\mathrm{i}\hbar t}{2m(\Delta x)_ {0}^{2}}} \right) $$

class FreeParticle(QuantumDynamics1D):

    def set_momentum(self, p0=1.0):
        self.p0 = p0
        
    def set_width(self, dx0=1.0):
        self.dx0 = dx0
        
    def get_wave_function(self, x, t):
        w = 1.0 + (1.0j*hbar*t) / (2.0*self.mass*self.dx0**2)
        ex = -x**2 / (4.0*self.dx0**2) \
            + 1.0j*(self.p0/hbar)*x \
            - 1.0j*(self.p0**2)*t / (2.0*self.mass*hbar)
        return (1.0 / np.sqrt(np.sqrt(2.0*np.pi*self.dx0**2))) \
                * (1.0/np.sqrt(w))*np.exp(ex/w)

    def propagate(self, time):
        nstep = time.shape[0]
        self.wavepackets = np.zeros((nstep, self.xgrid.shape[0]), dtype=complex)
        for i, t in enumerate(time):
            self.wavepackets[i, :] = self.get_wave_function(self.xgrid, t)         

m = 1.0
wp = FreeParticle(m)
wp.set_width()
wp.set_momentum()
xmin, xmax = -10.0, 60.0
ngrid = 128
tmin, tmax = 0.0, 20.0
nstep = 400

wp.setup_grid(xmin, xmax, ngrid)
time, dt = np.linspace(tmin, tmax, nstep, retstep=True)
wp.propagate(time)
anim_analytical = wp.animate(nstep=nstep, delay=50)

We can then compare the analytical result with the numerical one.

wp_num = QuantumDynamics1D(m)
wp_num.setup_grid(xmin, xmax, ngrid)
wp_num.set_potential(np.zeros(wp_num.xgrid.shape))
psi0 = wp.get_wave_function(wp_num.xgrid, 0.0)
wp_num.propagate(dt, nstep, psi0)
anim_numerical = wp_num.animate(nstep=nstep, delay=50)