Navn på gruppens medlemmer:

Om de numeriske øvingene i kvantemekanikk

Målet med de numeriske øvingene i kvantemekanikk er

  • å bruke numeriske metoder og visualiseringer for å forstå kvantemekanikk bedre.
  • å utvikle en generell "verktøykasse" du kan bruke også når du arbeider med andre øvinger og lærer nytt stoff i faget.
  • å se ligningene i kvantemekanikken og løsninger av disse fra et generelt perspektiv.
  • å utvikle ferdigheter innen programmering og numerikk.

Vi legger opp til bruk av Python sammen med bibliotekene numpy, scipy og matplotlib for numerikk, vitenskapelige beregninger og plotting. Bruk disse for alt de er verdt! God bruk av disse bibliotekene vil la deg uttrykke deg mer konsist enn om du kun bruker standard Python-funksjonalitet. Vi forventer at du selv finner fram til relevant funksjonalitet og dokumentasjon i bibliotekene. De er alle godt dokumentert på Internett.

I Jupyter Notebook og Jupyter Lab er det mulig å bruke ulike backends for plotting med matplotlib. Dette er grovt forklart forskjellige underliggende "motorer" som bestemmer utseende og funksjonaliteten til figurene som produseres. Avhengig av hvilket av programmene du bruker, vil du erfare at en backend fungerer bedre enn andre. I Jupyter Notebook fungerer notebook-backenden best og uten behov for installasjon av ekstra programvare, men i Jupyter Lab må tillegget jupyter-matplotlib installeres for å få den optimale widget-backenden til å fungere. Begge programmene støtter også inline-backenden uten behov for tilleggsprogramvare, men denne produserer mindre fleksible figurer og bør kun brukes som reserveløsning. Prøv deg selv fram med backendene foreslått under for å finne den som fungerer best for deg.

In [1]:
# uncomment ONE line to choose matplotlib backend
# if using Jupyter Notebook, use interactive "notebook" backend for best results
# if using Jupyter Lab, use interactive "widget" backend for best results
# if both fail, use static "inline" backend
%matplotlib notebook 
#%matplotlib widget 
#%matplotlib inline 

import numpy as np
import matplotlib.pyplot as plt

Kort om fysiske enheter i numeriske beregninger

I numeriske beregninger må vi ta hensyn til den begrensede presisjonen og størrelsen til flyttallene som datamaskinen bruker for å representere relle tall. I kvantemekanikken møter vi spesielt ofte på Plancks reduserte konstant $\hbar \approx 6.63 \cdot 10^{-34} \text{ Js}$ og gjerne kvadratet $\hbar^2$. Det er i utgangspunktet ingenting i veien for å benytte SI-enheter for alle størrelser som opptrer i disse øvingene. Et alternativ er å benytte atomære enheter. Dette enhetssystemet er skreddersydd for beregninger på atomært nivå. Her er blant annet $\hbar$, elektronmassen $m_e$ og elementærladningen $e$ definert til å ha tallverdi $1$. For eksempel uttrykkes energier som multipler av én hartree, $E_h = \hbar^2 / m_e a_0^2 \approx 4.36 \cdot 10^{-18} \text{ J}$. En annen løsning er å benytte enheter som $\text{nm}$ og $\text{eV}$ for lengder og energier. En hartree tilsvarer ca. 27.2 eV, dvs. det dobbelte av grunntilstandsenergien i hydrogenatomet.

Velg selv hensiktsmessige enheter til bruk i beregningene, men vær oppmerksom på den begrensede presisjonen til flyttall!

Vi bruker her atomære enheter.

In [2]:
hbar = 1

Numerisk løsning av den tidsuavhengige Schrödingerligningen

Å løse den tidsuavhengige Schrödingerligningen $$ \hat{H} \psi = -\frac{\hbar^2}{2 m} \psi'' + V \psi = E \psi, $$ dvs. å bestemme energiegenverdier $E$ og tilhørende energiegenfunksjoner $\psi(x)$ for et gitt potensial $V(x)$, er et sentralt problem i kvantemekanikken. Dette er ofte ingen enkel oppgave. Selv for potensialer som gir ligningen analytiske løsninger, kreves det ofte betydelig innsats og bruk av spesielle teknikker for å komme fram til disse. Vi skal her se på en elegant og generell teknikk for å løse ligningen numerisk for et vilkårlig potensial (i én dimensjon).

