Obligatorisk Ã¸ving i TFY4215 april 2018Â¶

LÃ¸sningsforslag til oppgave 2Â¶

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

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

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

Velger steglengden $dx = 1$ Ã…, samt bÃ¸lgepakkens senterenergi $E_0 = 1.0$ eV:

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$ , med en fri partikkel antagelse:

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

Som i oppgave 1 velger vi et omrÃ¥de med bredde $L = N_{\rm tot} \cdot dx$ mellom to harde vegger, med $N_{\rm tot} = 1000$ . For Ã¥ unngÃ¥ at bÃ¸lgepakken kommer i berÃ¸ring med de harde veggene, kan vi for eksempel velge et harmonisk potensial omkring midten av "totalomrÃ¥det", og "kontakter" med konstant potensial $V_0 \gg E_0$ pÃ¥ hver side. Og som i oppgave 1 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):

#Ntot = N*10 = totalt antall punkter langs x-aksen
N = 100
#Lager harmonisk potensial med kontakter med konstant potensial pÃ¥ hver side.
#Velger "kontakter" med 4*N punkter pÃ¥ hver side og V = V0 >> E0.
#Mellom de to "kontaktene" er potensialet harmonisk, symmetrisk og med
#minimumsverdi V = 0 i midten.
V0 = 50.0*1.6E-19
V = [V0]*4*N + [V0*((n-1.0*N)/(N*1.0))**2 for n in range(2*N+1)] + [V0]*4*N
#GjÃ¸r om listen til numpy array
V = np.asarray(V)
#Potensialet i enhet eV (i tilfelle plottebehov)
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 = totalt antall punkter langs x-aksen
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$ :

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()

Vi ser at for smÃ¥ verdier av $n$ har vi en lineÃ¦r sammenheng mellom $E_n$ og $n$ . Dette er som forventet for en harmonisk oscillator. For $E > V_0$ fÃ¥r vi en dispersjonsrelasjon lik den vi hadde i oppgave 1.

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

σ = 10 d x

$\sigma = 10 dx$ og

x_{0} = N t o t * d x / 2

$x_0 = Ntot*dx/2$ , som er omtrent midt i det harmoniske potensialet:

#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*10.0
#x0 = pakkens senterposisjon; start minst 5*sigma til hÃ¸yre for x=0
x0 = Ntot*dx/2
#
#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) = {âˆ‘}_{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} = {âˆ«}_{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 mÃ¥ nÃ¥ velge et passende tidssteg og et passende tidsintervall, beregne $\Psi(x,t)$ , og deretter tidsutviklingen av bÃ¸lgepakkens standardavvik i bÃ¥de posisjon og impuls, dvs $\Delta x(t)$ og $\Delta p(t)$ . Med senterenergi $E_0$ har vi de to klassiske vendepunktene for $x$ -verdier som oppfyller $V(x) = E_0$ . Potensialet er pÃ¥ formen $V(x) = \frac{1}{2} m \omega^2 (x - x_0)^2$ , med $x_0 = 5 N dx$ , slik at vendepunktene er ved $x_0 \pm \sqrt{2E_0/m\omega^2}$ . Med $V(x_0) = 0$ og $V(x_0 \pm N dx) = V_0$ blir $\omega = \sqrt{2V_0/m}/N dx$ , slik at svingeperioden i det harmoniske potensialet er $T = 2\pi/\omega = 2\pi N dx \sqrt{m/2V_0}$ . Med tallverdiene $N=100$ , $V_0=50$ eV og $dx = 1$ Ã… er $T = 15$ fs.

