The Propagator in the Eigenrepresentation
of the Free Harmonic Oscillator

In this section I discuss some of the aspects of the practical implementation of the iteration scheme outlined in subsection 2.1.3.

From a formal point of view the methods formulating ${\hat{U}}$ in the position representation (section 3.2) and the eigenrepresentation of the harmonic oscillator (subsection 2.1.3 and the present section) obviously are very similar -- for example, both equation (2.43) and equation (3.28) are explicit linear mappings of vectors of complex numbers.

What is more, the dimensions of the matrices $U$ and $U^{\mbox{\scriptsize (pr)}}$ used for the iterations are of the same order of magnitude, too. This can be seen as follows. The dimension $j_{\mbox{\scriptsize max}}$ to be used for $U^{\mbox{\scriptsize (pr)}}$ has been derived in the previous subsection; on the other hand, the dimension needed for $U$ is determined by the maximum index $m_{\mbox{\scriptsize max}}$ of the set $\big\{ \left\vert m \right> \big\vert \, \mbox{$0\leq m\leq m_{\mbox{\scriptsize max}}$} \big\}$ necessary for resolving quantum phase space structures at a distance of $x_{\mbox{\scriptsize max}}$ from the origin. Roughly speaking, $\vert\left< x \left\vert m \right> \right.\!\vert^2$ is peaked around the corresponding classical turning points at $\pm\sqrt{\hbar(2m+1)}$ and decays exponentially for larger values of $\vert x\vert$.^3.2This means that, as a rule of the thumb, the dimension of $U$ needed for describing the dynamics on the interval $[-x_{\mbox{\scriptsize max}},x_{\mbox{\scriptsize max}}]$ is not very much larger than $x_{\mbox{\scriptsize max}}^2/2\hbar$. This estimate is close to the value of $j_{\mbox{\scriptsize max}}$ determined in equation (3.32).