Numeriske løsningsmetoder innebærer alltid en viss avgrensning og diskretisering for å gjøre problemet endelig og håndterlig for en datamaskin. Vi avgrenser her delen av rommet vi ser på til å ligge mellom to endepunkter $x_0$ og $x_{N+1}$ og deler opp intervallet mellom dem i punktene $x_0, x_1, \ldots, x_N, x_{N+1}$ med lik avstand $\Delta x$ mellom hvert punkt. Utenfor dette området definerer vi potensialet til $V(x \leq x_0) = V(x \geq x_{N+1}) = \infty$, slik at $\psi(x \leq x_0) = \psi(x \geq x_{N+1}) = 0$ og det kun er bølgefunksjonens verdier på rutenettet $\boldsymbol{x} = [x_1, \ldots, x_N]^T$ som er ukjente og av interesse. Til hvert punkt tilordner vi verdiene $\psi_i = \psi(x_i)$ og $V_i = V(x_i)$ til energiegenfunksjonene og potensialet, og vi refererer til verdiene av funksjonene i alle punktene ved hjelp av vektorene $\boldsymbol{V} = [V_1, \ldots, V_N]^T$ og $\boldsymbol{\psi} = [\psi_1, \ldots, \psi_N]^T$.

En intuitiv og enkel tilnærming av den deriverte til en funksjon er den sentrale differansen $$ \psi'(x) = \frac{\psi(x + \Delta x / 2) - \psi(x - \Delta x / 2)}{\Delta x} $$ Om vi bruker denne tilnærmingen to ganger, kan vi også tilnærme den andrederiverte som $$ \psi''(x) = \frac{\psi'(x + \Delta x / 2) - \psi'(x - \Delta x / 2)}{\Delta x} = \frac{\psi(x + \Delta x) - 2 \psi(x) + \psi(x - \Delta x)}{\Delta x^2}$$

Ved å sette denne tilnærmingen inn i den tidsuavhengige Schrödingerligningen, kan vi tilnærme den numerisk som $$-\frac{\hbar^2}{2 m} \frac{\psi_{i+1} - 2 \psi_i + \psi_{i-1}}{{\Delta x}^2} + V_i \psi_i = E \psi_i \qquad \text{for}\,\, i = 1, 2, \ldots, N$$

Vi kan uttrykke denne ligningen på en elegant måte ved å innføre $N \times N$-Hamiltonmatrisen $H$ med elementer $$ H_{i j} = \begin{cases} \hbar^2 / (m \Delta x^2) + V_i & \text{for} \,\, i = j & \text{(på diagonalen)} \\ -\hbar^2 / (2 m \Delta x^2) & \text{for} \,\, i = j \pm 1 & \text{(på semidiagonalene)} \\ 0 & \text{ellers} \\ \end{cases} $$ og benytte oss av vektoren $\boldsymbol{\psi} = [\psi_1, \ldots, \psi_N]^T$. Den tar da formen $$ H \boldsymbol{\psi} = E \boldsymbol{\psi} $$ Energiene $E$ og energiegenfunksjonene $\boldsymbol{\psi}$ er dermed egenverdier og egenvektorer til matrisen $H$!

Numerikkbiblioteker har funksjonalitet for å finne egenverdier og egenvektorer til vilkårlige matriser. De har gjerne også spesialiserte funksjoner som gjør dette mer effektivt for matriser med en spesiell form, for eksempel som den tridiagonale (samt reelle og symmetriske) formen til matrisen $H$.

Skriv en funksjon som beregner og returnerer alle energiegenverdier $E$ og tilhørende energiegenfunksjoner $\boldsymbol{\psi}$ for en partikkel med masse $m$ som befinner seg i et gitt potensial $\boldsymbol{V}$ på rutenettet $\boldsymbol{x}$. Normér energiegenfunksjonene i forstanden $\int |\psi|^2 \mathrm{d}x$ = 1.

Vi oppretter matrisen $H$ som beskrevet over og benytter numpy.linalg.eigh() til å finne egenverdiene og egenvektorene til $H$. For å senere kunne iterere enkelt gjennom egenvektorene (for wave in waves: ...), transponerer vi også matrisen med egenvektorene. Vi normerer så alle egenvektorene ved å utnytte at $\psi / \sqrt{\int |\psi|^2 \mathrm{d} x}$ alltid er normert. Vi tilnærmer normeringsintegralene som $\int |\psi|^2 \mathrm{d}x \approx \sum_i |\psi_i|^2 \Delta x$.

In [3]:
def hamiltonian(x, v, m):
    dx = x[1] - x[0]
    diag = hbar**2 / (m*dx**2) + v
    semidiag = np.full(len(x)-1, -hbar**2 / (2*m*dx**2))
    H = np.diag(diag, k=0) + np.diag(semidiag, k=+1) + np.diag(semidiag, k=-1)
    return H

