10 Closed Quantum Systems in QuTiP: The Schrödinger equation
10.1 From wavefunctions to the Schrödinger equation in configuration space
In Chapter 7 we moved from Hamilton’s equations for a single point to Liouville’s equation for an ensemble. Quantum mechanics follows a similar path, replacing phase-space densities with wavefunctions (or density operators). We now show this parallelism in detail.
10.1.1 From a single wavepacket to the quantum continuity law
A non-relativistic quantum system in one dimension is described by a complex wavefunction \(\psi(x,t)\). Its probability density and current are
\[ \rho(x,t) \;=\; |\psi(x,t)|^2, \qquad j(x,t) \;=\; \frac{\hbar}{m}\,\Im\!\bigl[\psi^{\!*}(x,t)\,\partial_x \psi(x,t)\bigr]. \]
Starting from the time-dependent Schrödinger equation
\[ i\hbar\,\partial_t \psi(x,t) \;=\; \hat H\,\psi(x,t) \quad\text{with}\quad \hat H \;=\; -\frac{\hbar^2}{2m}\,\partial_x^2 + V(x) \tag{10.1}\]
and its complex conjugate, a short calculation yields the continuity equation
\[ \partial_t \rho + \partial_x j \;=\; 0, \tag{10.2}\]
which plays the same role as Liouville’s incompressibility condition: probability is transported in configuration space without being created or destroyed.
10.1.2 The Schrödinger equation as a linear ODE in Hilbert space
While Equation 10.1 is a partial differential equation, it is linear in \(\psi\), exactly like the matrix ODE \(\dot{\mathbf y}=A\mathbf y\) used for classical linear systems. Written in Dirac notation
\[ i\hbar\,\frac{d}{dt}\,|\psi(t)\rangle \;=\; \hat H\,|\psi(t)\rangle, \]
the resemblance becomes explicit: replace the classical matrix \(A\) by the quantum Hamiltonian operator \(\hat H\) and the state vector \(\mathbf y\) by a ket \(|\psi\rangle\).
For statistical mixtures one introduces the density operator \(\hat\rho\) and obtains the von Neumann equation
\[ i\hbar\,\partial_t \hat\rho \;=\; [\hat H,\,\hat\rho], \tag{10.3}\]
a direct quantum counterpart of the Liouville equation in Equation 7.5.
With these ingredients we now possess a one-to-one map between the classical Liouville formulation and the quantum Schrödinger formulation, both expressible as sparse linear ODEs ready for numerical treatment. Let’s consider a simple example of a harmonic oscillator, which is described by the Hamiltonian
\[ \hat{H} = \hbar \omega \hat{a}^\dagger \hat{a} \, , \]
where \(\hat{a}\) and \(\hat{a}^\dagger\) are the annihilation and creation operators, respectively, and \(\omega\) is the angular frequency of the oscillator.
In QuTip, the Schrödinger equation is solved with the function sesolve, which stands for “Schrödinger equation solver”.
10.2 Plot the expectation values
We can access to the expectation values with the command sol.expect[i], where i is the index of the \(i\)-th operator for which we want to calculate te expectation value as a function of time.
fig, ax = plt.subplots()
ax.plot(tlist, sol.expect[0], label=r"$\langle \hat{H} \rangle$", lw=2)
ax.plot(tlist, sol.expect[1], label=r"$\langle \hat{a} + \hat{a}^\dagger \rangle$", lw=2)
ax.legend()
ax.set_xlabel(r"$t$")
ax.set_xlim(tlist[0], tlist[-1])
ax.set_ylim(None, 16)
# Show in Quarto
plt.savefig("_tmp_fig.svg")
plt.close(fig)
SVG("_tmp_fig.svg")
10.3 Access directly to the evolution of the state
We can also access directly to the wavefunction at each tlist. This can be done by simply calling sol.states to the solver without the e_ops operators.
To check this feature, let’s see if after \(10\) cycles we get still the initial state, by calculating the fidelity
\[ \langle \psi \left(t = 10 T\right) \vert \psi \left( t = 0 \right) \rangle \]
where \(T = \frac{2 \pi}{\omega}\).
We then switch to the position eigenbasis. Thus, we need to diagonalize the position operator. This in general involves the Laguerre functions, but here we limit ourselves to numerically diagonalize the position operator.
We now show the wavefunctions in the position basis at three different times: \(t_0\), \(t_1\) and \(t_2\), with
\[ t_0 = 0 \, , \quad t_1 = \frac{\pi}{\omega} \, , \quad t_2 = \frac{2 \pi}{\omega} \, , \]
showing the exact periodicity of this system.
idx_t0 = 0
idx_t1 = np.where(tlist >= np.pi/w)[0][0]
idx_t2 = np.where(tlist >= 2*np.pi/w)[0][0]
psi0_x = U.dag() * sol.states[idx_t0]
psi1_x = U.dag() * sol.states[idx_t1]
psi2_x = U.dag() * sol.states[idx_t2]
# We define the potential for the harmonic oscillator
V = 0.5 * m * w**2 * xlist**2
fig, ax = plt.subplots()
ax.plot(xlist, 700 * np.abs(psi0_x.full())**2,
label=r"$\vert \langle \psi (t_0) \vert \psi (t_0) \rangle \vert^2$", lw=2)
ax.plot(xlist, 700 * np.abs(psi1_x.full())**2,
label=r"$\vert \langle \psi (t_1) \vert \psi (t_1) \rangle \vert^2$", ls="--", lw=2)
ax.plot(xlist, 700 * np.abs(psi2_x.full())**2,
label=r"$\vert \langle \psi (t_2) \vert \psi (t_2) \rangle \vert^2$", ls="-.", lw=2)
ax.plot(xlist, V, color="black", ls="--")
ax.legend()
ax.set_xlabel(r"$x$")
ax.set_xlim(xlist[0], xlist[-1])
ax.set_ylim(0, 120)
# Show in Quarto
plt.savefig("_tmp_fig.svg")
plt.close(fig)
SVG("_tmp_fig.svg")
We can also export an animation, showing more easily the time evolution of the state
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
plt.rcParams.update({'font.size': 8})
fig, ax = plt.subplots(figsize=(4.6, 2.8))
line, = ax.plot(xlist, 700 * np.abs(psi0_x.full())**2, lw=2)
ax.plot(xlist, V, color="black", ls="--")
ax.set_xlabel(r"$x$")
ax.set_xlim(xlist[0], xlist[-1])
ax.set_ylim(0, 120)
plt.close(fig) # Otherwise the static figure also appears
def update(frame):
psi_t = U.dag() * sol.states[frame]
line.set_ydata( 700 * np.abs(psi_t.full())**2 )
return line,
fps = 25
ani = FuncAnimation(fig, update, frames=len(tlist), blit=True, interval=1000 / fps)
HTML(ani.to_jshtml())