There are some decisive differences between the matrices used for the two methods, as well. First of all, with respect to identical dimensions of the matrices, storing the matrix elements of $U$ only requires a quarter of the memory needed for storing the matrix elements of $U^{\mbox{\scriptsize (pr)}}$: in practice, the kick matrix elements $K_{mm'}$ of equation (2.51) -- rather than the full FLOQUET matrix elements $U_{mm'}$ -- are stored in the computer memory, such that the selection rule $K_{mm'}=0$ for $m-m'$ odd and the symmetry property $K_{mm'}=K_{m'm}$ can be used -- see equation (2.51). Then, each time they are needed during the iteration, the $U_{mm'}$ are trivially calculated from the $K_{mm'}$ using the splitting (2.45). Storing the $K_{mm'}$ requires computer memory for roughly $m_{\mbox{\scriptsize max}}^2/4$ complex numbers, i.e. for the typical value of $m_{\mbox{\scriptsize max}}=6000$ some 138 MBytes are needed. The approximate CPU time per iteration of the quantum map then typically varies from 45 sec (on a Pentium II with 350 MHz) to 10 sec (Pentium 4 with 2 GHz).

Second, in contrast to the evaluation of (3.30), computation of the $K_{mm'}$ and thus of the $U_{mm'}$ is a nontrivial task. Equation (2.51) shows that the formula for the $K_{mm'}$ is composed of several terms (BESSEL functions, generalized LAGUERRE polynomials, factorials, exponentials) the absolute values of which can -- and in practice do -- differ greatly, such that all kinds of numerical cancellations, overflows and underflows are to be expected -- and indeed do occur -- when (2.51) is evaluated directly. Therefore, this calculation has to be implemented in a slightly more sophisticated way, basically by exploiting the recurrence relation

for the generalized LAGUERRE polynomials [AS72]. In this way, not only the abovementioned numerical problems are avoided, but also the computation of the whole matrix is sped up considerably.

Starting from the kick matrix elements in the form of equation (2.48) and using the definition

$$N_{mm'}(l) \; := \; \left< m \left\vert e^{ il{\hat{x}}} \r... ...+ \left< m \left\vert e^{-il{\hat{x}}} \right\vert m' \right>$$ (3.28)

I obtain

$$K_{mm'} \; = \; {\mbox{J}}_0\left( \frac{V_0}{\hbar} \right... ...l \, {\mbox{J}}_l\left(\frac{V_0}{\hbar}\right) N_{mm'}(l) .$$ (3.29)

The matrix elements $N_{mm'}(l)$ in this expression need to be studied with respect to their behaviour along the diagonals given by

$$\delta m \; := \; m-m' \; = \; \mbox{const.}$$ (3.30)

With equation (2.50), $N_{mm'}(l)$ can be rewritten in terms of a suitably normalized variant,

$$\tilde{{\mbox{L}}}_{m'}^{(\delta m)}(z) \; := \; \frac{1} ... ...lta m+m'-1)\cdots(m'+1)}} \; {\mbox{L}}_{m'}^{(\delta m)}(z),$$ (3.31)

of the generalized LAGUERRE polynomials:

$$N_{\delta m+m',m'}(l) \; = \; \displaystyle \left\{ \dis... ... m$\ odd. \rule[-0.3cm]{0.0cm}{0.3cm} } \end{array} \right.$$ (3.32)

In this expression, working with $\tilde{{\mbox{L}}}_{m'}^{(\delta m)}$ (rather than ${\mbox{L}}_{m'}^{(\delta m)}$) is advantageous, as it collects all $m'$-dependent terms. And using the indices $\delta m$, $m'$ rather than $m$, $m'$ is appropriate here, because in the next step the $N_{\delta m+m',m'}(l)$ are calculated by repeated application of a recurrence relation along the diagonals of the matrix, where $m'$ specifies the individual matrix elements on the diagonal determined by $\delta m$. Note that equations (3.37) and (3.38) hold for $\delta m\geq 0$ only, because the same is true for equation (2.50). For $\delta m<0$, by equation (2.51) the symmetry $K_{mm'}=K_{m'm}$ allows to determine the kick matrix elements with equations (3.35-3.38), too.

The $N_{\delta m+m',m'}(l)$ can now be calculated efficiently using the recurrence relation

which is obtained by evaluating the $\tilde{{\mbox{L}}}_{m'}^{(\delta m)}\big( \hbar l^2/2 \big)$ in equation (3.38) by means of the definition (3.37) and the recursion relation (3.33) of the generalized LAGUERRE polynomials.

The numerical problems mentioned above with respect to the direct evaluation of equation (2.51) are reduced to a minimum with this algorithm, because multiplying and adding up very large and very small numbers at the same time is avoided here as far as possible: the definition (3.37) avoids the factorials of equation (2.51), and the respective absolute values of the numerators and denominators in each of the fractions in the recursion (3.39a) are constructed to be as comparable as possible.

Similarly, when in the next chapters $\left< x \left\vert m \right> \right.$ needs to be evaluated numerically for larger values of $m$, this cannot be done using equation (2.39) directly. Rather, a scaled variant of the HERMITE polynomials is defined,

$$\tilde{{\mbox{H}}}_m(z) \; := \; \frac{1}{\sqrt{2^mm!}} \, {\mbox{H}}_m(z),$$ (3.32)

which can be evaluated efficiently and in a numerically sound way even for larger values of $m$ via the recurrence relation

that follows immediately from the conventional recurrence relation

for the HERMITE polynomials [AS72]. The harmonic oscillator eigenstates in the position representation are then obtained as

$$\left< x \left\vert m \right> \right. \; = \; \frac{1}{\sq... ...\tilde{{\mbox{H}}}_m \! \left( \frac{x}{\sqrt{\hbar}} \right).$$ (3.31)

With the $N_{\delta m+m',m'}(l)$ computed as described above, the series in equation (3.35) can be evaluated term by term. In the parameter range considered for the numerical calculations in the following chapters, the series can typically be truncated at values of $l$ not exceeding 200, while keeping the error induced by the truncation reasonably small. (For the maximum relative error due to this truncation, $10^{-14}$ was a somewhat typical value in the practical calculations.)

As in the previous sections, boundary conditions have to be imposed here, too. After each time step it is checked if the condition (3.24a) is still satisfied. In practice, for diffusing states violation of this boundary condition -- due to the implicit cut-off error as a result of the finite number $m_{\mbox{\scriptsize max}}$ of basis elements $\left\vert m \right>$ used for expanding the states -- is one of the two main sources of numerical error. The other is the unavoidable rounding error generated in the course of the numerical computation of the $K_{mm'}$.

Footnotes

....^3.2

Note that the same applies with respect to the $p$-direction. The quantum phase space portion covered by $\big\{ \left\vert m \right> \big\vert \, \mbox{$0\leq m\leq m_{\mbox{\scriptsize max}}$} \big\}$ is a circular disc with the approximate area $\pi\hbar(2m_{\mbox{\scriptsize max}}+1)$, i.e. it grows linearly with $m_{\mbox{\scriptsize max}}$, the growth rate being given by the square of the oscillator length ($\sqrt{\hbar}$ in the scaling used here).
See figures A.1-A.3 in appendix A for a graphical demonstration of this property of the harmonic oscillator eigenstates.