def normalize(x, wave):
    dx = x[1] - x[0]
    norm = np.sum(np.abs(wave)**2) * dx
    wave /= np.sqrt(norm)
    return wave

def get_stationary_states(x, v, m):
    H = hamiltonian(x, v, m)
    eigvals, eigvecs = np.linalg.eigh(H)
    energies, waves = eigvals, eigvecs.T
    
    for i in range(0, len(waves)):
        waves[i] = normalize(x, waves[i])
    
    return energies, waves

I resten av øvingen skal vi rett og slett bare bruke denne numeriske løsningsmetoden på en rekke forskjellige potensialer. I noen av eksemplene skal vi også sammenligne de numeriske verdiene med analytiske resultater. For å gjøre denne prosessen så enkel som mulig, foreslår vi at du her skriver én "ultimat" plottefunksjon som du kan gjenbruke i alle disse oppgavene.

Skriv en funksjon som framstiller potensialet $\boldsymbol{V}$, energiegenverdier $E$ og energiegenfunksjoner $\boldsymbol{\psi}$ (eller absoluttkvadratene $|\boldsymbol{\psi}|^2$) på rutenettet $\boldsymbol{x}$ grafisk. Funksjonen skal også kunne brukes til å sammenligne to sett med (numeriske og analytiske) energier og energiegenfunksjoner.

Gjør gjerne dette parallelt med resten av oppgavene, slik at du kan tilpasse framstillingen basert på behovene som oppstår. Se gjerne i forelesningsnotater, bøker og søk rundt på Internett for å få litt inspirasjon til hvordan framstillingen kan gjøres.

Denne oppgaven er åpen for kreativitet og kan gjøres på utallige måter.

Vi går for én "superfigur" hvor vi plotter potensialet sammen med de stasjonære tilstandene rundt deres korresponderende energier. Et eventuelt annet sett med energier og energiegenfunksjoner plotter vi enkelt og greit oppå det første, da vi forventer at de skal overlappe. Vi åpner også for å skalere bølgefunksjonene etter behov, for å unngå at hele figuren blir et stort kluss av bølgefunksjoner på kryss og tvers.

In [4]:
def plot_stationary_states(x, v, energies, waves, energies2=None, waves2=None, wave_scale=1):
    fig, ax = plt.subplots()

    ax.set_xlabel("$x / a_0$")
    ax.set_ylabel("$E / E_h$, $\psi$ or $|\psi|^2$ (qualitative)")
    
    # determine appropriate area of plot to view
    ymin = min(energies[0] + np.min(waves[0]), np.min(v))
    ymax = energies[-1] + np.max(waves[-1])
    dy = ymax - ymin
    ymin -= 0.05 * dy
    ymax += 0.05 * dy
    ax.set_ylim(ymin, ymax)
    
    # plot potential extended with (seemingly) infinitely high walls
    dx = x[1] - x[0]
    x_ext = np.concatenate(([x[0]-dx], x, [x[-1]+dx]))
    v_ext = np.concatenate(([ymax], v, [ymax]))
    ax.plot(x_ext, v_ext, color="black", linewidth=5, label="V")

    for i in range(0, len(energies)):
        energy = energies[i]
        wave = waves[i] * wave_scale
        
        color_index = i % 10 # there are only 10 colors
        color = f"C{color_index}"
        ax.axhline(energy, color=color, linestyle="dashed") # plot line for energy
        ax.plot(x, energy + wave, color=color) # plot wavefunction about energy line
        
        # optionally plot second set of energies and waves for comparison
        if energies2 is not None and waves2 is not None:
            energy2 = energies2[i]
            wave2 = waves2[i] * wave_scale
            ax.plot(x, energy2 + wave2, color=color)
        
    ax.legend()

Partikkel i boks

Et av de første kvantemekaniske problemene vi støter på er partikkel i boks. Her er potensialet, de normerte energiegenfunksjonene og energiegenverdiene $$ V(x) = \begin{cases}0 & \text{for}\,\, 0 \leq x \leq L \\ \infty & \text{ellers} \end{cases}, \quad \psi(x) = \sqrt{\frac{2}{L}} \sin{\frac{n \pi x}{L}}, \quad E = \frac{n^2 \pi^2 \hbar^2}{2 m L^2}, \quad \quad n = 1, 2, \ldots, $$

Sammenlign numeriske og analytiske verdier for noen energier og energiegenfunksjoner for et elektron i en boks grafisk.

Hvordan er spredningen i energinivåene?

