Tijdsafhankelijke Schrödingervergelijking

De volledige tijdsafhankelijke Schrödingervergelijking behorende bij een Hamiltoniaan \(\hat{H}\) luidt:

\[ i\hbar\frac{\partial \Psi}{\partial t} = \hat{H}\Psi \]

Deze vergelijking gaat op voor ieder quantummechanisch systeem. In dit document beperken we ons tot een deeltje met massa \(m\) in een ééndimensionale potentiaal \(V(x)\). In dat geval krijgt de vergelijking de specifieke vorm:

\[ i\hbar\frac{\partial \Psi}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2\Psi}{\partial x^2}+V\Psi \]

De oplossingen \(\Psi(x,t)\) van deze vergelijking zijn functies van zowel de positie \(x\) als de tijd \(t\). Om het overzichtelijk te houden, maken we in eerste instantie gebruik van scheiding van variabelen. We zoeken dan naar oplossingen die geschreven kunnen als product van een ruimteafhankelijke functie en een tijdafhankelijke functie:

\[ \Psi(x,t) = \psi(x)\varphi(t) \]

Let op

Het is belangrijk je te realiseren dat dit niet noodzakelijk geldt voor elke oplossing. Sterker nog, voor veruit de meeste oplossingen geldt het niet! Toch is het handig om op deze manier te werk te gaan. Uiteindelijk blijkt namelijk dat elke oplossing kan worden geschreven als lineaire combinatie van een aantal van bovenstaande productoplossingen.

Voor het ruimteafhankelijke deel van de golffunctie \(\psi(x)\) kunnen we nu een aparte tijdsonafhankelijke Schrödingervergelijking opstellen:

\[ -\frac{\hbar^2}{2m}\frac{d^2\psi(x)}{dx^2}+V\psi(x) = E\psi(x)\ , \]

of simpelweg \(\hat{H}\psi=E\psi\). De oplossingen \(\psi_n\) van deze vergelijking zijn stationair in de tijd en hebben een goed gedefiniëerde energie \(E_n\). Deze energie-eigentoestanden vormen samen een complete set: elke mogelijke ruimtelijke golffunctie kan geschreven worden als lineaire combinatie van toestanden \(\psi_n\).

Let op

De oplossingen \(\psi_n\) worden vaak kortweg aangeduid als de 'eigentoestanden' van het systeem. Dat is strikt genomen onvolledig: het zijn de eigentoestanden van de operator \(\hat{H}\), en dus eigenlijk 'energie-eigentoestanden'. Gebruik in het vervolg gerust de term 'eigentoestanden', maar wees je ervan bewust dat elke observabele eigentoestanden heeft en dat het hier dus specifiek gaat om eigentoestanden van de Hamiltoniaan.

Ook voor het tijdsafhankelijke deel \(\varphi(t)\) kan een aparte vergelijking worden afgeleid:

\[ i\hbar\frac{d\varphi(t)}{dt} = E\varphi(t) \]

De oplossing van deze vergelijking is heel eenvoudig:

\[ \varphi(t) = e^{-iEt/\hbar} \]

Voegen we dit samen met de oplossing van de tijdsonafhankelijke Schrödingervergelijking, dan vinden we voor de energie-eigentoestanden:

\[ \Psi_n(x,t) = \psi_n(x)\varphi(t) = \psi_n(x)e^{-iE_nt/\hbar} \]

Dit resultaat kunnen we als volgt samenvatten:

Stelling

De energie-eigentoestanden van een quantummechanisch systeem hebben een tijdsafhankelijke complexe fasefactor \(e^{-iE_nt/\hbar}\). Deze fase draait met een snelheid die lineair schaalt met de energie van de toestand.

Is dit niet in strijd met de eerdere opmerking dat de energie-eigentoestanden stationair zijn in de tijd? Nee, want op het moment dat we de waarschijnlijkheidsdichtheid \(|\Psi_n(x,t)|^2\) bepalen, valt de fasefactor weg:

\[ |\Psi_n(x,t)|^2 = \psi_n^*(x)e^{+iE_nt/\hbar}\ \psi_n(x)e^{-iE_n/t/\hbar} = |\psi_n(x)|^2 \]

