\documentstyle[epsf]{report} \setlength{\textwidth}{6.0in} \setlength{\textheight}{9.0in} \setlength{\oddsidemargin}{0.5in} \setlength{\evensidemargin}{0.5in} \setlength{\topmargin}{-0.5in} \makeindex % Re & Im for math \def\Re {\mbox{Re}\,} \def\Im{\mbox{Im}\,} \begin{document} \chapter{The Physics of Kaonic Hydrogen} \section{Introduction} Kaonic Hydrogen is an exotic, quasi-stable bound system in which an antikaon, $K^-$, orbits a proton nucleus. The physics of this `atom' proves very interesting to physicists as even near threshold energies it couples strongly to other strangeness -1 baryon bound and reaction states. Since the Bohr radii for kaonic hydrogen is about 1000 times smaller than in electronic atoms, the orbiting hadron lies close enough to the nucleus as to experience its nuclear strong force as well as the Coulomb interaction. Because of this proximity there can be a varying degree of overlap between the wavefunctions for the antikaon and proton which cause these additional channels to open: \begin{equation}\label{chan} K^-p \longrightarrow \left\{ \begin{array}{lll} K^-p&+8.6KeV&atom\; bound \\ \Sigma^{\pm,0} \pi^{\mp,0}&+100MeV&annihilation\; open \\ \Lambda^0 \pi^0&+180MeV&annihilation\; open \\ \overline{K}^0n&-5MeV&charge\; exchange\; closed \end{array} \right. \end{equation} To quantize the exact energies above, a program written specifically for this purpose was employed \cite{rlfor}. Containing terms representing both the Coulomb and nuclear force, the Coulomb force is specifically computed with a large radius cuttoff so that a nonsingular momentum space formulation can be dealt with. This momentum space solution was invoked with an ingenious subtraction technique pioneered by Lande, Kwon and Tabakin \cite{kt}. \section{Solving the Schr\"{o}dinger Equation} For any given system, eigenstates, and their accompanying energy eigenvalues, are traditionally retrieved via the solution of the time-independent Schr\"{o}dinger equation: \begin{equation}\label{tisch} H\Psi=E\Psi. \end{equation} Splitting the Hamiltonian into the accustomed kinetic and potential parts: \begin{eqnarray}\label{spsch} (H_o+V)\Psi&=&E\Psi,\\ H_o&=&{\mathbf p}^2/2m,\nonumber\\ {\mathbf p}&=&-i\hbar{\mathbf \nabla},\nonumber \end{eqnarray} allows for the rewriting of this operator equation as the well recognized differential one: \begin{equation}\label{dscheq} \left(-\frac{\hbar^2}{2m}{\mathbf \nabla}^2+V\right)\Psi=E\Psi. \end{equation} Multiplying through by $2m/\hbar^2$ and making the simple substitution: \begin{equation}\label{usubs} U=\frac{2mV}{\hbar^2}, \end{equation} while remembering that: \begin{equation}\label{esubs} E=\frac{\hbar^2k^2}{2m}, \end{equation} the Schr\"{o}dinger equation can be cast into the form of a nonhomogeneous Helmholtz equation: \begin{equation}\label{helmeq} \left({\mathbf \nabla}^2+k^2\right)\Psi=U\Psi. \end{equation} Nonhomogeneous because the term on the right side is nonvanishing $(U\not=0)$. The advantage here is that the Helmholtz equation is a special function, the solutions of which have been studied in great detail and well documented \cite{arf}. \subsection{Green's Function} If $\varphi_n$ is taken to be an eigenfunction satisfying the homogeneous Helmholtz equation: \begin{equation}\label{hhelmeq} \left({\mathbf \nabla}^2+k^2\right)\Psi=0, \end{equation} then the solution to the nonhomogeneous equation (\ref{helmeq}) may be written in terms of an integral equation of a Green's function $G$: \begin{equation}\label{inteqg} \Psi(k;{\mathbf r})=\varphi(k;{\mathbf r})+\int G(k;{\mathbf r} ,{\mathbf r}')U({\mathbf r}')\Psi(k;{\mathbf r}')\,\mathrm{d}{\mathbf r}. \end{equation} Here the Green's function is defined to be that symmetric function which satisfies the point source equation: \begin{equation}\label{gdef} \left({\mathbf \nabla}^2+k^2\right)G(k;{\mathbf r},{\mathbf r}') =\delta({\mathbf r}-{\mathbf r}'). \end{equation} Point source since $\delta$, the Dirac delta function, is nonzero only at the point ${\mathbf r}={\mathbf r}'$. Because the eigenfunctions of the hamiltonian, which is a hermitian operator, form a spanning set, the Green's function may be expanded in a series of these eigenfunctions of the homogeneous equation (\ref{hhelmeq}): \begin{equation}\label{gexpan} G({\mathbf r},{\mathbf r}')=\sum_{n=0}^\infty a_n({\mathbf r}')\varphi({\mathbf r}). \end{equation} The Dirac delta function is represented by its eigenfunction expansion: \begin{equation}\label{dexpan} \delta({\mathbf r}-{\mathbf r}')=\sum_{n=0}^\infty \varphi_n^*({\mathbf r}')\varphi_n({\mathbf r}). \end{equation} By inserting the preceding two expansions into the point source formulation (\ref{gdef}) and applying the orthogonality condition $\int \varphi_m^*\varphi_n=\delta_{mn}$, an expression for the coefficients $a_n$ may be isolated and substituted back into (\ref{gexpan}) to yield: \begin{equation}\label{gsum} G({\mathbf r},{\mathbf r}')=\sum_{n=0}^\infty \frac{\varphi_n^*({\mathbf r}')\varphi_n({\mathbf r})}{E-E_n}. \end{equation} As the continuum of scattering states where the potential vanishes are represented by plane waves, $\varphi_n$ may be so written. Replacing the sum in (\ref{gsum}) by an integral, the Fourier-integral representation of the Green's function takes form: \begin{equation}\label{gfir} G({\mathbf r},{\mathbf r}')=\frac{1}{{2\pi}^2}\int \frac{e^{i{\mathbf k}'\cdot({\mathbf r}-{\mathbf r}')}}{E-E_n}\,\mathrm{d}^3k'. \end{equation} Evaluation of this integral begins with replacement by a complex contour integration \begin{equation}\label{zintx} \oint_c f(z)\,\mathrm{d}z \longrightarrow \int_{-\infty}^{+\infty} f(x)\,\mathrm{d}x \end{equation} closed at infinity so that Cauchy's integral formula may be applied. The fact that the two poles lie directly on the real axis may be remedied by Feynman's prescription of assuming the energy to have a small positive imaginary part $+i\epsilon$: \begin{equation}\label{igfir} G({\mathbf r},{\mathbf r}')=\frac{1}{{2\pi}^2}\int \frac{e^{i{\mathbf k}'\cdot({\mathbf r}-{\mathbf r}')}}{E-E_n+i\epsilon}\,\mathrm{d}^3k'. \end{equation} This integral may now be equated to a principle value part in addition to a part proportional to the residues and ultimately substituted back into equation (\ref{inteqg}) for solutions. \subsection{Matrix Form} An alternate understanding can be gained from Heisenberg's matrix formulation of quantum mechanics. Here operators can be portrayed as simple matrices and state functions as vectors. Beginning with the matrix operator form of the Schr\"{o}dinger equation (\ref{tisch}): \begin{eqnarray}\label{mmeqsch} \left[H\right]_{nn}\left[\Psi\right]_{n1}&=&\left[E\right]_{nn}\left[\Psi\right]_{n1},\\ \left[H_o+V\right]_{nn}\left[\Psi\right]_{n1}&=&\left[E\right]_{nn}\left[\Psi\right]_{n1},\nonumber\\ (\left[H_o\right]_{nn}+\left[V\right]_{nn})\left[\Psi\right]_{n1}&=&\left[E\right]_{nn}\left[\Psi\right]_{n1},\nonumber\\ \left[H_o\right]_{nn}\left[\Psi\right]_{n1}+\left[V\right]_{nn}\left[\Psi\right]_{n1}&=&\left[E\right]_{nn}\left[\Psi\right]_{n1},\nonumber \end{eqnarray} a rearrangement is possible: \begin{eqnarray}\label{rsch} \left[V\right]_{nn}\left[\Psi\right]_{n1}&=&\left[E\right]_{nn}\left[\Psi\right]_{n1}-\left[H_o\right]_{nn}\left[\Psi\right]_{n1},\\ \left[V\right]_{nn}\left[\Psi\right]_{n1}&=&(\left[E\right]_{nn}-\left[H_o\right]_{nn})\left[\Psi\right]_{n1},\nonumber\\ \left[V\right]_{nn}\left[\Psi\right]_{n1}&=&\left[E-H_o\right]_{nn}\left[\Psi\right]_{n1}.\nonumber \end{eqnarray} Expanded, this reads as (note $\left[E\right]_{nn}\equiv E\left[1\right]_{nn}$ since energy is an eigenvalue): \begin{equation}\label{mrsch} \left[ \begin{array}{cccc} V_{11} & V_{12} & \cdots & V_{1n} \\ V_{21} & V_{22} & \cdots & V_{2n} \\ \vdots & \vdots & & \vdots \\ V_{n1}&V_{n2}& \cdots &V_{nn} \end{array} \right] \left[ \begin{array}{c} \Psi_1 \\ \Psi_2 \\ \vdots \\ \Psi_n \end{array} \right]= \left[ \begin{array}{cccc} E-H_{o_{11}} & H_{o_{12}} & \cdots & H_{o_{1n}} \\ H_{o_{21}} & E-H_{o_{22}} & \cdots & H_{o_{2n}} \\ \vdots & \vdots & & \vdots \\ H_{o_{n1}} & H_{o_{n2}} & \cdots & E-H_{o_{nn}} \end{array} \right] \left[ \begin{array}{c} \Psi_1 \\ \Psi_2 \\ \vdots \\ \Psi_n \end{array} \right]. \end{equation} Multiplying both sides through by the inverse of the operator (matrix) found on the right gives: \begin{equation}\label{dsch} [E-H_o]_{nn}^{-1}[V]_{nn}[\Psi]_{n1}=[\Psi]_{n1}. \end{equation} Finally, factoring forces this into the form of a system of linear homogeneous equations: \begin{equation}\label{fsch} [1-(E-H_o)^{-1}V]_{nn}[\Psi]_{n1}=0. \end{equation} The inverse of the operator $(E-H_o)$ is just the inverse of the matrix representing it and is identified as the Green's function presented earlier: \begin{equation}\label{gv} [1-GV]_{nn}[\Psi]_{n1}=0. \end{equation} The identity operator `1' is simply the identity matrix. Nontrivial solutions for a homogeneous system are possible only when the determinant vanishes: \begin{equation}\label{crit} \rm det (1-GV) \equiv \left|1-GV\right|=0. \end{equation} \subsection{Analytic Continuation} Post evaluation of the Green's function (independent of the Feynman's prescription) the domains of the final forms of the previous equations are extended into the complex plane. This process whereby the independent variable ($E$) is made complex ($E_r+iE_i$) is called analytic continuation \cite{arf}. By so assuming, additional information is provided through the extra degree of freedom; specifically on transitions or growth and decay of states which is related to twice the imaginary part of the energy satisfying (\ref{crit}) and which is denoted $\Gamma=2|E_i|$. More is said of this in the interpretation of results section. \section{The Code} \subsection{File Read Parameters} A majority of the parameters needed during runtime are read in from a text file. This allows the user to conveniently make changes between consecutive runs without altering the code. The input file must have the name {\tt runbs} and reside in the same directory as the executable. A typical example of {\tt runbs} would appear as follows: \begin{verbatim} -0.008715 -0.00080 -0. -0. -0. -0. -0. -0. 0.0001 64.0 0 1 40 -8 -2.4 00 100. 2400. 1.55 1.55 0.00 1.55 1.55 0.00 3.0 7. 1 1 10 +0 0 1 7 0 0 0 0 5 0 0 0 0 +2 0 0 0 0 0 002 +14 0. 0. 0. 0. \end{verbatim} representing the coded variables: \begin{verbatim} drr(1) eim drr(2) drr(3) drr(4) drr(5) drr(6) drr(7) tol aitmax nr lxmax ngp kode catom nnuc csize cnuc achp acmp wsp achn acmn wsn rcoul rcut nz na nifty(1) nifty(2) nifty(3) ... nifty(20) nes nwaves b0r b0i c0r c0i \end{verbatim} Decimal points that appear not to be needed actually indicate to the program that the number is to be interpreted as a floating point value. Some of the important entries above are: \begin{itemize} \item {\tt drr(1)}: Initial guess at the real part of the energy that the search starts with. Subsequent guesses for {\tt drr(n)} that are nonzero will cause the program to loop over these additional n-1 energy estimates. \item {\tt eim}: Initial guess at the imaginary part of the energy that the search starts with. \item {\tt tol}: Calculations are carried out to the desired tolerance for $\left|1-GV\right|=0$ search. \item {\tt aitmax}: Maximum number of iterations in search. If no convergence, break in flow with error code. \item {\tt ngp}: Number of grid points used in the solution of the integral equation. Determines size of array. \item {\tt catom, csize, cnuc}: Gauss integration points are distributed over the nuclear and Coulomb region using these values. The idea is to distribute one set of points in the low momentum, atomic, region where the Coulomb potential is thought to be more important and another set in the larger momentum, nuclear, region where the strong nuclear potential is most important. {\tt csize} determines the cutoff on the atomic region. {\tt catom} and {\tt cnuc} are the distribution parameters for the atomic and nuclear regions respectively. The sign on {\tt catom} acts as a switch to specifies a distribution method but does not otherwise effect the scaling value. \item {\tt achp, acmp, wsp}: Proton distribution size parameters. \item {\tt achn, acmn, wsn}: Neutron distribution size parameters. \item {\tt rcoul}: The size (radius in fm) of the charged sphere for the Coulomb potential. \item {\tt rcut}: Cut of radius that marks the boundary of the physical model. \item {\tt nz}: Nuclear charge number. \item {\tt na}: Nuclear mass number. \item {\tt nifty(1)}: Specifies type of particle. A value of `10' indicates antikaon-proton pair. \item {\tt nifty(3)}: Specifies quadrant in which to search for poles; determines bound or resonance states. \item {\tt nifty(6)}: Specifies forces included in calculation. \item {\tt nifty(15)}: Specifies the number of open channels. \end{itemize} \subsection{Execution} Before the code can be run, the source code must be compiled into object files and those object files must then be linked to form the executable program. This is often done in one step and loosely referred to as simply `compiling'. Typing the command:\\ \\ {\tt f77 pear.f -o bsprog}\\ \\ at the prompt from within the directory containing the source code {\tt pear.f} will form an executable program called {\tt bsprog}. A user defined name, such as {\tt bsprog}, is allowed when the switch {\tt -o} is used. The program may now be run by simply entering the text {\tt ./bsprog} at the prompt. A more popular method of compiling is to leave all the subroutines listed separately as differently named source files (this can be accomplished by applying the {\tt fsplit} command to {\tt pear.f}). Compile them all together forming object files, and then to link all the object files to form a load module. The advantage here is that if only a small change is made to one of the subroutines, then only that subroutine would have to be recompiled before linking. From within the directory containing all the separate source files typing:\\ \\ {\tt f77 -c *.f}\\ \\ at the prompt will make object files with the {\tt .o} extension. The {\tt -c} switch instructs the compiler to compile only without calling the linker and the * is a wildcard taking the place of all the filenames with the {\tt .f} extension. Linking is now accomplished with the command:\\ \\ {\tt f77 -o bsprog *.o}\\ For further options available under the {\tt f77} compiler, or instructions in general, students may access the manual by entering {\tt man f77} or consult additional reference material \cite{lgsc}. \subsection{Interpretation of Results} Upon termination of a run, a file called {\tt outbs} it written to disk. The top portion of the file is merely an echo of input parameters, some of which were discussed above. Echoing is a valuable technique that shows that these data were read in correctly as well as insuring that intermediate results appear as expected. It also provides a record of information that accompany results for later referral. The blocks in which {\tt zdet=} can be seen in the lower right show the progress of the search. Just to the left in these same blocks the energy iterations appear. The very bottom of {\tt outbs} will read as: \begin{verbatim} 0E(MeV)= -.008276 -.000030 eps,gamma(eV)= -336.2 60.1 Del E(ev)= 336.2 -30.1 \end{verbatim} {\tt 0E(MeV)=} are the real and imaginary parts of the energy, within tolerance, that satisfy $\left|1-GV\right|=0$. The {\tt Del E(ev)=} term gives the difference between the exact calculation offered by the code and the classically calculated Bohr energy. Notice the imaginary part here is redundant since the Bohr model has no imaginary part built into it. $\Gamma$ (gamma) is twice the imaginary part and taken to be always positive. Mathematically, both positive and negative complex parts are possible but the algorithm searches out only the negative values on physical grounds. Positive values will lead to states that grow with time as can be inferred from what follows. \begin{equation}\label{egama} E=E_r+iE_i=E_r \pm i \Gamma /2 \buildrel \rm neg \over \longrightarrow E_r - i \Gamma /2. \end{equation} $\Gamma$ is a byproduct of the analytic continuation and provides stability information, being inversely (for negative values) proportional to the mean lifetime of the state. These states for which the energy lies off the real axis in the complex plane ($\Gamma \not= 0$) and have a finite lifetime are referred to as {\tt resonance} states. From the uncertainty relationship: \begin{equation}\label{uncer} \Delta E \Delta t = \frac{\hbar}{2}, \end{equation} it is understood that the energy uncertainty varies inversely as the uncertainty in time. It is also true that the time dependance of a quantum state goes as: \begin{equation}\label{tdep} \Psi(t) \propto e^{-iEt/\hbar}=e^{iE_rt/\hbar}e^{-\Gamma t/2\hbar}, \end{equation} which leaves the probability density decaying in time: \begin{equation}\label{pdtdep} |\Psi|^2 \propto e^{-\Gamma t/\hbar}=e^{-t/\tau \hbar}, \end{equation} with a time constant, or lifetime: \begin{equation}\label{tconst} \tau=\frac{\hbar}{\Gamma}. \end{equation} Thus a larger $\Gamma$ corrosponds to a state with a shorter lifetime and thus a greater energy uncertainty by (\ref{uncer}) \cite{RLQMII}. \section{Exploration} \subsection{Source Code Libraries} Some common tasks come up repeatedly in scientific programming. Solving a system of equations, elementary function evaluation, finding eigenvalues, interpolation and extrapolation are just a few examples. Such tasks may be written once and recycled in effect. The idea is to have these tasks written as subroutines by experts so that they are well optimized and yield superlative results. The order in which a matrix is accessed in memory for instance (row or columns first) will have a bearing on how fast a routine runs as well as many other factors. These routines are then archived in a library for the public to access in include in there own larger programs. Some minor modifications may be necessary to weave these subroutines seamlessly into a programmers source code, but the savings can be significant in both programming and run time. It is not uncommon to have execution times dropped by an order or magnitude. That can mean running a day instead of a month! This also translates to a dollar amount if one is working for a company or renting computer time. The program {\tt pear.f} incorporates subroutines from the {\tt SLATEC} library. This is an acronym for Sandia, Los Alamos, Air Force Weapons Laboratory Technical Exchange Committee. These subroutines can be identified by searching through the source code for the string {\tt (sla)} which resides in the title block. The function of these routines as they exist within {\tt pear.f} is to aid in the advanced search for the zeros of the determinant $\left|1-GV\right|=0$ as well as providing error codes for information upon abnormal termination. Many such libraries, in addition to a wealth of additional information, can be accessed via the world wide web at {\tt www.netlib.org}. Choose the option of browsing the Netlib repository. There Slatec can be found; as well as Linpack, Lapack, Svdpack and many others. Explore these resources and then try downloading library routines for use in a program of your own design. \subsection{Find Your Way Home} The search that is employed by {\tt pear.f} is nonlinear. That is to say that initially, starting with an estimate, it will take large steps. As it draws nearer to the convergence point, the search will actually switch methods. This is all done in an effort to make the routine as efficient and timely as possible. This is a two parameter search as it must simultaneously seek out both the real and imaginary, $E=E_r+iE_i$, parts of the energy that satisfy (\ref{crit}) under analytic continuation. It is of interest to watch the convergence of the routine. To gain an understanding of its limitations. For some poor initial estimates, or other possible circumstances, the routine may not converge at all. To this end, the code {\tt pear.f} will be modified to dump the search results to a file for plotting. The following two changes are to be made in order so that the statement line numbers will agree with instructions. Changes should also always be applied to a saved copy of the source code to preserve the working original listing. \begin{enumerate} \item Insert after line 20: {\tt open(4,file='plot.out',status='unknown')} \item Insert after line 1236: {\tt write(4,*) zdet} \end{enumerate} After recompiling and running the code, a file called {\tt plot.out} should appear within the directory containing two columns of numbers. The first column represents consecutive iterations of the real part of the energy. The second column represents the imaginary part. These columns should each be plotted separately against line number. Try this several times for different initial guesses. What if multiple values for $E=E_r+iE_i$ are possible that satisfy (\ref{crit})? \subsection{Plain Old Hydrogen} When it seems that a scientific code has been completed and is ready for use, it is always a good idea to test the output of the program. The code should be modified as to calculate a known quantity and that result should then be compared against the accepted value. Needless to say, the values should agree to within an error that is determined by the model used and the computational means. With this agreement intact, an additional confidence in the code is instilled and subsequent results are more prone to be trusted. The code {\tt pear.f} may be modified to yield results for plain old garden variety hydrogen. Here a lone electron orbits a proton. Nuclear forces are ignored and additional channels are closed. The following changes reflect this modification. \begin{enumerate} \item {\tt drr(n)}=n$^{th}$ initial guess at energy \item {\tt eim}=0.0; Bohr model lacks imaginary energy \item {\tt catom}=.01; sets smaller scale factor for grid \item {\tt nifty(1)}=10; sets $k^-$ as orbiting particle \item {\tt nifty(6)}=6; uses only Coulomb force \item {\tt nifty(15)}=0; leaves only 1 channel open \item Replacement in line 759: {\tt 493.699d+00 $\rightarrow$ 0.511003d+00} \end{enumerate} Since no electron option is available, the switch {\tt nifty(1)} is set for an antikaon and its rest mass is adjusted to match that of an electron (0.511003 MeV). As the calculation is done in momentum space (Fourier transform of coordinate space taht scales inversely) the scale factor must be reduced to compensate for the larger orbital radius. The output values are then to be contrasted with the well known energies given by the Bohr model for Hydrogen (-13.60, -3.40, -1.51, -0.85 eV). Notice that the initial guesses that start the search must be adjusted `close' to these values so the search does not deviate to neighboring energies. What could be expected if a search for the n$^{th}$ arbitrary quantum level were carried out for large values of n? \subsection{Kaonic Hydrogen Hold the Nuclear} Making now the solitary change of bumping up the rest mass of the electron back up to that of an antikaon, the effects of transforming simple hydrogen into kaonic hydrogen can be studied one step at a time. Remembering that the particle type ({\tt nifty(1)=10}) is already specified as an antikaon, \begin{enumerate} \item restore original antikaon mass ({\tt 493.699d+00} MeV) in line 759. \end{enumerate} \subsection{Add the Strong Interaction} Because of the much tighter Bohr radius for the more massive antikaon, as compared to the electron in the hydrogen atom, it is expected that the incorporation of the nuclear strong force play a crucial role and should change results appreciably. \begin{enumerate} \item {\tt nifty(6)}=0; uses the strong in addition to the Coulomb force \end{enumerate} \subsection{Open Additional Channels} The term channel is an analogy that has crept into the vernacular of quantum mechanics. The parallelism is that the fluid carried by these channels is that of the probability density and must amount to a "volume" of unity. An additional channel opening represents an additional path for the probability fluid to follow and a different possible state or outcome. The state's chance for success is directly proportional to the fraction of the total amount of fluid flowing through that particular channel. The exchange of an electron leading to the formation of two different particles is one example of such a possible channel here and is also the last given in expression (\ref{chan}). To open additional channels within the computational model, the switch {\tt nifty(15)} needs to be adjusted. It can take on values {0,1,2,3} and different configurations should be tried. To see all the values that {\tt nifty(n)} can accept the reader is referred to the coded lines of the source numbering 761-785. Also the output file ({\tt outbs}) will list the choices made for {\tt nifty(n)} and a short abbreviation of its meaning. As the additional channels are opened look for a drop (more bound) in energy. \begin{thebibliography}{99} \bibitem{rlfor} R. H. Landau, \textsl{Coupled Bound and Continuum Eigenstates in Momentum Space} Phys. Rev., {\bf C 27}, 2191--2197 (May 1983). \bibitem{kt} Y. R. Kwon, F. Tabakin, Phys. Rev., {\bf C 10}, 391 (1978). \bibitem{arf} G. Arfken, \textsl{Mathematical Methods for Physicists Third Edition,} Academic Press Inc., San Diego, 915--920 (1985). \bibitem{RLQMII} R. H. Landau, \textsl{Quantum Mechanics II,} John Wiley \& Sons Inc., New York, 49--104 (1996). \bibitem{lgsc} R. H. Landau, P. J. Fink Jr., \textsl{A Scientist's and Engineer's Guide to Workstations and Supercomputers} John Wiley \& Sons Inc., New York, 179--217 (1993). \end{thebibliography} \end{document}