Split-Operator-Methode
# Split-Operator-Methode

Die Split-Operator-Methode (SOP) ist ein numerisches Verfahren mit dem die zeitabhängige Schrödingergleichung gelöst werden kann. Bei der Methode wird der Hamiltonoperator ${\hat {H}}$ in einen kinetischen Teil ${\hat {T}}$ (Impulsteil) und in einen Potentialteil ${\hat {V}}$ gespalten und einzeln angewendet. Dabei wird von der schnellen Fourier-Transformation (FFT) Gebrauch gemacht, um zwischen Impulsraum und Ortsraum zu wechseln.

## Die Schrödingergleichung

Die zeitabhängige Schrödingergleichung ist definiert als

    $i\hbar {\frac {\partial \psi (x,t)}{\partial t}}={\hat {H}}(t)\psi (x,t),$

wobei $\textstyle {\hat {H}}(t)=-{\frac {\hbar ^{2}}{2m}}\Delta +V(t)$ der Hamiltonoperator ist.

Die Wellenfunktion $\psi (x,t)$ wird im Ortsraum auf einem äquidistanten Gitter dargestellt. Als Startwerte werden die Werte von $\psi (x,t_{0})$ zur Zeit $t_{0}$ an den Gitterpunkten vorgegeben. Durch das Verfahren wird die Wellenfunktion zu einem späteren Zeitpunkt $t=t_{0}+\Delta t$ berechnet.

Die Wirkung des Hamiltonoperators $\textstyle {\hat {H}}(t)={\frac {{\hat {\mathbf {p} }}^{2}}{2m}}+{\hat {V}}(t)$ auf eine Wellenfunktion ${\hat {H}}\psi ={\hat {T}}\psi +{\hat {V}}\psi$ wird mit der schnellen Fourier-Transformation berechnet. Dazu wird neben dem Gitter im Ortsraum auch ein Gitter im Impulsraum benötigt. Die Auflösung im Impulsraum $\Delta k={\tfrac {2\pi }{L}}$ ist durch die Länge $L$ des Gitters im Ortsraum festgelegt. Es gilt $\Delta k\Delta x={\tfrac {2\pi }{N}}$, wobei $N$ die Anzahl der Gitterpunkte ist.

## Anwendung der diskreten Fourier-Transformation

Der Potentialoperator ${\hat {V}}$ besitzt im Ortsraum eine diagonale Matrixdarstellung und wirkt daher lokal auf jeden Gitterpunkt $x_{i}$:

    $\left({\hat {V}}(t)\psi \right)(x_{i})=V(x_{i},t)\cdot \psi (x_{i}).$

Genauso wird der kinetische Operator ${\hat {T}}$ mit seiner diagonalen Darstellung im Impulsraum berechnet. Für jeden Gitterpunkt $k_{i}$ gilt:

    $\left({\hat {T}}{\tilde {\psi }}\right)(k_{i}))={\frac {\hbar ^{2}}{2m}}{k_{i}}^{2}{\tilde {\psi }}(k_{i}).$

Dabei ist die diskrete Darstellung der Wellenfunktion ${\tilde {\psi }}(k_{i})$ im Impulsraum durch die diskrete Fourier-Transformation ${\hat {Z}}^{\dagger }$ gegeben:

    ${\tilde {\psi }}(k_{i})=\langle k_{i}|\psi \rangle =\sum _{x_{j}}\langle k_{i}|x_{j}\rangle \langle x_{j}|\psi \rangle ={\frac {\Delta x}{\sqrt {2\pi }}}\sum _{x_{j}}e^{-ik_{i}x_{j}}\psi (x_{j})$

In Vektorschreibweise lautet diese Gleichung

    ${\vec {\tilde {\psi }}}=c^{-1}{\hat {Z}}^{\dagger }{\vec {\psi }}$

mit

    ${\vec {\tilde {\psi }}}:=({\tilde {\psi }}(k_{0}),\dotsm ,{\tilde {\psi }}(k_{N-1}))^{T}$
    ${\vec {\psi }}:=(\psi (x_{0}),\dotsm ,\psi (x_{N-1}))^{T}$
    $Z_{ij}:=N^{-{\frac {1}{2}}}e^{+ik_{i}x_{j}}$
    $c:={\tfrac {\sqrt {2\pi N}}{L}}.$

Entsprechend erhält man für die Rücktransformation in den Ortsraum

    $\psi (x_{j})={\frac {\Delta k}{\sqrt {2\pi }}}\sum _{k_{i}}e^{+ik_{i}x_{j}}{\tilde {\psi }}(k_{i})$

beziehungsweise

    ${\vec {\psi }}=c{\hat {Z}}{\vec {\tilde {\psi }}}$

mit den Gitterschrittweiten $\Delta x={\tfrac {L}{N}}$ bzw. $\Delta k={\tfrac {2\pi }{L}}$. Hierbei ist $L$ die Länge des Gitters im Ortsraum und $N$ die Zahl der Punkte im Orts- und Impulsraum. Die Konstante $c$ wird nur benötigt, wenn die richtige Normierung der Funktion ${\tilde {\psi }}$ gewünscht wird. Die Fourier-Transformation erhält die Norm der Vektoren ${\vec {\tilde {\psi }}}$ und ${\vec {\psi }}$.

## Split-Operator-Methode

Die Berechnung der $e$-Funktion eines Operators wird in der Diagonaldarstellung des Operators besonders einfach. Die Split-Operator-Methode verwendet eine Zerlegung des Hamiltonoperators in die Operatoren für kinetische Energie ${\hat {T}}$ und für potentielle Energie ${\hat {V}}$, welche im Impuls- bzw. Ortsraum Diagonalform annehmen.

Der durch die Nicht-Vertauschbarkeit von ${\hat {T}}$ und ${\hat {V}}$ entstehende Fehler kann durch die symmetrische Aufspaltung

    $e^{-{\frac {i}{\hbar }}{\hat {H}}\Delta t}\approx e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}\cdot e^{-{\frac {i}{\hbar }}{\hat {V}}\Delta t}\cdot e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}$

auf Terme der Größenordnung ${\Delta t}^{3}$ reduziert werden: Mit ${\hat {X}}:=-{\tfrac {i}{\hbar }}{\hat {T}}\Delta t$ und ${\hat {Y}}:=-{\tfrac {i}{\hbar }}{\hat {V}}\Delta t$ erhält man für die rechte Seite

    $\exp \left({\frac {\hat {X}}{2}}\right)\exp \left({\hat {Y}}\right)\exp \left({\frac {\hat {X}}{2}}\right)=\exp \left({\frac {\hat {X}}{2}}+{\hat {Y}}+{\frac {\hat {X}}{2}}+\underbrace {{\frac {1}{2}}\left[{\frac {\hat {X}}{2}},{\hat {Y}}\right]+{\frac {1}{2}}\left[{\hat {Y}},{\frac {\hat {X}}{2}}\right]} _{0}+{\frac {1}{12}}\left[\left[{\frac {\hat {X}}{2}},{\hat {Y}}\right],{\hat {X}}+2{\hat {Y}}\right]+\dotsm \right).$

Der führende Fehlerterm ist somit proportional zu ${\Delta t}^{3}\left[\left[{\hat {T}},{\hat {V}}\right],{\hat {T}}+2{\hat {V}}\right]$.

## Diagonalform

Eine Koordinatentransformation ${\hat {Z}}^{\dagger }$ vom Orts- in den Impulsraum ermöglicht eine einfache Berechnung von

    $e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}\psi ={\hat {Z}}e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}{\hat {Z}}^{\dagger }\psi .$

Mit der diagonalen Darstellung des Operators der kinetischen Energie

    ${\hat {T}}={\hat {Z}}^{\dagger }{\hat {T}}{\hat {Z}}={\begin{pmatrix}{\frac {\hbar ^{2}}{2m}}{k_{0}}^{2}&0&\cdots &0\\0&\ddots &0&0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &{\frac {\hbar ^{2}}{2m}}k_{N-1}^{2}\\\end{pmatrix}}$

erhält man

    $e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}={\begin{pmatrix}\ddots &0&\cdots &0\\0&e^{-{\frac {i}{\hbar }}{\frac {\Delta t}{2}}{\frac {\hbar ^{2}}{2m}}k_{i}^{2}}&0&0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &\ddots \\\end{pmatrix}}.$

Die Koordinatentransformation erfolgt auf dem $N$-Punkt-Gitter $x_{0},\dotsm ,x_{N-1}$ mit Hilfe der diskreten Fourier-Transformation:

    ${\tilde {\psi }}(k_{i})=N^{-{\frac {1}{2}}}\sum _{j=0}^{N-1}\psi (x_{j})e^{-ik_{i}x_{j}}$  für $i=0,\dotsm ,N-1$

oder ${\vec {\tilde {\psi }}}={\hat {Z}}^{\dagger }{\vec {\psi }}$.

## Numerischer Algorithmus

Durch Zusammenfassen der aufeinanderfolgenden Terme ${\hat {Z}}e^{-{\frac {i}{\hbar }}{\frac {\hat {T}}{2}}\Delta t}{\hat {Z}}^{\dagger }$ zweier Zeitschritte lässt sich die Zahl der Fourier-Transformationen, d. h. der numerische Aufwand, reduzieren: ${\hat {Z}}^{\dagger }{\hat {Z}}=1$, und die beiden $e$-Funktionen mit ${\frac {\hat {T}}{2}}$ ergeben $e^{-{\frac {i}{\hbar }}{\hat {T}}\Delta t}$.

Die Wellenfunktion nach $n$ Zeitschritten erhält man also durch:

* Fourier-Transformation von $\psi _{0}$
* Multiplikation mit den Diagonalelementen $e^{-{\frac {i}{\hbar }}{\frac {\Delta t}{2}}{\frac {\hbar ^{2}}{2m}}k_{i}^{2}}$ (halber Zeitschritt)
* Rücktransformation
* Multiplikation mit den Diagonalelementen $e^{-{\frac {i}{\hbar }}V_{i}\Delta t}$
* Fourier-Transformation
* Multiplikation mit den Diagonalelementen $e^{-{\frac {i}{\hbar }}\Delta t{\frac {\hbar ^{2}}{2m}}k_{i}^{2}}$ (ganzer Zeitschritt)
* usw., bis beim letzten Schritt noch einmal eine Multiplikation mit halben Zeitschritt wie in der zweiten Zeile notwendig wird.
