Time Dependent Schroedinger Equation

This wouldn't be a quantum dynamics simulation without some time dependence!

In [1]:
from IPython.display import HTML
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
np.set_printoptions(precision=4, suppress=True)

from quantized.basis import HarmonicOscillator, EigenBasis
from quantized.operators import Hamiltonian, Overlap

from quantized.time_evolution import TimeEvolvingState, TimeEvolvingObservable

Time independent Solution

Hoisting the time independent solution from the previous tutorial:

In [2]:
def potential(x):
    return -0.5*x**2 + 0.04 * x**4 + 1

basis = [
    HarmonicOscillator(n=0, center=-2),
    HarmonicOscillator(n=1, center=-2), 
    HarmonicOscillator(n=0, center=2), 
    HarmonicOscillator(n=1, center=2),
]

H = Hamiltonian(potential).matrix(basis)
S = Overlap().matrix(basis)

eig_basis = EigenBasis.from_basis(basis, H, S)
eig_basis
Out[2]:
EigenBasis(states=(0.6992 * HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0)
 + -0.0869 * HarmonicOscillator(n=1, center=-2, mass=1.0, omega=1.0)
 + 0.6992 * HarmonicOscillator(n=0, center=2, mass=1.0, omega=1.0)
 + 0.0869 * HarmonicOscillator(n=1, center=2, mass=1.0, omega=1.0), 0.6967 * HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0)
 + -0.1280 * HarmonicOscillator(n=1, center=-2, mass=1.0, omega=1.0)
 + -0.6967 * HarmonicOscillator(n=0, center=2, mass=1.0, omega=1.0)
 + -0.1280 * HarmonicOscillator(n=1, center=2, mass=1.0, omega=1.0), -0.0577 * HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0)
 + -0.6608 * HarmonicOscillator(n=1, center=-2, mass=1.0, omega=1.0)
 + -0.0577 * HarmonicOscillator(n=0, center=2, mass=1.0, omega=1.0)
 + 0.6608 * HarmonicOscillator(n=1, center=2, mass=1.0, omega=1.0), 0.1599 * HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0)
 + 0.7476 * HarmonicOscillator(n=1, center=-2, mass=1.0, omega=1.0)
 + -0.1599 * HarmonicOscillator(n=0, center=2, mass=1.0, omega=1.0)
 + 0.7476 * HarmonicOscillator(n=1, center=2, mass=1.0, omega=1.0)), energies=(1.112830250685375, 1.133287513512955, 2.114702264414521, 2.4540685273585483), ao_S=array([[ 1.    ,  0.    ,  0.0183, -0.0518],
       [ 0.    ,  1.    ,  0.0518, -0.1282],
       [ 0.0183,  0.0518,  1.    ,  0.    ],
       [-0.0518, -0.1282,  0.    ,  1.    ]]), eigvecs=array([[ 0.6992,  0.6967, -0.0577,  0.1599],
       [-0.0869, -0.128 , -0.6608,  0.7476],
       [ 0.6992, -0.6967, -0.0577, -0.1599],
       [ 0.0869, -0.128 ,  0.6608,  0.7476]]), ao_basis=[HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0), HarmonicOscillator(n=1, center=-2, mass=1.0, omega=1.0), HarmonicOscillator(n=0, center=2, mass=1.0, omega=1.0), HarmonicOscillator(n=1, center=2, mass=1.0, omega=1.0)])

Time evolving states

To get a state evolving in time, we need the eigenbasis and an initial state.

In [3]:
initial = HarmonicOscillator(n=0, center=-2)
time_state = TimeEvolvingState(initial, eig_basis)
time_state
Out[3]:
TimeEvolvingState(initial_state=HarmonicOscillator(n=0, center=-2, mass=1.0, omega=1.0))

Time evolving states are simply functions of x and t, namely:

In [4]:
for t in range(4):
    val = time_state(0.0, t)
    print(f"State's wavefunction has value {val} at time {t}")
State's wavefunction has value (0.1016537883064175+0j) at time 0
State's wavefunction has value (0.009977734082449926-0.08967532301214656j) at time 1
State's wavefunction has value (-0.05664250439519702-0.01945510321670284j) at time 2
State's wavefunction has value (-0.027585791118214437+0.010539634672632739j) at time 3

To get a better feel for this, we can plot the wave function over time:

In [5]:
%%capture  
# matplotlib nonsense

from quantized.plotting import animate_state

fig, ax = plt.subplots()
ax.set_xlim(-4, 4)
ax.set_ylim(-1, 4)

a = animate_state(fig, ax, initial, time_state, potential, nframes=200, interval=20)
In [6]:
HTML(a.to_html5_video())
Out[6]:

Here we can even get a glimpse at some tunneling! If you squint really hard, you can see some movement in the right well.

Now would be an interesting time to look at what happens if we change up the initial state.

In [7]:
%%capture
# matplotlib nonsense again
fig, ax = plt.subplots()
ax.set_xlim(-4, 4)
ax.set_ylim(-1, 4)

initial_state_2 = HarmonicOscillator(n=1, center=-2)
time_state_2 = TimeEvolvingState(initial_state_2, eig_basis)
anim = animate_state(fig, ax, initial_state_2, time_state_2, potential, nframes=200, interval=20)
In [8]:
HTML(anim.to_html5_video())
Out[8]:

Whoa! A big difference! Notice how the state is much more mobile. Since we started off in the n=1 state rather than the n=2, it has much more energy to cross the barrier.

Time evolving Observables


So we've seen how to look at the time evolution of some initial state in a potential. Next we want to look at how observables change over time. Let's use the second of our initial states, it should be more interesting. If we look at the time evolution of the density in the left and right wells, we can get a feel for the oscillation of the particle's density over time.

In [9]:
positive_density = time_state_2.observable(Overlap(lower_limit=0.0), hermitian=True)
neg_density = time_state_2.observable(Overlap(upper_limit=0.0), hermitian=True)

fig, ax = plt.subplots()

times = np.linspace(0, 40, 400)
ax.plot(times, [neg_density(t) for t in times], label="Negative x")
ax.plot(times, [positive_density(t) for t in times], label="Positive x")

_ = ax.legend()

Notice how the density oscillates between the left well and the right well. The period is roughly 20 atomic units of time, which if you look again at the second animation above (which was animated for 40 atu), you can see the pattern repeating itself.

Since this is the type of thing we're doing in the analysis, it's been made easy to do. The code below does the same exact thing we just did.

In [10]:
from quantized.hopping_matrix import OccupancyProbabilites

occprobs = OccupancyProbabilites.from_1d_state(time_state_2, borders=(-np.inf, 0.0, np.inf))

fig, ax = plt.subplots()
occ_probs_over_time = zip(*[occprobs(t) for t in times])
for i, ot in enumerate(occ_probs_over_time):
    ax.plot(times, ot, label=f"Region={i}")
_ = ax.legend()