Obligatorisk øving i TFY4215 april 2018¶

Løsningsforslag til oppgave 1¶

Vi starter som vanlig med å importere numpy, samt matplotlib for å lage figurer:

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import display, Image, HTML
print(animation.writers.list())

['ffmpeg', 'ffmpeg_file', 'html']

Vi trenger $\hbar$ og elektronmassen $m$ ; og bruker ganske enkelt SI-enheter:

hbar=1.05E-34
m=9.11E-31

Velger tallverdier for (det konstante) potensialet $V_0 = 0$ , steglengden $dx = 1$ Å, samt bølgepakkens senterenergi $E_0 = 1.0$ eV:

V0=0.0
dx=1.0E-10
E0=1.0*1.6E-19

Fra $E_0$ beregnes tilhørende bølgetall $k_0$ , impuls $p_0$ og hastighet $v_0$ :

k0=np.sqrt(2.0*m*E0)/hbar
p0=hbar*k0
v0=p0/m

Vi gir potensialboksen bredde $L = N \cdot dx$ , som betyr at $V = 0$ på $N$ punkter, her med $N = 1000$ . Deretter konstrueres Hamilton-matrisen, med diagonalelementer $V_i + \hbar^2/m (dx)^2$ og med $H(i, i \pm 1) = - \hbar^2/2 m (dx)^2$ . Egenverdier og egenvektorer til matrisen $H$ returneres fra numpy-funksjonen linalg.eigh(H):

#N = antall punkter langs x-aksen
N = 500
#Potensial for fri partikkel, V=0 i alle punkter
V = [V0]*N
#Lager numpy array
V = np.asarray(V)
#Potensialet i enhet eV
V_eV = V/1.6E-19
#d = liste med diagonalelementer i Hamiltonmatrisen H
d = [v + hbar**2/(m*dx**2) for v in V]
#e = verdi til ikke-diagonale elementer i H, dvs H(i,i+-1) = -e
e = - hbar**2/(2*m*dx**2)
#Ntot = antall elementer i V
Ntot=len(V)
#Initialisering av matrisen H: Legger inn verdi 0 i samtlige elementer
H = [[0]*(Ntot) for n in range(Ntot)]
#Dobbel for-løkke som lager den tridiagonale matrisen H
for i in range(Ntot):
    for j in range(Ntot):
        if i==j:
            H[i][j]=d[i]
        if abs(i-j)==1:
            H[i][j]=e

#Gjør om H fra liste til numpy array (muligens uten betydning?)
H = np.asarray(H)
#Finner energy = egenverdiene og psi_matrix = egenvektorene til matrisen H           
energy,psi_matrix = np.linalg.eigh(H)
#energy_eV = liste med energiegenverdier i enheten eV
energy_eV = energy/1.6E-19

Vi kan nå plotte de beregnede energiegenverdiene som funksjon av "nummer" $n$ . Vi vet fra før at både bølgetallet og impulsen til partikkelen er proporsjonal med (kvantetallet) $n$ , så denne figuren illustrerer essensielt "dispersjonsrelasjonen" $E(p)$ .

n_axis = [n for n in range(Ntot)]
plt.figure('Energiegenverdier')
#plt.plot(n_axis,energy_eV,'bo')
plt.plot(n_axis,energy_eV)
plt.title('Energiegenverdier',fontsize=16)
plt.xlabel('$n$',fontsize=16)
plt.ylabel('Energi (eV)',fontsize=16)
plt.show()

Her er det allerede et par saker en kan merke seg: For det første ser vi at vi har fått like mange energiegenverdier som vi har punkter langs $x$ -aksen, eller som størrelsen $N$ på Hamilton-matrisen. Dette var selvsagt ikke så overraskende: En $N \times N$ -matrise vil ha $N$ egenverdier, og vi visste fra før at vi ikke har degenerasjon i et slik endimensjonalt system. For det andre ser vi at $E(p)$ er kvadratisk for små verdier av $n$ (dvs små verdier av $p$ ), som for en fri partikkel. Men for økende verdier av $n$ ser vi at det blir stadig større avvik fra en kvadratisk dispersjonsrelasjon. Hele kurven $E(n)$ minner mistenkelig om en halv periode av en harmonisk funksjon, og det er nettopp det den er. Vi har fått et energibånd med tillatte stasjonære tilstander, mellom $E=0$ og $E = 2.4204 \cdot 10^{-18}$ J = 15.128 eV. Det er ingen tilfeldighet at denne båndbredden tilsvarer nøyaktig fire ganger absoluttverdien av de ikke-diagonale elementene $-\hbar^2/2m(dx)^2$ i Hamilton-matrisen, men vi overlater til emnet TFY4220 Faste stoffers fysikk å vise hvorfor det er nettopp slik! Her konkluderer vi i første omgang med at dersom vi ønsker å beskrive en tilnærmet fri partikkel, bør bølgepakkens (senter-)energi velges så liten at den befinner seg i området med kvadratisk dispersjon. Ovenfor valgte vi $E_0 = 1.0$ eV, noe som ser ut til å være brukbart inne i området med kvadratisk $E(n)$ .