Vi deler opp et intervall langs $x$-aksen i mange punkter og setter potensialet til $V = 0$ i alle disse. Vi gir så disse posisjonene og potensialverdiene til funksjonen vi laget tidligere for å finne numeriske verdier for energiene og energiegenfunksjonene. Vi definerer et par funksjoner for å beregne analytiske energier og energiegenfunksjoner fra uttrykkene over. Til slutt bruker vi plottefunksjonen vi laget tidligere til å sammenligne resultatene.

Fra både det analytiske uttrykket for energien og figuren ser vi at energiene øker kvadratisk (dvs. med $n^2$), slik at avstanden mellom energiene øker med energien.

In [5]:
def infbox_wave(n, x):
    L = x[-1] - x[0]
    return np.sqrt(2/L) * np.sin(n*np.pi*x/L)

def infbox_energy(n, x):
    L = x[-1] - x[0]
    return n**2 * np.pi**2 * hbar**2 / (2 * m * L**2)

m = 1
N = 10 # compare N states with lowest energy
x = np.linspace(0, 10, 1000)
v = 0 * x

energies, waves = get_stationary_states(x, v, m)
energies, waves = energies[:N], waves[:N]

# compare with analytic results
energies2 = [infbox_energy(n, x) for n in range(1, N+1)]
waves2 = [infbox_wave(n, x) for n in range(1, N+1)]

plot_stationary_states(x, v, energies, np.abs(waves)**2, energies2, np.abs(waves2)**2, wave_scale=0.5)

Harmonisk oscillator

Et annet standard kvantemekanisk problem er den harmoniske oscillatoren med $$ V(x) = \frac{1}{2}m \omega^2 x^2, \quad \psi(x) = \frac{1}{\sqrt{2^n n!}} \cdot \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \cdot \exp{\left(-\frac{m \omega x^2}{2 \hbar}\right)} \cdot H_n\left(\sqrt{\frac{m \omega}{\hbar}}x\right), \quad E = \left(n+\frac{1}{2}\right)\hbar \omega, \quad \quad n = 0, 1, 2, \ldots $$ Funksjonene $H_n(y)$ (med dimensjonsløs $y$) kalles (fysikerens) Hermitepolynomer. De er tilgjengelige i numerikkbiblioteker, men kan også beregnes fra rekursjonsrelasjonen $$ H_n(x) = 2 x H_{n-1}(x) - 2 (n-1) H_{n-2}(x), \quad H_0(x) = 1, \quad H_1(x) = 2 x$$

Den harmoniske oscillatoren er spesielt interessant i topartikkelsystemer, der et problem med for eksempel to atomer med masse $m_1$ og $m_2$ i et toatomig molekyl reduseres til et ekvivalent enpartikkelproblem med redusert masse $m = m_1 m_2 / (m_1 + m_2)$. Sammen med konstanten $\omega$ utgjør denne et mål på en fjærkonstant som beskriver vibrasjonsbevegelsen mellom de to atomene.

Sammenlign numeriske og analytiske verdier for noen energier og energiegenfunksjoner i en harmonisk oscillator grafisk.

Hvordan er spredningen i energinivåene?

Vi bruker den oppgitte rekursjonsrelasjonen for å beregne Hermitepolynomene, men kunne like gjerne gjort dette med scipy.special.hermite().

Fra både det analytiske uttrykket for energien og figuren ser vi at avstanden mellom energinivåene er konstant.

In [6]:
def hermite(n, x):
    if n == 0:
        return np.ones(len(x))
    if n == 1:
        return 2 * x
    return 2 * x * hermite(n-1, x) - 2 * (n-1) * hermite(n-2, x)

def harmosc_wave(n, x, w):
    return 1/np.sqrt(2**n * np.math.factorial(n)) * (m*w/(np.pi*hbar))**(1/4) * np.exp(-m*w**2*x**2/(2*hbar)) * hermite(n, np.sqrt(m*w/hbar)*x)

def harmosc_energy(n, w):
    return hbar*w*(n+1/2)

w = 1
m = 1
x = np.linspace(-10, +10, 1000)
v = 1/2*m*w**2*x**2

N = 15
energies, waves = get_stationary_states(x, v, m)
energies, waves = energies[:N], waves[:N]

# compare with analytic results
energies2 = [harmosc_energy(n, w) for n in range(0, N)]
waves2 = [harmosc_wave(n, x, w) for n in range(0, N)]

plot_stationary_states(x, v, energies, np.abs(waves)**2, energies2, np.abs(waves2)**2, wave_scale=1.5)

Partikkel i endelige brønnpotensialer