Een andere manier om hiernaar te kijken is dat energie altijd relatief is. De energie van een energie-eigentoestand is weliswaar bepaald, maar omdat we vrijheid hebben om het energie-nulpunt te kiezen, kunnen we de energie van alle toestanden samen altijd verschuiven met een willekeurige constante.

Dat betekent dat we aan een enkele energie-eigentoestand in principe iedere energie kunnen toekennen, dus ook \(E_n = 0\). Alles wat fysisch meetbaar is, zoals \(|\Psi|^2\), zou niet van deze keuze mogen afhangen. De factor \(e^{-iE_nt/\hbar}\) wordt daarom ook wel een globale fasefactor genoemd. Deze heeft geen fysische betekenis en kan altijd worden weggedeeld.

Stelling

Voor een enkele energie-eigentoestand heeft de tijdsafhankelijke fasefactor geen fysische betekenis. De bijbehorende waarschijnlijkheidsdichtheid is constant in de tijd.

Om deze stelling te visualiseren beschouwen we als voorbeeld de oneindig diepe potentiaalput. In onderstaande grafiek zie je hoe de eerste aangeslagen toestand \(\Psi_2\) variëert in de tijd, terwijl de waarschijnlijkheidsdichtheid \(|\Psi_2|^2\) constant blijft.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
plt.rcParams["animation.html"] = "jshtml"

n = 256
x = np.linspace(0,1,n)
y2 = np.sin(2*np.pi*x)

fig, axs = plt.subplots(2, 1, figsize=(5, 4))

axs[0].set_yticks([])
psi2, = axs[0].plot(x, y2.real, 'tab:blue')
psi2i, = axs[0].plot(x, y2.imag, 'tab:blue', alpha=0.5)
axs[0].legend(('Re [ $\Psi_2$ ]', 'Im [ $\Psi_2$ ]'), loc='lower left')
axs[0].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

axs[1].set_xlabel ('$x/a$')
axs[1].set_yticks([])
prob, = axs[1].plot(x, (y2*np.conj(y2)).real, 'tab:blue')
axs[1].fill_between(x, (y2*np.conj(y2)).real, color='tab:blue', alpha=0.25)
axs[1].legend(('$|\Psi_2|^2$',), loc='upper left')
axs[1].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

plt.close()

def animate(i):
    yt2 = np.exp(2*(i / 60)*np.pi*1j)*y2
    psi2.set_ydata(yt2.real)
    psi2i.set_ydata(yt2.imag)
    prob.set_ydata((yt2*np.conj(yt2)).real)

ani = animation.FuncAnimation(fig, animate, interval=40, frames=60)
ani

Once Loop Reflect

De situatie wordt iets lastiger als we gaan kijken naar oplossingen van de tijdsafhankelijke Schrödingervergelijking die niet overeenkomen met energie-eigentoestanden. De energie-eigentoestanden vormen een complete set, dus iedere willekeurige toestand kan worden geschreven als lineaire combinatie hiervan:

\[ \Psi(x,t) = \sum_{n=1}^\infty c_n\Psi_n(x,t) = \sum_{n=1}^\infty c_n\psi_n(x)e^{-iE_nt/\hbar} \]

In dit geval heeft iedere term een verschillende fasefactor. Deze kan dus niet worden weggedeeld. Omdat de fasefactoren met verschillende snelheden (energieën) draaien, bouwen de verschillende termen binnen de lineaire combinatie over de tijd een faseverschil op met elkaar. Termen die elkaar op het ene moment versterken, kunnen elkaar even later uitdoven.

Als voorbeeld kijken we opnieuw naar de oneindig diepe potentiaalput. Nu beschouwen we een lineaire combinatie tussen de grondtoestand \(\Psi_1\) en de eerste aangeslagen toestand \(\Psi_2\) (waarbij geldt dat \(E_2 = 4 E_1\)):

\[ \Psi(x,t) = \frac{1}{\sqrt{2}}\psi_1(x) e^{-iE_1t/\hbar} + \frac{1}{\sqrt{2}}\psi_2(x) e^{-iE_2t/\hbar} \]

In onderstaande grafiek ze we dat de waarschijnlijkheidsdichtheid nu wel degelijk variëert met de tijd. Het deeltje begint met grote waarschijnlijkheid in de linkerhelft van de put, waarna de waarschijnlijkheid om het deeltje aan te treffen heen en weer 'klotst' van links naar rechts.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
plt.rcParams["animation.html"] = "jshtml"