Vi lager i neste omgang starttilstanden $\Psi(x,0)$ som en "minimum uskarphet" (normert) gaussisk bølgepakke:

Ψ (x, 0) = \frac{1}{(2 π σ^{2})^{1 / 4}} \exp [- \frac{(x - x_{0})^{2}}{4 σ^{2}} + i k_{0} x] .

$\begin{equation} \Psi(x,0) = \frac{1}{(2\pi\sigma^2)^{1/4}}\,\exp\left[-\frac{(x-x_0)^2}{4\sigma^2} + i k_0 x\right]. \end{equation}$ Med denne funksjonen har starttilstanden et standardavvik i posisjon lik

σ

$\sigma$ , "tyngdepunktet" i posisjon

x_{0}

$x_0$ , og en senterimpuls

p_{0} = ℏ k_{0}

$p_0 = \hbar k_0$ . I eksemplet nedenfor brukes tallverdiene

σ = 50 d x

$\sigma = 50 dx$ og

x_{0} = 10 σ

$x_0 = 10 \sigma$ :

#x = liste med posisjonsverdier (i enheten m)
x = np.asarray([dx*n for n in range(Ntot)])
#x_nm = liste med posisjonsverdier i enheten nanometer
x_nm = x * 1.0E9
#INITIALTILSTANDEN Psi(x,0)
#sigma = initialtilstandens standardavvik i posisjon
sigma = dx*30.0
#x0 = pakkens senterposisjon; start minst 5*sigma til høyre for x=0
x0 = 5.0*sigma
#
#Lager starttilstanden Psi(x,0)
normfactor = (2*np.pi*sigma**2)**(-0.25)
gaussinit = np.exp(-(x-x0)**2/(4*sigma**2))
planewavefactor = np.exp(1j*k0*x)
Psi0 = normfactor*gaussinit*planewavefactor
#Plotter sannsynlighetstettheten for initialtilstanden
plt.figure('Initiell sannsynlighetstetthet')
plt.plot(x_nm,np.abs(Psi0**2))
plt.title('Initiell sannsynlighetstetthet',fontsize=16)
plt.xlabel('$x$ (nm)',fontsize=16)
plt.ylabel('$|\Psi(x,0)|^2$',fontsize=16)
plt.show()

Nå kan vi beregne utviklingskoeffisientene $c_n$ i utviklingen av $\Psi(x,0)$ i stasjonære tilstander:

Ψ (x, 0) = \sum_{n = 0}^{N - 1} c_{n} ψ_{n} (x),

$\begin{equation} \Psi(x,0) = \sum_{n=0}^{N-1} c_n \psi_n(x), \end{equation}$ som betyr at

c_{n} = \int_{0}^{L} ψ_{n}^{*} (x) Ψ (x, 0) d x .

$\begin{equation} c_n = \int_0^L \psi_n^*(x)\,\Psi(x,0)\,dx. \end{equation}$

#Lager Psi(x,0) som lineærkombinasjon av egentilstander
#Lager utviklingskoeffisientene c[0] ... c[Ntot-1]
#Omgjør egenfunksjonene fra reelt til komplekst representerte tall
psi_matrix_complex = psi_matrix*(1.0 + 0.0j)
c = np.zeros(Ntot, dtype = np.complex128)
for n in range(Ntot):
    c[n] = np.vdot(psi_matrix_complex[:,n],Psi0)    
#Plotter |c[n]|**2
#Dette representerer på et vis impulsfordelingen i bølgepakken
n_axis = [n for n in range(Ntot)]
plt.figure('Utviklingskoeffisientene')
plt.plot(n_axis,np.abs(c)**2)
plt.title('Utviklingskoeffisientene',fontsize=16)
plt.xlabel('$n$',fontsize=16)
plt.ylabel('$|c(n)|^2$',fontsize=16)
plt.show()

Som en sjekk på at programmet fungerer som det skal, rekonstruerer vi starttilstanden med utgangspunkt i de beregnede utviklingskoeffisientene $c_n$ . Deretter plotter vi både den "originale" og den rekonstruerte sannsynlighetstettheten $|\Psi(x,0)|^2$ i samme figur.

