Time Independent Schroedinger Equation
In this tutorial, we're gonna solve the time independent schroedinger equation for a non trivial potential system.
import numpy as np
from matplotlib.pyplot import subplots
np.set_printoptions(precision=4, suppress=True)
from quantized.basis import HarmonicOscillator, EigenBasis
from quantized.operators import Hamiltonian, Overlap
Using Arbitrary Potentials and a Basis¶
A good example to start is a symmetric double well potential. Let's create one. This also demonstrates the ability of the library to work on arbitrary potential functions
def potential(x):
return -0.5*x**2 + 0.04 * x**4 + 2
fig, ax = subplots()
x = np.linspace(-4, 4, 300)
ax.plot(x, potential(x))
ax.set_xlabel("x")
_ = ax.set_ylabel("V(x)")
We can't solve this exactly, so we need to introduce a basis set. A reasonable choice would be to use a few harmonic oscillator functions in each of the two wells.
basis = [
HarmonicOscillator(n=0, center=-2),
HarmonicOscillator(n=1, center=-2),
HarmonicOscillator(n=0, center=2),
HarmonicOscillator(n=1, center=2),
]
fig, ax = subplots()
ax.plot(x, potential(x), "k--")
ax.set_xlabel("x")
ax.plot(x, basis[0](x) + basis[0].energy, color="C0")
ax.plot(x, basis[1](x) + basis[1].energy, color="C0")
ax.plot(x, basis[2](x) + basis[2].energy, color="C1")
_ = ax.plot(x, basis[3](x) + basis[3].energy, color="C1")
And that's basically it! We have a basis ready to go, and we've got our potential. Let's get the Hamiltonian!
Getting the Hamiltonian and Overlap Matrices¶
H = Hamiltonian(potential).matrix(basis)
H
Normally, we'd diagonalize the hamiltonian and be done. However, we have a non-orthogonal basis, so we need to account for that. We have to find an eigenbasis of both the hamiltonian and overlap operators. This means, we need the overlap matrix, and we need to solve the generalized eigenvalue problem, $\mathbf{HC} = \mathbf{ESC}$
S = Overlap().matrix(basis)
S
As you can see, we're justified in needing to diagonalize the overlap matrix. The basis functions within each well are indeed orthonormal, but the ones in separate wells are not.
The Eigenstates¶
Fortunately, it's easy to get the eigenstates from the $\mathbf{H}$ and $\mathbf{S}$ matrices:
eig_basis = EigenBasis.from_basis(basis, H, S)
print(f"EigenEnergies = {eig_basis.energies}")
print(f"EigenVectors = \n{eig_basis.eigvecs}")
The eigen states can be represented as linear combinations of the original basis.
$\phi(x)_j = \sum_i \chi_i \mathbf{C}_{ji}$
Where $\phi(x)_j$ is the $j$'th eigen state, $\chi_i$ is the $i$'th basis state, and $\mathbf{C}_{ji}$ is the $i, j$'th element of the eigenvector matrix.
The eigenstates are available from EigenBasis
for i, state in enumerate(eig_basis.states):
print(f"State {i}\n")
print(state)
print("-"*80)
And each of these eigenstates is available as a plain old function in one dimension
fig, ax = subplots()
ax.plot(x, potential(x), "k--")
ax.set_xlabel("x")
for state, e in zip(eig_basis.states, eig_basis.energies):
ax.plot(x, state(x) + e)
There we have it! The solutions to the double well potential. In real systems we will want to use a larger basis set, but this suffices for demonstration purposes.
Transforming from Original to Eigenbasis¶
Often it will be considerably easier to compute a quantity in the original basis. It can be transformed into the eigenbasis by a similarity transform,
$ \mathbf{X}_{eigen} = \mathbf{C}^{-1} \mathbf{S}^{-1} \mathbf{X}_{original} \mathbf{C}$
Where $\mathbf{S}^{-1}$ is the overlap matrix and $ \mathbf{C}$ is the eigenvector matrix. For example, we may want to validate that our eigenbasis diagonalizes the hamiltonian matrix:
eig_basis.transformed(H)
We also might want to make sure that we've diagonalized the overlap matrix as well:
eig_basis.transformed(S)
Remember, this transormation can be done to any matrix. And before we've shown how we can transform any of our Operators into matrix form. So if we make an operator, we can easily get its form in the eigen basis. To demonstrate, we can get the overlap in the $-x$ portion of the potential with ease:
eig_basis.transformed(Overlap(-np.inf, 0).matrix(basis))