n = 256
x = np.linspace(0,1,n)
y1 = np.sin(1*np.pi*x)
y2 = np.sin(2*np.pi*x)

fig, axs = plt.subplots(3, 1, figsize=(5, 6))

axs[0].set_yticks([])
psi1, = axs[0].plot(x, y1.real, 'tab:brown')
psi2, = axs[0].plot(x, y2.real, 'tab:red')
axs[0].legend(('Re [ $\Psi_1$ ]','Re [ $\Psi_2$ ]'), loc='lower left')
axs[0].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

axs[1].set_ylim((-2,2))
axs[1].set_yticks([])
psi, = axs[1].plot(x, (y1 + y2).real, 'tab:blue')
psii, = axs[1].plot(x, (y1 + y2).imag, 'tab:blue', alpha=0.5)
axs[1].legend(('Re [ $\Psi_1+\Psi_2$ ]','Im [ $\Psi_1+\Psi_2$ ]'), loc='lower left')
axs[1].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

plt.close()

def animate(i):
    yt1 = np.exp(1*(i / 60)*np.pi*1j)*y1
    yt2 = np.exp(4*(i / 60)*np.pi*1j)*y2

    psi1.set_ydata(yt1.real)
    psi2.set_ydata(yt2.real)
    psi.set_ydata((yt1 + yt2).real)
    psii.set_ydata((yt1 + yt2).imag)

    axs[2].clear()
    axs[2].set_xlabel ('$x/a$')
    axs[2].set_ylim((-0.2,3.3))
    axs[2].set_yticks([])
    prob, = axs[2].plot(x, ((y1 + y2)*np.conj(y1 + y2)).real, 'tab:blue')
    prob.set_ydata(((yt1 + yt2)*np.conj(yt1 + yt2)).real)
    axs[2].fill_between(x, ((yt1 + yt2)*np.conj(yt1 + yt2)).real, color='tab:blue', alpha=0.25)
    axs[2].legend(('$|\Psi_1+\Psi_2|^2$',), loc='upper left')
    axs[2].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

ani = animation.FuncAnimation(fig, animate, interval=40, frames=120)
ani

Once Loop Reflect

Hoewel we de fasefactoren niet helemaal kunnen elimineren, kan het wel inzichtelijk zijn om één van de fasefactoren (bijvoorbeeld die van de grondtoestand) als globale fasefactor weg te delen:

\[ \Psi(x,t) = \frac{1}{\sqrt{2}}\psi_1(x) + \frac{1}{\sqrt{2}}\psi_2(x) e^{-i(E_2-E_1)t/\hbar} \]

Het gevolg is dat nu één van de twee golffuncties stilstaat, waardoor het makkelijker is om de interferentiecriteria te volgen over de tijd. Merk op dat dit voor de waarschijnlijkheidsdichtheid geen verschil maakt. Deze mag immers nooit afhangen van de keuze van het energie-nulpunt.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
plt.rcParams["animation.html"] = "jshtml"

n = 256
x = np.linspace(0,1,n)
y1 = np.sin(1*np.pi*x)
y2 = np.sin(2*np.pi*x)

fig, axs = plt.subplots(3, 1, figsize=(5, 6))

axs[0].set_yticks([])
psi1, = axs[0].plot(x, y1.real, 'tab:brown')
psi2, = axs[0].plot(x, y2.real, 'tab:red')
axs[0].legend(('Re [ $\Psi_1$ ]','Re [ $\Psi_2$ ]'), loc='lower left')
axs[0].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

axs[1].set_ylim((-2,2))
axs[1].set_yticks([])
psi, = axs[1].plot(x, (y1 + y2).real, 'tab:blue')
psii, = axs[1].plot(x, (y1 + y2).imag, 'tab:blue', alpha=0.5)
axs[1].legend(('Re [ $\Psi_1+\Psi_2$ ]','Im [ $\Psi_1+\Psi_2$ ]'), loc='lower left')
axs[1].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

plt.close()

def animate(i):
    yt1 = np.exp(0*(i / 60)*np.pi*1j)*y1
    yt2 = np.exp(3*(i / 60)*np.pi*1j)*y2

    psi1.set_ydata(yt1.real)
    psi2.set_ydata(yt2.real)
    psi.set_ydata((yt1 + yt2).real)
    psii.set_ydata((yt1 + yt2).imag)

    axs[2].clear()
    axs[2].set_xlabel ('$x/a$')
    axs[2].set_ylim((-0.2,3.3))
    axs[2].set_yticks([])
    prob, = axs[2].plot(x, ((y1 + y2)*np.conj(y1 + y2)).real, 'tab:blue')
    prob.set_ydata(((yt1 + yt2)*np.conj(yt1 + yt2)).real)
    axs[2].fill_between(x, ((yt1 + yt2)*np.conj(yt1 + yt2)).real, color='tab:blue', alpha=0.25)
    axs[2].legend(('$|\Psi_1+\Psi_2|^2$',), loc='upper left')
    axs[2].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

ani = animation.FuncAnimation(fig, animate, interval=40, frames=120)
ani

Once Loop Reflect

Dit brengt ons tot de volgende conclusie:

Stelling

In een lineaire combinatie van energie-eigentoestanden bouwen de verschillende termen een faseverschil op met elkaar. Dit faseverschil draait met een snelheid die schaalt met het energieverschil tussen de termen. De waarschijnlijkheidsdichtheid fluctueert hierdoor als functie van de tijd.

Let op

Deze tijdevolutie geldt voor een Hamiltoniaan die zelf onafhankelijk is van de tijd. Beschouw je een Hamiltoniaan die wel afhangt van de tijd, dan is het dus belangrijk onderscheid te maken tussen tijdevolutie die voortkomt uit de tijdsafhankelijkheid van de Hamiltoniaan enerzijds, en de hier beschreven intrinsieke tijdevolutie anderzijds.

Het spreekt voor zich dat het aantal mogelijke lineaire combinaties eindeloos is. Mocht je zelf verder willen experimenteren, gebruik dan gerust de Pythoncode van bovenstaande grafiek:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
plt.rcParams["animation.html"] = "jshtml"

n = 256
x = np.linspace(0,1,n)
y1 = np.sin(1*np.pi*x)
y2 = np.sin(2*np.pi*x)

fig, axs = plt.subplots(3, 1, figsize=(5, 6))

axs[0].set_yticks([])
psi1, = axs[0].plot(x, y1.real, 'tab:brown')
psi2, = axs[0].plot(x, y2.real, 'tab:red')
axs[0].legend(('Re [ $\Psi_1$ ]','Re [ $\Psi_2$ ]'), loc='lower left')
axs[0].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

axs[1].set_ylim((-2,2))
axs[1].set_yticks([])
psi, = axs[1].plot(x, (y1 + y2).real, 'tab:blue')
psii, = axs[1].plot(x, (y1 + y2).imag, 'tab:blue', alpha=0.5)
axs[1].legend(('Re [ $\Psi_1+\Psi_2$ ]','Im [ $\Psi_1+\Psi_2$ ]'), loc='lower left')
axs[1].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

plt.close()

def animate(i):
    yt1 = np.exp(0*(i / 60)*np.pi*1j)*y1
    yt2 = np.exp(3*(i / 60)*np.pi*1j)*y2

    psi1.set_ydata(yt1.real)
    psi2.set_ydata(yt2.real)
    psi.set_ydata((yt1 + yt2).real)
    psii.set_ydata((yt1 + yt2).imag)

    axs[2].clear()
    axs[2].set_xlabel ('$x/a$')
    axs[2].set_ylim((-0.2,3.3))
    axs[2].set_yticks([])
    prob, = axs[2].plot(x, ((y1 + y2)*np.conj(y1 + y2)).real, 'tab:blue')
    prob.set_ydata(((yt1 + yt2)*np.conj(yt1 + yt2)).real)
    axs[2].fill_between(x, ((yt1 + yt2)*np.conj(yt1 + yt2)).real, color='tab:blue', alpha=0.25)
    axs[2].legend(('$|\Psi_1+\Psi_2|^2$',), loc='upper left')
    axs[2].axhline(y=0,linewidth=0.5, linestyle='dashed', color='grey')

ani = animation.FuncAnimation(fig, animate, interval=40, frames=120)
ani