#Rekonstruerer Psi(x,0) fra utviklingskoeffisientene
Psi0rekon = np.zeros(Ntot,dtype=np.complex128)
for n in range(Ntot):
    Psi0rekon = Psi0rekon + c[n]*psi_matrix_complex[:,n]

#Plotter sannsynlighetstetthet for rekonstruert initialtilstand
plt.figure('Rekonstruert og original initialtilstand')
plt.plot(x_nm,np.abs(Psi0rekon**2),x_nm,np.abs(Psi0**2))
plt.title('Rekonstruert og original initialtilstand',fontsize=16)
plt.xlabel('$x$ (nm)',fontsize=16)
plt.ylabel('$|\Psi(x,0)|^2$',fontsize=16)
plt.show()

Fornøyd noterer vi oss at de to kurvene overlapper perfekt.

Vi velger et passende tidssteg og et passende tidsintervall, beregner $\Psi(x,t)$ , og deretter tidsutviklingen av bølgepakkens standardavvik, $\Delta x(t)$ :

#Setter tidssteget til 1/100 av tiden det tar for bølgepakken å
#flytte seg en lengde lik lengden av hele x-området
tidssteg = N*dx/(v0*100)
#Initialiserer x2mean=<x**2> og xmean2=<x>**2 og tidsaksen
Ntid = 500
x2mean = np.zeros(Ntid)
xmean2 = np.zeros(Ntid)
tid = np.zeros(Ntid)
deltax_analytisk = np.zeros(Ntid)
#Lar bølgepakken kollidere et par ganger med de harde veggene
for i in range(Ntid):
    tid[i] = i*tidssteg
    deltax_analytisk[i] = np.sqrt(sigma**2+hbar**2*tid[i]**2/(4*m**2*sigma**2))
    #Lager Psi(x,t). Nullstiller først 
    Psi_t = np.zeros(Ntot,dtype=np.complex128)
    for n in range(Ntot):
        Psi_t = Psi_t + c[n]*psi_matrix_complex[:,n]*np.exp(-1j*energy[n]*tid[i]/hbar)
    
    rho_t = np.abs(Psi_t)**2
    x2mean[i] = dx*np.dot(x**2,rho_t)
    xmean2[i] = dx*dx*np.dot(x,rho_t)**2

deltax = np.sqrt(x2mean-xmean2)

plt.figure('Standardavvik i posisjon')
plt.plot(tid,deltax,tid,deltax_analytisk)
plt.title('Standardavvik i posisjon',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$\Delta x$ (m)',fontsize=20)
plt.show()

Vi observerer at den numeriske løsningen overlapper med den analytiske løsningen i starten. Når bølgepakken kolliderer med en hard vegg, presses den sammen og får redusert standardavvik i posisjon. Vi observerer dessuten et systematisk avvik fra den analytiske løsningen etter hvert som tiden går. Dette skyldes at dispersjonsrelasjonen ikke er perfekt kvadratisk. For økende energier nærmer $E(p)$ seg en lineær funksjon, og som kjent gir lineær dispersjon ingen forbredning av bølgepakken. Dette kommer meget tydelig fram hvis du går tilbake til starten av programmet og setter for eksempel $E_0$ lik 10 eV.

#ANIMASJONSDELEN
#
#Forbereder figur, akser, og plotteelementet som skal animeres:
rhomax = np.max(np.abs(Psi0rekon)**2)
#fig = plt.figure('Bølgepakkeanimasjon',figsize=(16,8))
fig, ax = plt.subplots()
ax.set_xlim(( 0, N*dx))
ax.set_ylim((0, 1.5*rhomax))
line, = ax.plot([], [], lw=2)

#Initialisering; plotter bakgrunnen:
def init():
    line.set_data([], [])
    return line,
    
#Animasjonsfunksjon; kalles sekvensielt:
def animate(i):
    t = i*tidssteg
    if (i % 100 == 0):
       print(i)
       
    #Lager Psi(x,t) 
    Psi_t = np.zeros(Ntot,dtype=np.complex128)
    for n in range(Ntot):
        Psi_t = Psi_t + c[n]*psi_matrix_complex[:,n]*np.exp(-1j*energy[n]*t/hbar)
    
    rho_t = np.abs(Psi_t)**2
    line.set_data(x, rho_t)
    return line,

#Kaller animatoren.
#frames=antall bilder, dvs maxverdi for i
#interval=varighet av hvert bilde i animasjonen målt i millisekunder
anim = animation.FuncAnimation(fig, animate, init_func=init, repeat=False,
                               frames=Ntid, interval=50, blit=True)
#Viser animasjonen
HTML(anim.to_html5_video())

0
100
200
300
400