Et tredje velkjent eksempel er enkeltbrønnen $$V(x) = \begin{cases} -V_0 & \text{for}\,\, 0 < x < w \\ 0 & \text{ellers} \end{cases}$$ med bredde $w$ og brønndybde $V_0 > 0$. I dette potensialet finnes ingen analytiske løsninger for de bundne stasjonære tilstandenes energiegenverdier.

Enkeltbrønnen kan generaliseres til et potensial bestående av $N_w$ slike enkeltbrønner plassert ved siden av hverandre med en fast avstand $g$ mellom hver brønn. Med første brønn i $x = 0$, kan vi uttrykke det sammensatte brønnpotensialet stykkevis som

$$V(x) = \begin{cases} 0 & \text{for}\,\, x < 0 \,\, \text{og} \,\, x > N_w (w + g) & \text{(utenfor brønnområdet)} \\ -V_0 & \text{for}\,\, \frac{x}{w+g} - \left\lfloor \frac{x}{w+g} \right\rfloor < \frac{w}{w+g} & \text{(i brønnene)} \\ 0 & \text{ellers} & \text{(mellom brønnene)} \\ \end{cases}$$

Dette er en enkel modell for det periodiske potensialet som et elektron opplever i et fast stoff med en regulær krystallinsk struktur, for eksempel et metall.

Framstill de bundne tilstandene for et elektron både i en enkeltbrønn og i et sammensatt potensial bestående av mange brønner grafisk. Legg inn et passe stort område med $V = 0$ på begge sider av brønnområdet.

Hvordan distribueres energinivåene i potensialet bestående av mange brønner sammenlignet med enkeltbrønnen? Kan du ut fra dette forklare hva vi mener med energibåndstrukturen til et fast stoff ved hjelp av begrepene båndbredde og båndgap?

I enkeltbrønnen fordeler energinivåene seg nokså likt som for en partikkel i boks. Med 5 brønner ved siden av hverandre, derimot, samler energinivåene seg i "kvasi-kontinuerlige" energibånd med forbudte gap imellom. Dette utgjør energibåndstrukturen til det faste stoffet.

In [7]:
def wells(Nw, depth, width, gap, pad, points):
    x1 = -pad
    x2 = Nw * (width + gap) - gap + pad
    x = np.linspace(x1, x2, points)
    v = np.empty(len(x))
    for i in range(0, len(x)):
        if x[i] < 0 or x[i] > Nw * (width + gap):
            v[i] = 0 # outside well area
        elif x[i]/(width+gap) - np.floor(x[i]/(width+gap)) < width/(width+gap):
            v[i] = -depth # in well
        else:
            v[i] = 0 # between wells
    return x, v

m = 1
for Nw in (1, 5):
    x, v = wells(Nw, 3, 5, 1, 20, 1000)

    N = 40
    energies, waves = get_stationary_states(x, v, m)
    energies, waves = energies[:N], waves[:N]

    plot_stationary_states(x, v, energies, waves, wave_scale=0.2)

Partikkel i Deltabarriere eller Deltabrønn

Et fjerde velkjent eksempel er Deltabrønnen $$ V(x) = -\alpha \delta(x) $$ der $\delta(x)$ er Diracs deltafunksjon og $\alpha > 0$ er en konstant. Deltapotensialet er, grovt sagt, ikke annet enn et endelig brønnpotensial der $w \rightarrow 0$ og $V_0 \rightarrow \infty$, men slik at produktet $\alpha = V_0 w$ er endelig.

Framstill noen energier og energiegenfunksjoner for et elektron i en Deltabrønn grafisk. Approksimér Deltabrønnen som en smal og dyp enkeltbrønn.

Hvor mange bundne tilstander ser det ut til å være i Deltabrønnen?

Vi prøver oss fram med en gradvis smalere og dypere enkeltbrønn. Det ser ut til at det er nøyaktig én bunden tilstand i en Deltabrønn.

In [8]:
m = 1
x, v = wells(1, 50, 0.1, 0, 10, 1000)

N = 10
energies, waves = get_stationary_states(x, v, m)
energies, waves = energies[:N], waves[:N]

plot_stationary_states(x, v, energies, waves)

Avsluttende ord

Implementeringen av den numeriske metoden har nok vært krevende, men vi håper du setter pris på sluttresultatet. Selv om en hvilken som helst numerisk teknikk nødvendigvis fører med seg numeriske feil og har visse begrensninger, er det på den andre siden forhåpentligvis tilfredsstillende å kunne bruke én generell metode til å utforske så mange problemer.

Bruk gjerne den numeriske metoden som et verktøy i framtidige situasjoner der det kan være nyttig å få et raskt overblikk over stasjonære tilstander eller energinivåer!

In [ ]: