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 ${\displaystyle {\hat {H}}}$ in einen kinetischen Teil ${\displaystyle {\hat {T}}}$ (Impulsteil) und in einen Potentialteil ${\displaystyle {\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

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

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

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

Die Wirkung des Hamiltonoperators ${\displaystyle \textstyle {\hat {H}}(t)={\frac {{\hat {\mathbf {p} }}^{2}}{2m}}+{\hat {V}}(t)}$ auf eine Wellenfunktion ${\displaystyle {\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 ${\displaystyle \Delta k={\tfrac {2\pi }{L}}}$ ist durch die Länge ${\displaystyle L}$ des Gitters im Ortsraum festgelegt. Es gilt ${\displaystyle \Delta k\Delta x={\tfrac {2\pi }{N}}}$, wobei ${\displaystyle N}$ die Anzahl der Gitterpunkte ist.

Anwendung der diskreten Fourier-Transformation

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

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

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

${\displaystyle \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 ${\displaystyle {\tilde {\psi }}(k_{i})}$ im Impulsraum durch die diskrete Fourier-Transformation ${\displaystyle {\hat {Z}}^{\dagger }}$ gegeben:

${\displaystyle {\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

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

mit

${\displaystyle {\vec {\tilde {\psi }}}:=({\tilde {\psi }}(k_{0}),\dotsm ,{\tilde {\psi }}(k_{N-1}))^{T}}$

${\displaystyle {\vec {\psi }}:=(\psi (x_{0}),\dotsm ,\psi (x_{N-1}))^{T}}$

${\displaystyle Z_{ij}:=N^{-{\frac {1}{2}}}e^{+ik_{i}x_{j}}}$

${\displaystyle c:={\tfrac {\sqrt {2\pi N}}{L}}.}$

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

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

beziehungsweise

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

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

Split-Operator-Methode

Die Berechnung der ${\displaystyle 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 ${\displaystyle {\hat {T}}}$ und für potentielle Energie ${\displaystyle {\hat {V}}}$, welche im Impuls- bzw. Ortsraum Diagonalform annehmen.

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

${\displaystyle 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 ${\displaystyle {\Delta t}^{3}}$ reduziert werden: Mit ${\displaystyle {\hat {X}}:=-{\tfrac {i}{\hbar }}{\hat {T}}\Delta t}$ und ${\displaystyle {\hat {Y}}:=-{\tfrac {i}{\hbar }}{\hat {V}}\Delta t}$ erhält man für die rechte Seite

${\displaystyle \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 ${\displaystyle {\Delta t}^{3}\left[\left[{\hat {T}},{\hat {V}}\right],{\hat {T}}+2{\hat {V}}\right]}$.

Diagonalform

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

${\displaystyle 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

${\displaystyle {\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

${\displaystyle 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 ${\displaystyle N}$-Punkt-Gitter ${\displaystyle x_{0},\dotsm ,x_{N-1}}$ mit Hilfe der diskreten Fourier-Transformation:

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

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

Numerischer Algorithmus

Durch Zusammenfassen der aufeinanderfolgenden Terme ${\displaystyle {\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: ${\displaystyle {\hat {Z}}^{\dagger }{\hat {Z}}=1}$, und die beiden ${\displaystyle e}$-Funktionen mit ${\displaystyle {\frac {\hat {T}}{2}}}$ ergeben ${\displaystyle e^{-{\frac {i}{\hbar }}{\hat {T}}\Delta t}}$.

Die Wellenfunktion nach ${\displaystyle n}$ Zeitschritten erhält man also durch:

• Fourier-Transformation von ${\displaystyle \psi _{0}}$
• Multiplikation mit den Diagonalelementen ${\displaystyle e^{-{\frac {i}{\hbar }}{\frac {\Delta t}{2}}{\frac {\hbar ^{2}}{2m}}k_{i}^{2}}}$ (halber Zeitschritt)
• Rücktransformation
• Multiplikation mit den Diagonalelementen ${\displaystyle e^{-{\frac {i}{\hbar }}V_{i}\Delta t}}$
• Fourier-Transformation
• Multiplikation mit den Diagonalelementen ${\displaystyle 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.

Literatur

• I. N. Bronstein, K. A. Semendjajew, G. Musiol, H. Muehlig: Taschenbuch der Mathematik. Deutsch Harri GmbH, 2008.
• T. Fließbach: Quantenmechanik: Lehrbuch zur Theoretischen Physik III. 5. Auflage, Spektrum Akademischer Verlag, 2008, ISBN 978-3-8274-2020-6.
• Herbert Sager: Fourier-Transformation. vdf Hochschulverlag, Zürich 2012, ISBN 978-3-7281-3393-9.
• A. Askar, A. S. Cakmak: Explicit integration method for the time‐dependent Schrodinger equation for collision problems. In: Journal of Chemical Physics. Band 68, Nr. 6, 1978, S. 2794–2798, doi:10.1063/1.436072. 
• J. B. Delos: Theory of Electronic Transitions in Slow Atomic Collisions. In: Physical Review. Band 176, Nr. 1, 1968, S. 141–150, doi:10.1103/PhysRev.176.141. 
• Juha Javanainen, Janne Ruostekoski: Symbolic calculation in development of algorithms: split-step methods for the Gross–Pitaevskii equation. In: Journal of Physics A. Band 39, 2006, S. L179–L184, doi:10.1088/0305-4470/39/12/L0. 
• Michael Hintenender: Propagation von Wellenpaketen. In: MPQ-Berichte. MPQ163. Garching 1992 (online).