Nuclear Physics A491 (1989) 525-540 North-Holland, Amsterdam
THE
BOUND
STATES
MAGNETIC
OF AN ELECTRON MONOPOLE
R. SINGER Insriru?
and
AND
IN THE
FIELD
OF A
A PROTON
D. TRAUTMANN
fiir Physik, Universitiit Basel, CH-4056 Basel, Switzerland G. BAUR
Institut fi& Kernphysik,
Kernforschungsanlage
Jiilich, D-51 70 Jiilich, BRD
Received 14 April 1988 (Revised 24 August 1988) Abstract: Bound states of an electron in the field of a proton and a magnetic monopole, separated by a distance R, are studied in the non-relativistic Pauli approximation. The wave functions and the energy levels are obtained by a finite difference analysis of the corresponding Schrodinger equation. The agreement of the results with previous calculations is discussed.
1. Introduction Over
fifty years
ago,
Dirac
introduced
the magnetic
monopole
according
to
theoretical arguments ‘). Already Dirac ‘) and Banderet ‘) recognized that there exists no bound state of the monopole-charge system in the non-relativistic approximation. There are a number of papers which investigate the bound states of a charged particle in the field of a monopole and an additional charge Ze [refs. ‘-‘)I. Most of these works treat the monopole and the charge Ze as dyon, that is the particles are both situated in the origin. Recently, the observation of a possible has intensified the interest in this particle.
candidate
for the magnetic
monopole
“)
In the past few years, new experimental methods for detecting magnetic monopoles were invented ‘-12). They are based on the energy loss of a slowly moving monopole in matter (H- or He-gas). To understand the mechanism of energy loss one has to investigate the possibility of binding chemically a monopole to an atom or a molecule. The first to treat this problem was Malkus ‘). He estimated the energy of an electron in the field of a monopole and a proton, separated by a distance parameter R, by interpolating between two extreme cases: (i) Explicit calculation of the electron energy in the field of the dyon (i.e. R = 0). (ii) Calculation of the diamagnetic and paramagnetic energy of the electron for large R. Malkus suggested that chemical binding might be possible. Five years ago Drell et al. investigated the same problem i3*i4). He obtained his results by piecing together two solutions in various regions: 03759474/89/%03.50 @ Elsevier Science Publishers (North-Holland Physics Publishing Division)
B.V.
R. Singer et al. / Bound .staw.s of’ an electron
526
(i) For R = 0 the energy levels are explicitly levels are calculated
by a multipole
(ii) For large R the Zeeman Drell
also concluded
chemically
splitting that
known 3). For small values of R the
expansion. calculation
a neutral
is used.
magnetic
monopole
could
be bound
to atoms or molecules.
In order to find the stability of these chemically bound states against transition to the direct monopole-nucleus bound states it is necessary, in the simplest case, to be able to calculate exactly the bound states of an electron in the field of a monopole and a proton separated by a distance parameter R. Ruijgrok was the first who approached this problem by a numerical method “.16). Based on the wave functions of Kazama Dirac equation numerically and obtained
et al. “) he integrated the corresponding in this way the ground state for various
parameters R. But the Dirac theory has an enormous disadvantage: In the ground state the hamiltonian is no more well-behaved; there exist no nonvanishing solutions of the Dirac equation that fulfil the finite boundary condition at the origin I’). Kazama et al. resolve the problem by adding an infinitesimally small extra magnetic moment K
to the electron.
This
effect
prevents
the
electron
from
passing
through
the
monopole I’). Ruijgrok et al. obtained their results for positive as well as for negative values of K [ref. “)I. The shape of the monopole-atom potential for the lim K +O is quite obvious. Furthermore, Goldhaber ‘*) and Callias “) found that there exists a oneparameter ambiguity in the boundary condition of the ground state. To avoid these problems in this paper we make use of the non-relativistic Pauli approximation. Thus, we will neglect the spin-orbit interaction and all terms of higher order, but we will not neglect the electron spin. We examine the finite difference method to calculate in principle the complete discrete spectrum of the corresponding
Schrodinger
equation.
2. Formulation
of the problem
We are seeking to calculate some of the lowest bound states with the corresponding eigenfunctions of an electron in the field of a monopole and a proton. Let the proton p be fixed in the origin and the monopole Q at a distance R on the negative z-axis. Denote by r, the proton-electron distance and by rz the monopole-electron distance (fig. 1). The Schrodinger equation due to the non-relativistic Pauli approximation reads explicitly
(1)
R. Singer et al. / Bound states of an electron
527
Fig. 1. The geometrical arrangement of the problem in Cartesian and spherical coordinates. We use Born-Oppenheimer approximation, that is the proton p and the monopole Q are fixed at their positions.
In the singularity-free spherical coordinates
formalism of Wu and Yang 20,21) the vector ri, 8, p is given by
(A,), = (A,), = 0, (ArIp = (&t)p = 0,
(A,),
=t
R,:
a
1 regions
os?Y<;n+ti,
defined
O
&r-tk<+T,
O
we can require
in R,,
in R,,
A, in
(24
W)
by
o
OCcp<2T, which
E=o,i,1,5 examinations
r2+f;;;s
(A,),+
g is the magnetic charge of the monopole charge e by the quantisation rule ‘)
For further
a
1
where R, and Rp are the two overlapping R,:
r2-;;;;s
potential
o
must be related
to the electric
.
,... rotational
symmetry
around
the z-axis and
make the ansatz
(f(rl,fi)e-‘“‘P >
*=
a,bEN+{O},
g( r, , 6) emibq ’
(3)
meaning 9 = L1
f(r,, 6) e-iO~O g(r,
,
>
6) ebibmv
in R,
(34
in RD.
(3b)
and tie
=
f(rl, 6) eP@’ g( rl , 6) eeib@
>
R. Singer
528
et al. / Round States of an electron
In order to cancel out the Q-dependence a and b must be related
in the two resulting
Schrodinger
equations,
by a-b=l.
Now we have to deal with the boundary
(4) conditions
for f and g. First we make
use of the symmetry of the wave function $. After eliminating the dependence on Q we mention that f and g in region II must be either symmetric or antisymmetric to f and g in region I (fig. 1). Thus we only need to solve the Schrodinger equations in the region 0 enclosed of the boundary a0 defined by 6=O,r;
r1= 10,031
6 = [O,?T];
9
r,=cD.
(5)
According to the singularity of the Laplace operator A on the z-axis this boundary is no more well-behaved: The wave equation is singular on a part of a0 given by r, ,Z 0 -9 = 0, 7r. To avoid this problem we make for f and g the ansatz
f(r,, 6) =
-&F”,,w,
Because f and g must be finite on fi the boundary given by
conditions
for F and G are now
F(O,6)=F(r,,O)=F(r,,n)=O, G(O,6)=G(r,,O)=G(r,,rr)=O,
(7a)
where -9 E [0, P] and r, E [O,OO] . Moreover,
it is easy to show that F and G must tend to zero at infinity:
lim F( r, ,9) = 0,
r, -02
lim G(r,,
,,-02
6)=0.
(7b)
Instead of this boundary condition which is not practicable a modification of (7b) defined by F(r,6)=G(r,6)=0, By setting
(3) and (6) into (1) we obtain
on computers,
we choose
for l
(7c) equations
for F(r,,
6)
529
R. Singer et al. / Bound states of an electron
and G(r,,
6) as
-${;}-+${;}+~$(GF)
m
G3-1 +*eg
+
[ rf
hc r2(r2+r1 cos 6+R)
sin’ 6
eg R + rl cos 6 x
ri
where the distances
2
(H
+5X
2 --+s rl
hc
r,(r,+r,
r,
sin6 cos 6+R)
r,, r2 and R are to be taken in units of the Bohr radius
2
1
a, and
the energy E is measured in units of the ground state of the neutral H-atom (= 13.6 eV). s > 0 is a shift parameter to make the eigenvalues A positive so that we can apply the ideas described in sect. 3. 3. Numerical In order to determine
the eigenvalues
treatment
and the corresponding
eigenvectors
of the
Schrodinger equation (8) we follow the idea by Neuberger and Noid 22). Suppose A is a symmetric positive definite linear transformation on a real finite dimensional inner product space U. Now we use a power method to find the largest eigenvalues with the corresponding eigenvectors of A [ref. “)I. Choosing an arbitrary vector u,, E U, u0 # 0, we can define the sequence {u,}z, by the relation (9) Thus we find u, 0~ A$, . Now we can expand u0 in terms of a complete
set of eigenvectors
{a,}~=, of A:
n
Ug =
C
(10)
(Yiai
i=l
so that u, is given by (apart f
from the normalisation
~,I:ai=h;[n,a*+~2~i(~)‘ai].
factor) (11)
i=l
If IhJ > )A212 /A313 * . * 2 IA,], then we have a, cc lim 24, s-00 and therefore,
because
(12)
of (9), (13)
If A, is found we choose u. orthogonal to an eigenvector associated with A, and repeat process (9). The limits (13) and (12) then yield the second largest eigenvalue A2 and the corresponding eigenvector a2 (or again Al if it is degenerate). This process
530
R. Singer et al. / Bound states of an
electron
can be repeated as often as is needed to calculate in principle the complete discrete spectrum of A. In most physical problems there are the smallest eigenvalues which are of interest and not the larger ones. Now the smaller eigenvalues of A are equal to the reciprocal values of the larger ones of A-‘. But the finite difference approximation of A is mostly a sparse matrix whereas the inverse of such a matrix is typically full. Now the number of grid points is in the range of about 10’ so that it is clear that matrices corresponding to a finite difference approximation of A-’ cannot be processed on currently available computers. We handle this problem in the following way: Starting from the eigenvalue equation Au = hu with certain boundary conditions, let G., be the grid of the corresponding finite difference approximation. Denote by U the vector space of all real-valued functions on Gfd which fulfil the boundary conditions on boundary points and by Afd the finite difference approximation of A. Thus we have in the grid point P the equation (A~~~~~)(P)= (h~~d)tP)=: w(P),
(14) where U, w(P) E U and PE Gfd. Now in order to compute Azw(P) we solve for ufd in (14) by use of the successive over-relaxation method 24). In what follows we describe how to solve the two coupled partial differential equations (8) on fi with the boundary conditions (7a) and (7~) on $0 by the finite difference method. We pick a radiai grid size &, and an angular grid size S6 and denote by Gfd the corresponding grid on d (fig. 2). According to (14) the finite difference approximation of (8) can be written as uJd
E
-{
oFjr+l,j-j GF}j_l,j-‘~{ ;:GF}i,j-df3j{ ;}j,j={ ::}i,jy GFJ,+l-“t(
+{
GF}i,j-1
(15)
where
1
*
cos 29
c~‘=(i-I)2~~2*2(i-I)2~~sin~~ eg Srfr, sin 6
4.j = hc
r:
r,= (It-i- R2+2r,R
{
3 cos 6)“‘,
::}i,j=srq gF)Ij*
r,=(i-l)Sr,,
s=(j-l)s?Y,
(16)
531
R. Singer et al. / Bound states of an electron
Fig. 2. The grid G, on 6. The index i counts in radial direction whereas the index j runs in angular direction. For numerical stability I%, should be chosen in such a way that one has R =$2~ + l)Sr,; n=0,1,2 ,...
Using
the successive
over-relaxation
method
in solving
for n@ in (14) yields the
relations
where at each step the value of each term is the best current the acceleration parameter and ranges from 1 to 2.
approximation.
w is
4. Results of calculations In order to obtain
satisfactory
results by using the finite difference
described above, two main questions arise: (i) How large does one have to choose r in the boundary
approximation
condition
(7c)?
(ii) How small does one have to take the grid size Sr, and c%?? Numerical examination yields that for r C=lSa, the results did not change further with varying r. Thus, we choose r = 15ao. The same condition leads to the choice of Sr, = O.O6a, and 66 = q/92. The value for r turns out to be nearly independent of the separation parameter R as one expects because the binding of the electron to the proton is much stronger than the influence of the magnetic monopole on the electron. Moreover, we know the approximate value for the ground state, which is s -1. Thus, in order to make all eigenvalues positive, we set the shift parameter s equal to 2. For the acceleration parameter o the best value turns out to be 1.962. With this choice of w the computing time could be reduced by a factor of up to ten compared with the simple Gauss-Seidel method (w = 1).
R. Singer et al. / Bound states of an electron
532
The results lowest
for the bound
energy
magnetic electron
levels
charge angular
the electron
states are shown
as function
in fig. 3. This figure exhibits
of the separation
parameter
the five
R for the smallest
g given by eg/ hc = 5. The assignment
of the z-components
of the
momentum
to Drell’s reflections
about
spin-flip
(m,) was done
when a monopole
For R = 0 we obtain, as already the field of a dyon given by
found
according
passes
from z = -cc to z = +co [ref. “)I.
by Malkus 3), the states of an electron
in
(18) where
n = 1, A, 2,. . . . The first three excited
states (n = Ji)
are degenerate.
For large R (R 5 lOa,) the influence of the magnetic monopole is expected to nearly vanish. We actually find the energy levels to be very close to those of the neutral H-atom. Moreover, for large values of R the electron undergoes a magnetic field which is nearly homogeneous and parallel to the z-axis. In this case the Zeeman
-0.2
-0.4
-
0 E2
m. z-1 I
.
H-
2
E2
G
-0.6
t
\
i
-0.0
-1.0
-1.2
Fig. 3. The ground state and a proton separated energy (- 13.6 eV) and left-hand side (ED) are the levels of the neutral
and the first four excited states of an electron in the field of a magnetic monopole by a distance R. The energy E is applied in units of the H-atom ground-state the separation R is measured in units of the Bohr radius a,. The states on the the energy levels of an electron in the field of a dyon, on the right-hand side H-atom are shown (EH). The dotted curves exhibit the Zeeman splitting of the ground-state energy of the neutral H-atom (E 7).
R. Singer et al. / Bound states of an electron
energy
shift is easily calculated /jE
z
by 1 71 I R
l+~(~+l)+s(s+l)-l(r+l)
=tnti
{
Nj+l)
where R must be given again in units of the Bohr radius levels calculated
533
in this way are found to join smoothly
(19)
a,. For large R the energy the results of our numerical
calculations. For example, in the ground state the magnetic field of the monopole can be treated locally as homogeneous down to a distance of R = 3a,; the deviation is less than 1 p.c. (fig. 3). We find a minimum in the ground
state of E = -1.135
at a separation
of R = 0.95ao.
This is much less than the rough estimate of Malkus, who mentioned a minimum of E = -1.515 at a separation of R = 0.56a. [ref. 3)]. The agreement of the present ground state with that calculated by Ruijgrok is excellent: He computed a minimum of E = -1.13 at a separation of R = 0.9a, [refs. is*‘“)]. Thus, we mention that a magnetic monopole can be bound chemically to a neutral H-atom in the ground state. A further bound state which is indeed very weak can occur in the fourth excited state at a, c R =S2a, (fig. 3). In fig. 4 the vicinity of the dyon energy levels Ef and Ez are shown in a more detailed
representation.
From
an extrapolation
in R (R + 0) we can calculate
the
derivative of E as a function of R for R = 0. For the ground state this derivative vanishes (fig. 4, left). This result coincides exactly with that of the corresponding perturbative calculation. For the first three excited states the extrapolation yields for the triplet of levels E,, with rnj = +$, -i, -3 in the vicinity of EJ”, (fig. 4, right) E t1/2
=
-0.191R
-0.5,
E-I,, = -0.5 ) E-3,2
-+O.l94R
-0.5.
(20)
-1.00
1 m.:-l 1
E
-1.05 C
,
0.2
2
Ii
Fig. 4. The vicinity of the dyon energy levels EF (left) and EAT (right). The energy E is measured in units of the H-atom ground-state energy (13.6 eV), the separation R in units of the Bohr radius a,.
R. Singer et al. / Bound States of an electron
534
These expressions dipole approximation
are in good agreement and first-order
with those
degenerate
given by Drell 13). Using the
perturbation
theory he derived
the
expressions
E flj2
1 =
-
R-QS=-O.l93R-0.5,
2(4-h) E--1,2= -O.SR, E--3,2= +
1
R -0.5= +O.l93R-0.5.
2(4-h)
(21)
Now let us turn to the wave functions. Figs. 5-9 exhibit the unnormalized squared wave functions of the ground state and the first four excited states for R = 0 (Dyon states). The values of the pair of quantum numbers (a, 6) are (0, -1) for the ground state, the second excited and the fourth excited states, (1,0) for the first excited state and (-1, -2) for the third excited state. The ground state is central symmetric with an exponential suppressed r, dependence (fig. 5). The first and the third excited states are symmetric with respect to a coordinate transformation 19+ v - 6 (figs. 6, 8). This corresponds mathematically to an exchange of the sine- and cosine functions in the corresponding expressions for the wave function. The second excited state is again central symmetric but indicates an additional angular dependence (fig. 7). That means that the wave function contains only symmetric sine- and cosine terms. The r, dependence for these three states is equal. For small I, it is dominated by a power factor rl v whereas a=0 b
R=O
z-1
Fig. 5. Unnormalized squared wave function of the ground state for R = 0. In this and all the following wave function figures, two different kinds of representation are displayed: On the left-hand side a side a level three-dimensional representation in the r, , b-plane is exhibited whereas on the right-hand line representation in the upper half plane is shown.
R. Singer et al. / Bound states of an electron a=1
535
A=0
b=O
Fig. 6. Unnormaiized squared wave function of the first excited state for R = 0.
a=
0
R=O
b z-1 1 “i’_Z //------A...
_.
,/H----‘..
Fig. 7. Unno~aiized
squared wave function of the second excited state for R = 0.
a z-1
R=O
b z-2 3 “j’-Z
Fig. 8. Unnormalized squared wave function of the third excited state for R = 0.
536
R. Singer et al. / Bound states of an electron
for large values of rt the wave function is also suppressed exponentially. The maximum one finds at 0.5~1ok ri s O&Z,. The fourth excited state is independent of 6 but exhibits a zero at rr = 2~2, and a first maximum at rI = 4a. (fig. 9). These resuits agree exactly with those mentioned by Malkus 3), who obtained for the squared dyon wave functions /+rj/i(2the expressions
l+?‘*l*~exp(-2r,l%J , sin’ $29cos2 $8 + (3 15/~212~~(7-2P~
\Ji;;‘2/2= ((7 -2pj
8)’ cos’ j@)r~@-r) exp
sin’ $6 cos* $8 + (3 -p)’
(
-2
>
>
sin’ $Y}r:(fitr) exp ---.-2r* (
&a, ) )
l+kI’~‘12cc( I-$J2exp (-2))
(22)
where /.!I= 2 - +‘?. Thus the maximum in the triplet wave function actually occurs at rl = (2 - V~)LX,= 0.586~~ whereas the zero and the first maximum in the fourth excited state are to be found exactly at r, = 2a, and rl = 4a,. In figs. lo-12 we find a phenomenological explanation for the possibility of a strongly bound configuration in the ground state at R = 0.95~2~.These figures show the development of the corresponding wave function with increasing parameter R, We mention that, as the monopole moves out from the origin, the wave function is a=0
R=O
b z-1
Fig. 9. Unnormaljzed squared wave function of the fourth excited state for R =il,
537
R. Singer et al. 1 Bound states of an electron
a =O
R = 0.5
a,
b =-I 1
mj =-Z
Fig. 10. Unnormalized squared wave function of the ground state for increased at the inner regions with r, < O.Sa, whereas at r, b
a
=o
R = 0.95
R = OSa,. 0.50,
The wave function is
it is suppressed.
a0
b z-1 mj =-f
Fig. Il. Unnormalized squared wave function of the ground state for R = 0.95a,. The suppression of the wave function in the vicinity of the monopole and at the outer regions reaches its maximum.
increased at the inner shells (t, < R), whereas in the vicinity of the monopole and at the outer regions ( rI 3 R) the wave function decreases (fig. IO). This compression against the centre causes a shrinkage of the average electron orbital radius which means that the binding energy gets smaller (fig. 4). Taking into account the groundstate electron orbital radius of the neutral H-atom (r, = 1 a,) this effect reaches its maximum at R =50.95uo (fig. ll), whereas for R > 0.95~~ the effect turns over and
R. Singer et al. 1 Bound states of an electron
538
R qlOa,
a=0 b z-1
m.=-l I
2
Fig. 12. Unnormalized squared wave function of the ground state for R = lOa,. The influence of the monopole is very small and the wave function is practically equal to the corresponding one of the neutral H-atom.
the average radius gets larger again, causing the binding energy to rise (fig. 3). If the separation parameter is large enough (R 2 lOa,) the influence of the monopole is nearly negligible and thus the wave function is nearly equal to that of the neutral H-atom in the ground state (fig. 12). The decreasing effect of the wave function in the vicinity of the monopole is observed at any state. For example, we further examine the third excited state. Figs. 8 and 13 show the outward movement of the maximum in the wave function with increasing
R. This is again the reason
for the rising energy of the third excited
0
z-1
R = 0.5 aQ
b z-2
rnj=-3 2
Fig. 13. Unnormalized
squared
wave function
of the third excited
state for R =0.5a,.
state
539
R. Singer et al. / Bound states of an electron R =lOa,
a z-1
b 5-2
Fig. 14. Unnormalized squared wave function of the third excited state for R = lOa,. The wave function is practically equal to that of the first excited state of the neutral H-atom with quantum numbers n = 2, /=l, m=l.
(fig. 4). For R 3 lo&, H-atom with quantum
the state has turned into the first excited numbers n = 2, I= 1, m = 1 (fig. 14).
state of the neutral
5. Conclusion We have presented
a numerical
method
to calculate
in principle
the complete
spectrum of an electron in the field of a magnetic monopole and a proton. We have demonstrated excellent agreement between our results and those of previous works. Two results are of main importance: First, we confirmed the possibility of a strongly bound state configuration with E = -1.135 at R = 0.95a. and, second, we obtained the same results as Drell for the energy levels in the vicinity of R = 0. The present method is completely independent of the coordinate system, but it has the problem that higher eigenvalues are increasingly difficult to obtain. Moreover, for large separation
parameters
certainly in the case of excited close together.
R (R > lOa,), the method
is no more practicable,
states,
levels are situated
because
the energy
very
We would like to express our appreciation to Prof. T.H. Seligman, Dr. P. Zimak and T. Kauer for a number of very helpful discussions. We would also like to thank Prof. J.W. Neuberger, Prof. D.W. Noid and Prof. F. Rose1 for their help with the numerical calculations.
References 1) P.A.M. Dirac, Proc. Roy. 8.0~. London Al33 (1931) 60 2) P.P. Banderet, Helv. Phys. Acta 19 (1946) 503 3) W.V.R. Malkus, Phys. Rev. 83 (1951) 899
540 4)
5) 6) 7) 8) 9) 10)
11) 12) 13) 14) 15) 16) 17) 18) 19) 20) 21) 22) 23) 24)
R. Singer
et
al.
/ Bound
stares
of
an
electron
P. Osland and T.T. Wu, Nucl. Phys. B247 (1984) 421 P. Osland and T.T. Wu, Nucl. Phys. B247 (1984) 450 SK. Bose, J. of Phys. Al8 (1985) 1289 S.K. Bose, J. of Phys. Cl2 (1986) 1135 B. Cabrera. Phys. Rev. Lett. 48 (1982) 1378 F. Kajino, S. Matsuno, Y.K. Yuan and T. Kitamurd, Phys. Rev. Lett. 52 (1984) 1373 T. Hara, M. Honda, Y. Ohno, N. Hayashida, K. Kamata, T. Kifune, G. Tanahashi. M. Mori. Y. Matsubara, M. Tcshima, M. Kobayashi, T. Kondo, K. Nishijima and Y. Totsuka, Phys. Rev. Lett. 56 (1986) 553 DE. Groom, Phys. Reports 140 (1986) 323 G.E. Masek, L.M. Knapp. E. Miller, J.P. Stronski, W. Vernon and J.T. White, Phys. Rev. D35 (1987) 2758 S.D. Drell. N.M. Kroll, M.T. Mueller, S.J. Parke and M.A. Ruderman, Phys. Rev. Lett. 50 (1983) 644 N.M. Kroll, V. Ganapathi, S.J. Parke and S.D. Drell, Proceedings of Monopole 83, ed. J.L. Stone (Plenum, New York, 1984) p. 295 Th.W. Ruijgrok, J.A. Tjon and T.T. Wu, Phys. Lett. B129 (1983) 209 Th.W. Ruijgrok, Acta Phys. Pol. BlS (1984) 305 Y. Kazama, C.N. Yang and A.S. Cioldhaber, Phys. Rev. D15 (1977) 2287 A.S. Goldhaber, Phys. Rev. Dl6 (1977) 1815 C.J. Callias, Phys. Rev. D16 (1977) 3068 T.T. Wu and C.N. Yang, Phys. Rev. Dl2 (1975) 3845 T.T. Wu and C.N. Yang, Nucl. Phys. B107 (1976) 365 J.W. Neuberger and D.W. Noid, Chem. Phys. Lett. 104 (1984) 1 J.H. Wilkinson, The algebraic eigenvalue problem (Clarendon, Oxford, 19653 p. 570 G.D. Smith, Numerical solution of partial differential equations: Finite ditlerence methods (Clarendon, Oxford, 1985) p. 260