Numerical Evaluation of the Propagator in the Position Representation

In [BR95] a numerical ``brute force approach'' to determining the kicked harmonic oscillator dynamics is described. Using the present terminology and scaling, this approach can be formulated as a discretization of the integral form of the quantum map (2.37) in the position representation,

$$\big< x \big\vert \psi_{n+1} \big> \; = \; \int\limits _{-\i... ...style -\frac{i}{\hbar}V_0\cos x' } \Big\vert \psi_n \! \Big> ,$$ (3.21)

expressing the time evolution over one period using the position representation propagator

$$G_{\mbox{\scriptsize free}}\big(x,T;x',0\big) \; = \; \fra... ...{ i \, \frac{(x^2+x'^2)\cos T-2xx'} {2\hbar\sin T} \right\}$$ (3.22)

of the unkicked harmonic oscillator (2.29a) [FH65]. In contrast to the methods discussed in section 3.1 these expressions by construction take into account the FLOQUET nature of the system.

Using the $x$-discretization as given by equations (3.10), the integral expression (3.25) can be approximated by

$$\big< x_j \big\vert \psi_{n+1} \big> \approx \frac{x_{\mbo... ...}{\hbar}V_0\cos x_{j'} } \big< x_{j'} \big\vert \psi_n \big> ,$$ (3.23)

which is once again denoted more concisely in vectorial form,

$$\vec{\psi}_{n+1} \; = \; U^{\mbox{\scriptsize (pr)}} \, \vec{\psi}_n ,$$ (3.24)

with the vectors

$$\vec{\psi}_n \; := \; \left( \begin{array}{c} \big< x_1 \... ...\scriptsize max}}} \big\vert \psi_n \big> \end{array} \right)$$ (3.25)

and the $(j_{\mbox{\scriptsize max}},j_{\mbox{\scriptsize max}})$ propagator matrix $U^{\mbox{\scriptsize (pr)}}$ in the position representation, with the matrix elements

$$U_{jj'}^{\mbox{\scriptsize (pr)}} := \frac{x_{\mbox{\script... ...0\cos x_{j'} }, \; 1\leq j,j'\leq j_{\mbox{\scriptsize max}}.$$ (3.26)

$U^{\mbox{\scriptsize (pr)}}$ takes the same role with respect to the position representation as the propagator $U$ of equation (2.42) does with respect to the eigenrepresentation of the harmonic oscillator. Note that $U^{\mbox{\scriptsize (pr)}}$, which describes the propagation over a full period of the excitation, is the same for all $n$, such that it suffices to determine its matrix elements once and store them for repeated use.

In order to obtain a numerically stable algorithm, $U^{\mbox{\scriptsize (pr)}}$ must be required to be unitary,

$$\sum_{j'=1}^{j_{\mbox{\scriptsize max}}} \, \left\vert U_{jj... ...)}} \right\vert^2 \; \stackrel{!}{=} \; 1 \qquad \forall j,$$ (3.27)

whence with equations (3.26) and (3.30)

$$\frac{(x_{\mbox{\scriptsize max}}-x_{\mbox{\scriptsize min}}... ... \stackrel{!}{=} \;\; 2\pi\hbar \left\vert \sin T \right\vert$$ (3.28)

can be concluded; the parameters $x_{\mbox{\scriptsize min}}$, $x_{\mbox{\scriptsize max}}$, $j_{\mbox{\scriptsize max}}$ of the numerical algorithm, on the left hand side of equation (3.32), have to be chosen as specified by the physical parameters $\hbar$, $T$ on the right hand side. The kick amplitude $V_0$ does not interfere with the unitarity of $U^{\mbox{\scriptsize (pr)}}$.

In addition to the condition (3.32), a boundary condition of the type (3.24a) also needs to be satisfied, i.e. the interval $[x_{\mbox{\scriptsize min}},x_{\mbox{\scriptsize max}}]$ must be chosen large enough to allow the discrete mapping (3.27) to approximate the original (3.25) well enough. This is granted if the wave function $\left< x \left\vert \psi_n \right> \right.$ decays rapidly enough with $x$ reaching out to the boundaries of the interval considered.

These two restrictions combined are the reason why the algorithm given by the mapping (3.27) is of limited practical use when the long-time dynamics of spreading states is to be studied: for delocalized states spreading widely along the $x$-axis, very large values of $j_{\mbox{\scriptsize max}}$ are required by the condition (3.32), in particular when $\hbar$ takes on small values. This means that the memory requirements on the computer grow rapidly, and the evaluation of each iteration of (3.27) quickly becomes more time-consuming. For example, the values $-x_{\mbox{\scriptsize min}}=x_{\mbox{\scriptsize max}}=30$, $\hbar =0.01$ and $T=\pi /2$ lead to $j_{\mbox{\scriptsize max}}\approx 57000$, making it practically impossible to store the matrix elements $U_{jj'}^{\mbox{\scriptsize (pr)}}$ on a typical workstation and slowing down the speed of the computation considerably. It is important to keep in mind this limitation of the algorithm when working with it in practice.

Conversely, if the dynamics is followed for times $nT$ not too large, such that the initial wave packets have not spread too much, then equation (3.27) provides a simple and efficient way to evaluate the quantum map (2.37), especially for larger values of $\hbar$.