#Potensialet er V(x) = (1/2)m*w^2*(x-x0)^2. OmrÃ¥det for det harmoniske potensialet
#har utstrekning 2*N*dx. Perioden blir da T = 2*pi/w = 2*pi*N*dx*sqrt(m/2*V0).
#Med tallverdiene N=100, dx=1Ã…, V0=50eV og m=elektronmassen finner vi
#T = 1.5E-14 s = 15 fs.
#Da er det rimelig Ã¥ velge tidssteget med basis i perioden T, f eks som
#tidssteg = T/100.
#Med M=200 tidsverdier studerer vi dermed ca to hele svingeperioder, dvs fire
#kollisjoner med de myke veggene
#
T = 2*np.pi*N*dx*np.sqrt(m/(2*V0))
tidssteg = T/100
M = 200
xmean = np.zeros(M)
x2mean = np.zeros(M)
xmean2 = np.zeros(M)
pmean = np.zeros(M)
pmean2 = np.zeros(M)
p2mean = np.zeros(M)
norm = np.zeros(M)
halfhbar = np.zeros(M)
tid = np.zeros(M)
null = 0.0 + 0.0j
for i in range(M):
    halfhbar[i] = hbar/2
    tid[i] = i*tidssteg
    #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]*tid[i]/hbar)
    
    Psi_t_conj = np.conjugate(Psi_t)
    d_Psi_t = np.append(np.diff(Psi_t)/dx, null)
    d2_Psi_t = np.append(np.diff(d_Psi_t)/dx, null) 
    rho_t = np.abs(Psi_t)**2
    norm[i] = 0.0
    for n in range(Ntot):
        norm[i] = norm[i] + rho_t[n]

    norm[i] = norm[i]*dx
    xmean[i] = dx*np.dot(x,rho_t)
    x2mean[i] = dx*np.dot(x**2,rho_t)
    xmean2[i] = xmean[i]**2
    pmean[i] = dx*hbar*np.imag(np.dot(Psi_t_conj,d_Psi_t))
    p2mean[i] = dx*hbar*hbar*np.dot(np.abs(d_Psi_t),np.abs(d_Psi_t))
    pmean2[i] = pmean[i]**2

deltax = np.sqrt(x2mean-xmean2)
deltap = np.sqrt(p2mean-pmean2)

Vi plotter utregnede stÃ¸rrelser som funksjon av tiden:

Normeringsintegralet
Midlere posisjon
Midlere impuls
$\langle p^2 \rangle$ og $\langle p \rangle^2$
$\Delta x(t)$
$\Delta p(t)$
$\Delta x \cdot \Delta p$

plt.figure('Normering')
plt.plot(tid,norm)
plt.title('Normering',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('Normering',fontsize=20)
plt.show()


plt.figure('Midlere posisjon')
plt.plot(tid,xmean)
plt.title('Midlere posisjon',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$x$ (m)',fontsize=20)
plt.show()

plt.figure('Midlere impuls')
plt.plot(tid,pmean)
plt.title('Midlere impuls',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$p$ (kg m/s)',fontsize=20)
plt.show()

plt.figure('Midlere kvadratiske impuls og middelimpuls i kvadrat')
plt.plot(tid,p2mean,tid,pmean2)
plt.title('Midlere kvadratiske impuls og middelimpuls i kvadrat',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$p^2$',fontsize=20)
plt.show()


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

plt.figure('Standardavvik i impuls')
plt.plot(tid,deltap)
plt.title('Standardavvik i impuls',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$\Delta p$ (kg m/s)',fontsize=20)
plt.show()

plt.figure('Uskarphetsprodukt')
plt.plot(tid,deltax*deltap,tid,halfhbar)
plt.title('Uskarphetsprodukt',fontsize=20)
plt.xlabel('$t$ (s)',fontsize=20)
plt.ylabel('$\Delta x \cdot \Delta p$',fontsize=20)
plt.show()

Noen kommentarer:

BÃ¸lgepakken holder seg normert.
Midlere posisjon og midlere impuls svinger harmonisk, som forventet.
Standardavvik i posisjon oscillerer og nÃ¥r sine minimumsverdier der bÃ¸lgepakken snur.
Standardavvik i impuls oscillerer ogsÃ¥; her omtrent i motfase med standardavviket i posisjon.
Uskarphetsproduktet blir ikke mindre enn $\hbar/2$ .