- Email: [email protected]
Variational buildup of nuclear shell model bases
ANNALS
OF PHYSICS:
51,
lm-123
Variational
(1969)
Buildup
of Nuclear
Shell Model
Bases
T. SEBE* AND J. NACHAMKIN+ Theoretical Physics Branch, Chalk River Nuclear Laboratories, Chalk
River, Ontario, Canada
A method is proposed for accurately calculating a few of the low-lying states belonging to each nuclear spin within the shell model. A properly chosen initial vector is operated upon iteratively, generating a basis in which the Hamiltonian matrix is tridiagonal. The process may be halted when low-lying eigenvalues, up to a preset number, have been determined to a required accuracy, implying a truncation of the complete vector space. Simple formulae are given for the errors induced by this truncation, both for eigenvalues and eigenfunctions. Illustrative calculations are shown for the Of (T = 0), 2+ (T = 0), 4+ (T = 0), 2+ (T = 1) and 3f (T = 1) states of 20Ne. By taking the leading SUs states as initial vectors, quite accurate low-lying eigenvalues and eigenfunctions are obtained with a small-sized truncation space. The T = 0 states mentioned allow a three dimensional truncation to yield lowest eigenvalues with an error of less than 0.05 MeV as well as corresponding eigenfunctions having an overlap of greater than 99.5 ‘A with the exact ones. Special eigenvalue-problem techniques are discussed in the appendix.
1. INTRODUCTION
Although many complete shell-model calculations have been successfully performed for lp-shell nuclei (I), only limited successhas been met by employing the same methods in the 2s, Id shell. This is mainly due to explosive increases in size of the shell-model basis. Where, especially, more than four particles or four holes are involved most calculations employ a basiswhich has been truncated to include only what is thought to be the most important states of somesubmodel, such as the j-coupling or SU, model (2). Towards the beginning of the shell the SU, model seemsto provide a workable approximation (3) while near the other end, where holes are under consideration, the &coupling model is more often used (4). There seems to be no satisfactory truncation which works well in both regions. This last fact points out the need for exact calculations. * Permanent address: Hosei University, Koganei, Tokyo, Japan. + NRC Post-Doctoral fellow, address after June 1968: Los Alamos Scientific Laboratories, Los Alamos, New Mexico. 100
NUCLEAR
SHELL
MODEL
101
Recently Halbert et al (5) have performed exact shell-model calculations for as many as six particles in the nuclear 2s-Id shell, using forces given by Kuo (6) and Kuo & Brown (7). The dimensions of the shell-model spaces become very large reaching 537 for the 3+ (T = 1) states of 22Ne and 22Na. Although the accomplishment of this feat is encouraging it is still strongly manifested that a refinement of shell-model calculating techniques is required. Reasons for this are abundant and clear as soon as one starts to do large shell-model calculations. Unlike the atomic case, the nuclear residual force has not been very well parameterized so that a great deal of data fitting by repeated diagonalization of the Hamiltonian is inevitable. At best this process is quite time consuming. Then too, in the middle of the 2s-Id shell one encounters bases with definite spin and isospin containing more than 6000 functions (8). In our present state of computer technology this problem cannot be solved efficiently. In each nucleus the number of observed levels having the same isospin T and spin J is small so that only a few of the many shell-model states with the same T and J can be compared with experiment. The method proposed here permits us to investigate the properties of these low-lying states to any degree of accuracy. Our method is based on the Lanczos process of “minimized iterations” (9). This algorithm starts from an initial function or vector $1 and converts the original eigenvalue problem into one of finding the eigensolutions of a tridiagonal matrix. Lanczos applied the method to differential and integral operators, showing that the eigenvalue problem associated with the lowest eigensolutions of these operators may be solved very accurately by diagonalizing small matrices, if the initial function is chosen properly. The main aim of this work is to show that the nuclear Hamiltonian may be reduced to a small tridiagonal matrix with very little resulting error so long as we confine ourselves to studying low-lying states. The basic formulation of the method will be given in Section 2 and an error analysis in Section 3. Numerical calculations were applied to the Of (T = 0), 2+ (T = 0), 4+ (T = 0), 2+ (T = 1) and 3+ (T = 1) states of 20Ne in order to know how effective the method is. The results will be discussed in Section 4.
2. MATHEMATICAL
FORMULATION
The first step in our procedure is to operate with the Hamiltonian H on some unit vector +I (henceforth called the initial vector) to generate another unit vector, c#J~, orthogonal to r& . A new unit vector, &, orthogonal to @I and C& is derived next by the action of H on +2 , and so forth. Once the initial vector is chosen the sequence of operations can usually proceed until the complete space of mutually independent vectors is exhausted. The operations needed to generate
102
SEBE
AND
NACHAMKIN
the sequence of vectors +i in an N-dimensional the following equations H+,
Hilbert
space are described by
= h,+, + K,+, . . . . . . . . .
H9, = Kz--lb-l + h&v + K&c+1 . . . . . . . . . . . . . . . . H%v = &--14+-I
(1)
+ h&v,
where hi = (+i I H I 44,
(24
Ki = <9i+1 I H I +i)-
(2b)
and
It should be noted here that if +1 is orthogonal to any eigenvector of H, then the series of operations described by Eqs. (1) and (2) terminates after m steps, i.e. for some m < N, Km = 0. An important case for which this may happen is when H has some degenerate eigenvalues. In this event a new unit vector may be chosen orthogonal to the +i’s already obtained, and the operations indicated by Eqs. (1) and (2) are continued and/or repeated. We shall discuss this problem in more detail in the next section so, for the time being, it will be assumed that none of the Ki vanishes. We are now interested in solving the matrix equation HJ,=X+
(3)
for a few of the lower values of h. Here, the representation of the Hamiltonian, in the basis +r T& ,..., +N , is the N x N trigiagonal matrix
H=
h, KI 0 . . .
Kl 0 0 hs KS 0 KS h, KS . . . . . . . . .
: : . . . .
: : . . .
: : . . .
* 6-1
H,
: : . &?-I h,
We may now define the truncation of the basis of N vectors to a basis of n vectors by letting the new basis consist of the first n vectors &, &, ,..., +, . The Hamiltonian in our truncated basis is obviously represented by Hen), the n x n matrix in the upper-left corner of H. Henceforth truncation of a basis will have the above meaning.
NUCLEAR
SHELL
MODEL
103
The determinant det / h * 1 - H 1, which yields the characteristic polynomial of H, may be generated in an especially simple manner. Consider the sequence of functions {F’,(h)}, where F_,(X) = 0, F,(h) = 1, and F,(h) = (A - 111)F,-,(h)
- K~~lF&z(x).
(5)
Then it is easy to see by sequentially expanding with respect to the last row of the truncated Hamiltonian Hfn) that det / h . 1 - Htn) 1 = F,(X).
(6)
The functions Fl(A) are discussed again in Appendix I but we may sum up their properties by saying that the sequence {F&X)} is a Sturm Sequence (IO). Two important properties of (F&i)} are: (A)
The nth member of {F&l)) is the characteristic polynomial belonging to H(“), the Hamiltonian in the n-dimensional truncated space.
(B)
Between any two successive of F,(A).
zeros of
F,+,(X) there is exactly one zero
Property (B) follows because FK(X) is a Sturm sequence. Its proof may be found in standard texts, such as that by Wilkinson (10). Properties (A) and (B) are very important because they guarantee that the eigenvalues of the Hamiltonian in our truncated space all have the variational property. That is, we may arbitrarily cut off the sequence of operations defined by Eqs. (1) and (2) after obtaining n vectors +$ , generating Hen). At this point we are assured that the eigenvalues of H (“1, taken in ascending order, are greater than the corresponding exact eigenvalues of the Hamiltonian in the complete space. We may sharpen this result further by denoting the eigenvalues of Hfm), in ascending order, by @), q:“‘,..., v:,“‘. It then follows from (A) and (B) that qt.“’ z > 7)y+l’ > ... > Ai where Ai is the ith eigenvalue of the Hamiltonian, taken in ascending for each i < m the sequence of values Q(m) decreases monotonically and is bounded below by hi .l The possibility definitely exists, several of the lowest r]$ will have converged to a high degree of before the complete Hilbert space is exhausted. 1 This is merely a restatement of the well-known separation theorem (IO).
(7) order. That is, as m increases therefore, that accuracy long
104
SEBE
AND
NACHAMKIN
As a simple example of this phenomenon consider a Hamiltonian H with negative eigenvalues, the lowest of which, h, , is distinct from the others. Choosing an initial vector, +1 , which is not deficient in the ground state of this Hamiltonian we may expand cpl in terms of exact eigenvectors, +z”, and then operate on & with increasing powers of H obtaining
Eventually the term a,,&,“+~ will overwhelm all the others; the value of r for which the normalized right hand side of Eq. (8) is considered a good approximation to (Lr will depend only weakly on the size of the Hilbert space and will mainly depend on the separation of h, from the other eigenvalues. This value of r would, in fact, be simply related to the dimension of the truncation space needed to describe the ground state of H. To see this, it should merely be noticed from Eqs. (1) and (2) that
where the b’s are constants. Thus from Rayleigh’s principle, if an Z-dimensional truncation is required for the ground state then Z < r. By plotting the approximate desired eigenvalues versus the degree of the truncated space the calculation may be halted when, say, the greatest observed root changes by less than some arbitrary small number. Once the iterations have been stopped there still remains the question of just how accurate the approximate solutions are. Because of the simple form of the tridiagonal matrix it is very easy to apply the formulae derived for finding lower bounds on the eigenvalues in the Ritz variational procedure (11). These formulae as well as an analysis of the purity of the truncated eigenfunctions will be given in the next section. The accuracy of our approximation depends somewhat on the choice of an initial vector. We believe that some submodel is sure to give us a good initial vector but we leave further discussion of this point for the following sections. The fact that we may, in principle, solve the energy eigenvalue problem very accurately without a full knowledge of the Hamiltonian matrix is one of the distinct merits of our method. Sometimes, however, it might be more convenient to set up the entire Hamiltonian matrix in some well-known shell-model basis and then perform a transformation to tridiagonal form. In such cases there are two other familiar algorithms which will bring about this transformation, those of Givens and of Householder (10). A comparison of the Lanczos, Givens and Householder methods is made in Appendix II.
NUCLEAR
SHELL
3. ACCURACY
OF
MODEL
THE
105
METHOD
One of the major defects of our truncation procedure is that there is the possibility of missing a low-lying level. Although the probability of this occurring in a realistic problem is very small, it is easy to see how it may come about. Consider the initial vector, +I , and expand it in terms of the exact (unknown) eigenvectors of H, (Ly, where H+z” = hi+:“:
Repeated operation of H on +I yields WT) = Hr+,
= g /l,ra,+;x i=l
so that no new components appear in each successive @MT).Since the new basis functions, the +‘s, are linear combinations of the @fr)‘s, if +1 were orthogonal to a low-lying solution +,8” then our method would not yield a state which corresponds to $2;“. If a degeneracy exists, then for some i and j Xi = Xi and the combination (9) will always appear in each Cptr). Thus +rij will act as though it were a single vector and another linear combination of +,i”” and +j”” will be required. These defects can easily be overcome, however, by repeating the calculation using another initial vector and comparing the results with the original. Since the initial vectors are chosen from physical considerations it is unlikely that both initial vectors would be deficient in some low-lying state. In the case of a degeneracy the new initial vector will, in general, contain a different linear combination of the degenerate-state vectors, +ij, = bi +F + bj Jr,‘“. Solving with the new initial vector will then yield another independent linear combination of the degenerate vectors. If higher degeneracies are expected this process can be repeated, as in the example given in the next section. It should be realized that the diagonalization of the Hamiltonian in a truncated space defines a variational calculation. Formulae have been established for such calculations which yield upper and lower bounds on the eigenvalues as well as an estimate of the accuracy of the eigenvectors. As we have mentioned before, when there are degeneracies our truncated basis is divided into subspaces containing no degeneracies. The Hamiltonian may, in fact, be diagonalized separately in each of these subspaces. We will now assume
106
SEBE AND NACHAMKIN
that none of the eigenvalues we seek is degenerate.2 It is convenient to define what turns out to be an important concept. Denote the exact eigenvalues in the full space of dimension N by Xi and approximate ones in the truncated space of dimension I by Q, where
Now if
for i < A4 and M < Z then the sequence of eigenvalues h, , A2 ,..., & will be called separated with respect to the truncation to the space of dimension I. For separated eigenvalues we can make use of the analysis given by T. Kato (II), which reduces to a simple form in our case. Let the Z(N) approximate (exact) eigenvectors in the truncated (complete) space be denoted by ai( where (oJ+:“) belongs to an approximate (exact) eigenvalue Q(&). The quantity Q is an upper bound for Xi as mentioned in Section 2. A lower bound for hi is given in terms of the quantity Q (30) where ci2 = (4$ 1(H - qJ2 I mi)
(124
= (Qi 1H2 1Qi) - vi2.
(12b)
To obtain ci2 we first expand Qi in terms of the vectors & defined in Eq. (1). (13) k=l
The Hamiltonian can be expressed as the sum of a part acting only in the truncated space and a part which takes only the last vector, +r, out of this space. More explicitly, the first part is diagonal in the truncated space and yields viQIi when H is applied to mi while the second part transforms +r into K,+,+, , where the quantity K1 is the off-diagonal element defined in Eq. (2b). * Wilkinson (10) considers tridiagonal matrices where some of the eigenvalues are pathologically close to each other even though none of the off-diagonal matrix elements are small with respect to the diagonal elements. This situation, when it obtains, is always detectable because the Sturm sequence property (see Appendix II) allows us to count the roots lying between two numbers, no matter how close the roots are to each other. A simple remedy, should this happen, is to choose a new initial vector and re-tridiagonalize H. Almost invariably a small off-diagonal matrix element will appear, which may be set equal to zero with small resulting error.
NUCLEAR
SHELL
107
MODEL
Thus
and therefore (Qi 1H” 1 *J
= 7; + 1uIiKI j2.
(15)
From Eqs. (12) and (15) we obtain ei = j ariKr I.
(16)
If Eq. (11) holds for vi we may introduce the new quantity relation -qi < @h< hit,.
(Y~which satisfies the (17)
Kato shows now that (II) (18) The most precise (or greatest) lower bound is obtained by replacing oli by its maximum value, h,+l . When ci is very small it becomes permissible to take (~l~= qi+l instead of X,+l . In this case
Ai 3 xi = .qi -
Q2 . %+1- qi
(19)
If it happens that
then Ai will be called well separated with respect to the truncation hi 3 Yi = 7)i -
Ei s
and (21)
The quantity Y, gives a rather crude lower bound since Yi < Xi when ci is small. It should be noted that if a required eigenvalue is not well separated with respect to the truncation the quantity Yi as well as Xi may not yield a true lower bound for it. This may be seen in Figures l-5 where the lower bounds Yi are plotted as dotted curves. For small truncation spaces some of the apparent lower bounds for the higher eigenvalues are at first greater than the exact eigenvalues but after the eigenvalues become well separated the calculated lower bounds become true lower bounds and start to increase. The quantity qi is always an upper bound for Ai . Variations of Yi as the degree of the truncation space increases are discussed more precisely in Appendix III.
108
SEBE AND
NACHAMKIN
A measure of the accuracy of the eigenvectors in the truncated space is given by the quantity cio, where
The formula for the upper limit of l io as given by Kato (11), when the errors in approximate eigenvalues are sufficiently small, reduces to
where xi is the smaller of the two quantities (~+r - Q) and (Q - qiml). In conclusion, for a truncation to the space containing I vectors a complete variational error analysis may be performed with the knowledge of KI .
4. APPLICATION
TO
THE
NUCLEAR
SHELL
MODEL
Although our main interest lies in developing a method for attacking largerdimensional shell-model problems, it is perhaps best at this point to illustrate the main features of the present method with some simple examples. For these we have taken the 0+ (isobaric spin T = 0), 2+ (T = 0), 4+ (T = 0), 2f (T = 1) and 3+ (T = 1) states in 20Ne. For the first three sets the low-lying level spacings are fairly large while for the last two they are rather small. The dimensions of the untruncated 2~Id shell-model spaces for the above states are 21, 56, 44, 66 and 69 respectively. The nuclear Hamiltonian is taken to be (24)
where the first term is the single-particle spin-orbit interaction, the second expresses the splitting of the single-particle 2s and Id levels, and the last term is the residual two-body interaction. Values for [ and 77 are chosen so as to reproduce the observed energies in 17F and 170 (22), e.g. 5 = -2.03 MeV and 7 = - 1.27 MeV. The interaction Vij is assumed to be composed only of central forces having a gaussian radial form with range parameter equal to 0.7 (23), multiplied by an exchange mixture determined by p
zz.z0 ,
V13 = -60
MeV 9
V31 = -45 MeV ,
VW = 22.5 MeV.
(25)
This Hamiltonian was used by Akiyama et al (14) and yields results which are in reasonable agreement with observations in the 2s-ld shell.
NUCLEAR SHELL MODEL
109
Figs. 1,2 and 3 show the variations of the approximate eigenvaluesvi (solid lines) and the lower bounds, Yi , calculated from Eq. (21) (dashedlines), asthe dimensions of the truncated spaces increase. The graphs correspond to the four lowest 0+ (T = 0), 2+ (T = 0) and 4+ (T = 0) levels. For the sake of comparison the results of a full shell-model calculation and two other methods of truncation, to be explained later, are labelled (A), (B) and (C), respectively, and included in each figure. Using the leading SU, states 1[4] (80) J = L) as initial vectors, surprisingly rapid convergence is noted for the lowest level (first state) belonging to each spin; three-dimensional spaces allow us to obtain eigenvalues with an error of less than 0.1 MeV, and eigenfunctions whose overlaps with the exact ones
2 4 6 6 IO lZ;l4 I6 I8 m, n (DIMENSION OF TRUNCATION
SPACE)
FIG. 1. Approximate eigenvalues vi (solid lines) and calculated bounds Yj (dotted lines) for the first four O+ (T = 0) levels in 20Ne, plotted against increasing dimension of the truncation space. Levels labeled (A) are the exact eigenvalues hi . The sets (B) and (C) correspond to truncations arising from the supermultiplet and &coupling models, as explained in the text.
2
4
6
8
IO I2 I4 16 I8 20 n
FIG. 2. Same as Fig. 1 for 2+ (T = 0) levels.
110
SEBE
AND
NACHAMKIN
(T=O)
4+
711 5 kY,
I
,
f
I
2
4
6
8
FIG.
;-
+--J
3.
I
I
,
I
I
,
(A)
(8)
6)
IO 12 14 16 18 20
n
Same as Fig.
1 for 4+ (T = 0) levels.
exceed 99.5 %. For the second states the above accuracy was reached at n = 7 for O+, II = 11 for 2+ and n = 8 for 4+. Two more basis functions are necessary in each case to determine the third states to the same accuracy. As far as the first three levels with these spins are concerned, convergence does not seem to be strongly dependent upon the size of the complete shell-model space. The fourth 2+ level, however, converges more slowly than the fourth O+ or 4+ states, probably because the spacing between the fourth and fifth 2+ levels is smaller. In general the values for Yi do not start out as true lower bounds but show a decreasing trend as n increases. Eventually Yi becomes a true lower bound for the ith level and begins to increase as n increases, finally coming very close to its corresponding eigenvalue ~7~. In practise, then, we determine convergence by investigating the Yi’s as well as the 7;s. In order to display and clarify some of the characteristics of our method of solution, the values of ci in Eq. (16) are given in Table I for the lowest two Of (T = 0), 2+ (T = 0) and 4+ (T = 0) states for 12= 5, 10, 15 and 20. TABLE
I
VALUES OF THE QUANTITY ci, IN MEV, AS GIVEN BY Eqs. (12a) AND (16), FOR T = 0 STATES IN 20Ne. THE TRUNCATION-SPACE DIMBNSION IS DENOTED BY n
n
5 10 15 20
01+ 2.0 3.5 3.4 9.3
x x x x
10-l lo-’ 10-a 10-l’
a+
?A+ 2.6 3.1 x 10-Z 1.4 x 10-b 1.5 x 10-10
1.5 8.5 1.7 3.5
x x x x
10-l 10-a 10-e 10-e
A+ 3.0 4.5 x 10-l 8.9 x 10-a 2.3 x 1O-4
41+ 1.5 6.1 8.7 3.8
x x x x
10-l lo-* lo-’ lo-lo
4s+ 1.9 1.5 x 10-l 3.4 x 10-a 1.5 x 10-S
NUCLEAR
SHELL
MODEL
111
As already mentioned, the quantities ei give us measures of the accuracy of the approximate solutions. For the ten-dimensional truncation, that is n = 10, values of l i for the second 2+ and 4+ states have fallen to 0.45 MeV and 0.15 MeV, respectively. Using these quantities in Eqs. (18) and (23) one finds that for the second 2+ (T = 0) state: 9.619 MeV < X, < 9.730 MeV, (CD21i/y)
> 0.9725
(264 (26b)
and for the second 4+ (T = 0) state: 11.963 MeV < h, < 11.972 MeV, @j2 I t&> > 0.9976.
(274 W’b)
These values should be compared respectively with the exact values; 9.715 MeV, 0.9834 for 2+ (T = 0) and 11.970 MeV, 0.9999 for 4’ (T = 0). Thus we may conclude that the more accurate formula, Eq. (18), is useful for sufficiently small values of ei . (In the practical examples above, .45 MeV and .15 MeV proved to be sufficiently small.) Let us now return to the other methods of truncation referred to in Figs. 1, 2, and 3. Levels are labelled (B) and (C) corresponding to Wigner supermultiplet theory (1.5) and jj coupling, respectively. If only the basis functions belonging to the maximum symmetry [4] states were used in supermultiplet theory we would incur fairly large errors in the wave functions, exceeding 6 y0 even for the ground state of 20Ne. For this reason we follow the procedure of Inoue et al (3) and include also basisfunctions belonging to [3, I] symmetry. Functions with this symmetry interact with those of [4] symmetry through the spin-orbit interaction. Now the dimensions of the Of (T = 0), 2+ (T = 0) and 4+ (T = 0) spacesare 8, 22 and 18, respectively. This improvement of the supermultiplet model yields very good results for the lowest two states of each nuclear spin but, as seenin Figs. 1, 2 and 3, it begins to break down for the third states. Comparing with the exact wavefunctions the third O+, 2f and 4+ levels have overlaps of 0.992, 0.987 and 0.759, in that order. The fourth levels, however, cannot be trusted at all in this approximation. The truncated jj-coupling model we considered embraces only configurations of the form (I&$ (2s& (l&)7, where r < 1. In other words no configuration has more than one particle in the Id,,, orbit. Dimensions are now 10for 0+ (T = 0), 29 for 2+ (T = 0) and 25 for 4+ (T = 0). From the figures it is apparent that this approximation is rather a poor one compared with the supermultiplet theory;
-~__.
112
SJ3BE AND
NACHAMKIN
2+(T=
2-4682
II
14 16 18 20
n
FIG. 4. Same as Fig. 1 for 2+ (7’ = 1) levels. Greek letters AI and AZare omitted for clarity. 3+ (T= I 1
2
4
6
8
1
IO I2 I4 I6 I8 20 ”
FIG. 5. Same as Fig. 1 for 3+ (T = 1) levels. Most of the dotted Y, curve has been omitted for clarity. Also, for reasons of clarity, greek letters ha and h, have been omitted.
e.g. the overlaps of the lowest eigenfunctions with the exact ones are only 0.878 for O+,0.911 for 2+ and 0.918 for 4+. The eigenvalues rli and lower bounds Yi for 2f (T = 1) and 3+ (T = 1) states in z”Ne are given in Figs. 4 and 5, with the same notation as in Figs. l-3. An interesting phenomenon is displayed in Figs. 4 and 5. The lowest 3+ level converges faster than the lowest 2+ level even though the complete space for 3+ functions is larger than the complete 2+ space. Initial vectors were chosen here to correspond to the states 1[31] S = 0 J = L), simply for convenience in computing. The convergencies are rather slow as compared to the T = 0 cases.This probably arises, in part, becausethe spacings between adjacent levels are smaller in the T = 1 case. The lack of a proper initial vector is, however, also a cause of slow convergence.
NUCLEAR
SHELL
113
MODEL
As an example of the latter cause consider the four states belonging to the representation [3,1] (6,1), i.e. / S = 0, L = J), I S = I, L = J - I), 1S = I, L = J), j S = 1, L = J + 1). The expectation value of the Hamiltonian is lowest in the state with S = 1, L = J - 1. Taking an initial vector corresponding to this state, instead of the one actually used in Figs. 4 and 5, allows an accurate truncation with as few as seven vectors to describe both the 2f and 3+ lowest states. If all four lowest levels are required for both spins then it turns out that there is no remarkable difference between the two different initial vectors. Truncating at n = 20 yields overlaps between the exact and approximate eigenfunctions of greater than 99.5 ‘A for the levels shown in Figs. 4 and 5. The values of l i for T = 1 are tabulated in Table II and show slower decreasesagainst increasing n than do those for T = 0. The supermultiplet and jj models are again labelled (B) and (C), respectively, in Figs. 4 and 5. The only difference here is that only states with symmetry [3, I] are included in the supermultiplet scheme.We did this becausestatesof symmetries [2,2] and [2, 1, l] may couple directly with [3, I] through spin-orbit and two-body forces, making it difficult to obtain a truncation of proper size for a comparison with our calculations. In this model both the 2+ and 3+ spacesare 24-dimensional. Employing this supermultiplet model the lowest 2f state has an overlap of 0.951 with the second lowest of the exact calculation while the second lowest has an overlap of 0.921 with the lowest exact state. This means that the approximation predicted an inverted order for the lowest two levels. For 3+, the second lowest eigenfunction has only a 0.776 overlap with the exact one. Taking the samejj configurations as for the T = 0 casethe dimensions are 36 for the 2+ spaceand 38 for the 3+. In spite of the large spacesin this approximation, it does not yield very good results; overlaps with the exact eigenfunctions do not exceed 0.97 for any of the 2f or 3’ levels shown. TABLE SAME
2,+
n
5
AS TABLE
I FOR
2,+
II T =
29+
1 STATES
3,+
IN
20Ne
3,+
3,+
7.9 x 10-l
1.9
2.6
10
5.8
1.5
1.2
1.2 7.4 x 10-l
2.9 8.3 x 10-l
4.0 3.9 8.5 x 10-l
x IO-’
15
2.3 x IO-’
2.6 x 10-l
1.1
7.9 x 10-z
2.6 x 10-l
20
1.1
x
1.7 x 10-S
1.7 x 10-l
3.7 x 10-a
4.0 x 10-z
2.2 x IO-1
25
1.1
x 10-a
2.1 x 10-a
3.9 x 10-a
4.1 x 10-z
9.9 x 10-S
8.4 x 1O-s
10-Z
114
SEBE AND
NACHAMKIN
We would again like to point out that, by using our method, the size of the truncated space needed to obtain a good approximation depends mainly on the level structure of the states of interest and only weakly on the full size of the shell-model basis, while other methods of truncation yield bases which rapidly increase in size as particle number increases. As a further example we consider a case wherein a degeneracy arises, by taking our Hamiltonian to consist merely of the single-particle spin-orbit force. Four different initial vectors were chosen for 2+ (7’ = 0) of 20Ne, namely 1 [4] (8, 0) S = 0 L = 2), 1 [4] (0, 4) s = 0 L = 2), j [3, l] (4, 2) s = 1 L = 3), and I [3, l] (6, 1) S = 1 L = 3). Lanczos iterations were stopped in each case when an off-diagonal matrix element became very small, and the resulting tridiagonal matrix was solved for its lowest eigenstate. In all cases this state turned out to have the same energy as the configuration (&j4. Calling the four states so obtained $1 , I&, I,& and z,& the 4 x 4 symmetric matrix with elements Mij = (J,& j 4,) was constructed: 0.8723 1 -0.7658 Quoting
0.1531 -0.0640 1 i -0.3451
I.
(28)
results to four decimal places detM} = 0.0000, det(M(4 14)) = 0.0788 # 0,
(294 (29b)
where M(4 x 4) is the same as M but with the fourth row and fourth column deleted. Eq. (29a) is the criterion that +I ,..., $4 are linearly dependent while Eq. (29b) tells us that I+&& and z,& are linearly independent. Deleting other row-column combinations also gives a non-vanishing determinant. This simple but important example shows how a degeneracy, in this case of third order, may be resolved by assuming new initial vectors. (Even though it is highly improbable that a degeneracy would not be resolved in the way described above we feel that the method of choosing new initial vectors could well become an art based on intuition.) 5. CONCLUDING
REMARKS
As has already been shown our method has several characteristics which make it suitable for solving for the lower-lying states in large shell-model calculations.
NUCLEAR
SHELL
115
MODEL
Restated, they are: (1)
The special iterational character, which suggests that the size of the truncated space will be almost independent of the size of the original shell-model space;
(2)
A very simple set of error formulae, safely accurate truncation;
(3)
The time consuming eliminated.
diagonalization
which
allows
an economical
of large matrices
yet
is virtually
Another question is what is the best method of constructing the tridiagonal Hamiltonian matrix directly. As mentioned before, it is not really necessary to generate the entire matrix in the complete shell-model basis. In principle the I x I tridiagonal matrix of H may be expressed in terms of the expectation values of H, HZ,..., Hz’ evaluated in the initial state. T.D. Newton attacked the problem of generating the tridiagonal matrix directly, but since he was attempting full calculations using a Slater-determinant basis, the large problems were out of his reach (26). In all the examples we have given, the entire Hamiltonian was set up and then tridiagonalized. We feel that this may not be the best way to advance to larger-dimensional problems. At present work is proceeding on problems in the 2s-ld shell which involve more than four particles or holes. Finally but perhaps most importantly, we have seen that the properties of many low-lying states may be grasped by using very small truncation spaces. Since convergence in these spaces depends, to a great extent, on the submodel from which the initial vectors are chosen much can be added to the physical interpretation of these states. Choosing an initial vector from, say, the SU, model, the jj-coupling model, the Nilsson model (17) or the model proposed by Harvey & Sebe (18), and watching for the rapid convergence of the resulting spaces is a fairly stringent test of these models’ efficacy. In short our method supplies a testing ground for the various nuclear models.
The authors would like to express their thanks to Dr. M. Harvey for many stimulating discussions and to Prof. T. D. Newton for some incisive remarks. We are also indebted to Mr. H. H. Clayton for his kind interest in our work. One of us (T. S.) wishes to thank Mr. G. C. Hanna for the hospitality extended to him at the Chalk River Nuclear Laboratories.
116
SEBE AND NACHAMKIN
APPENDIX
I.
EIGENSOLUTIONS
OF TRIDIAGONAL
MATRICES
In this appendix we consider the problem of numerically computing the eigenvalues, hi, and the eigenvectors, Gi, of a real, symmetric, N x N tridiagonal matrix. We will be making use of the following property of Sturm Sequences (see Section 2 for properties (A) and (B)): (C) For any pO the number of times the Sturm Sequence {F&J) changes sign is equal to the number of eigenvalues of H which are greater than p,, . Starting from initial upper and lower bounds each hi is derived by the bisection method (10). We feel that the method of bisections is of enough interest to be explained in detail. The following steps define an algorithm which will isolate the eigenvalues of a tridiagonal matrix quickly and with high accuracy.
(1) Set initial upper and lower bounds on the eigenvalue desired, defining an interval. (2) Compute the Sturm sequence, using Eq. (5), at the midpoint of this interval, counting the number of times it changes sign. (3) From the results of step (2) the eigenvalue is determined to be in either the upper or the lower half interval, so that a new interval is defined which is half the size of the current interval. (4) If the new interval is smaller than some predetermined size then the sought-after eigenvalue may be set equal to the number corresponding to the middle of this interval otherwise go back to step (2). Step (1) involves finding upper and lower bounds on given eigenvalues without solving for them. In many physical problems this might be fairly easy but we make use of a simple infallible method to obtain reasonable upper and lower bounds for each eigenvalue of a tridiagonal matrix. To do this we make use of Gershgorin’s theorems which state that (IO): I) For a matrix H with elements Hii , each interval (or disk in the complex plane, if H has complex eigenvalues) defined by I Hii - P I < Gi ,
U-1)
where Gi = c I HijI, j=i contains at least one eigenvalue of H, and II) Each eigenvalue of H lies in at least one Gershgorin disk.
U-2)
NUCLEAR
SHELL
117
MODEL
If H is tridiagonal Gi never involves summing more than two terms. For ease of computation in the real, symmetric tridiagonal case we choose the off-diagonal matrix element of largest absolute size, calling it bm,, . Next, the diagonal elements are reordered according to algebraic size so that the nth one, a, , determines the interval an - 2bmax < p < a, + 2bmax,
(1-3)
and Gershgorin’s theorem guarantees our finding the nth eigenvalue in this interval. As the calculation for the eigenvalues proceeds (Steps (l)-(4) above), new bounds on subsequent eigenvalues may be obtained if the so-defined intervals overlap. Once the desired eigenvalues of H are evaluated we next seek the corresponding eigenvectors. If we know some eigenvalue, Xi, exactly and compute the Sturm sequence {F&J} exactly then it is easy to verify that the vector Gi, with components 01~, 01~,..., 01~, where =
FN-l(xi),
OIN-1
=
-HN.N--1FN--2@i)v
OIiV-7:
=
t-1”
OIN
(1-4)
HN,N-lHN-l.N-2
*‘*
HN-k+l,N-kFN-k-l(xi),
satisfies the matrix equation (Ai + 1 - H) +,i = 0
(1-5)
In the general case we cannot obtain the exact value Xi but must be satisfied with an approximation qi , for which there is no nontrivial solution of
We may, however,
(Q * 1 - H) tj+ = 0.
(I-5a)
(pi . 1 - H) +i = Xi
V-6)
solve
where xi is nonnull, obtaining the approximate
eigenvector
& = (Q . 1 - H)-l
xi
& , where (I-7)
If we now take xi = & and solve Eq. (I-6) again the new +i will be an even better approximation to the exact eigenvector +,i, if Q is close to hi. To see this, expand xi in terms of the exact eigenvectors, Xi =
C a&
W)
118
SEBE AND NACHAMKIN
Then (Q * 1 - H)-7 xc = c a&/(~
* 1 - Ai)
(I-9)
and Gs will dominate the series if, for i f j, I qi - hi I Q I r)e - h I.
(I-10)
Repeating the process purifies the eigenvector very quickly if Eq. (I-10) holds and xi is not deficient in +i . Following Wilkinson fairly closely the purification process described above is performed four times (Wilkinson recommends three), solving the linear equations by elimination and back substitution (IO). For approximate eigenvalues correct to four or five significant digits this process yielded eigenvectors, for all the physical cases tested, which were so accurate that the eigenvalues could be recalculated to ten or twelve significant digits. APPENDIX
II.
RELATIVE MERITS OF TRIDIAGONALIZATION
PROCESSES
Since the Hamiltonian is quite conveniently expressed when in tridiagonal form we feel that more should be said about reducing it to this form. The three methods usually employed for the transformation are the Givens (29), Householder (19), and Lanczos algorithms (9). The first two require the Hamiltonian to be in matrix form from the outset, that is, a “parent” matrix must be known. Theoretically, the last method is capable of generating the tridiagonal “daughter” matrix directly by knowing the expectation value of powers of the Hamiltonian in an initial state. Most frequently the parent matrix is generated. In realistic shell-model problems this is usually a large sparse matrix whose low-lying eigensolutions are required. There seems to be no doubt that if the parent matrix is specified, the Householder method is superior to both the Givens and Lanczos methods in achieving total tridiagonalization, for reasons of storage and speed. Accuracy of all three is high enough to be of minor consideration (10). The Givens method proceeds via two-dimensional rotations, the (i, j)-th element being eliminated by a rotation in the (i + 1, j)-th plane. An additional matrix, the same size as the parent matrix, is needed just to store the details of the transformation. Thus, even if the N x N parent matrix is real and symmetric, a total of N(N - 1) + 1 storage locations is needed to completely tridiagonalize it and retain all the information concerning the transformation. At each point of a Lanczos matrix iteration a new N-dimensional vector must be stored so that a total of N(N + 2) - 1 storage locations must be allocated for total tridiagonalization. Since Householder’s tridiagonalization requires only $N(N + 3) locations, it can accommodate much larger matrices.
NUCLEAR
SHELL
119
MODEL
As far as speed is concerned, the Householder & Givens processes require respectively approximately $V3 and 4N3 multiplications as well as N and $(N - l)(N - 2) square roots for tridiagonalization, while the Lanczos method requires a reorthogonalization process at each step and thus requires a total of approximately as many multiplications and square roots as the Givens method. After the tridiagonal matrix has been solved for eigenvectors the final eigenvectors of the parent matrix are found by the Givens method at only half the speed of either the Householder or Lanczos method. Even when an exotic unit vector + is to be used as an initial vector in the Lanczos algorithm, it is far more economical to first transform H by the unitary orthogonal matrix P where P = 1 - (+ - eJ(+ - elY/[l
- (e,, +)I,
(II-l)
where XT is the transpose of X, and H’ = P-lHP
= PHP.
(11-2)
Each step of the Householder method applied to II’ is equivalent to the corresponding Lanczos iteration. It should be noted that P was chosen in that form because P+ = e, (11-3) and a Householder routine may easily be adapted to perform this transformation. Only N extra storage locations are required if this recommendation is followed. If the tridiagonal representation of the Hamiltonian is to be generated directly by means of Lanczos iterations the exciting possibility exists that, given an appropriate initial vector, the number of iterations needed for very accurate determination of the low-lying states may be very small. In the future if a variational truncation based on the Lanczos procedure is not used, there may be no alternative to diagonalizing unreasonably large matrices in shell-model calculations. We have already mentioned that the actual numerical errors encountered by tridiagonalizing a matrix by the methods mentioned were very small. To get an idea of the errors introduced we will consider the percentage mean-square deviation, x, (11-4) X = Cithi l)J2/Cihi2~ where hi and qi are the exact eigenvalues of the parent matrix and the corresponding exact eigenvalues of the daughter tridhgonal matrix. Since the computer performing the algorithm does not usually work with infinite accuracy, x f 0.
120
SEBE AND NACHAMKIN
Wilkinson gives what turn out to be overly large upper bounds for x in the Givens & Householder reductions (10). They are, respectively, x(Givens) < 12 x Naj2 x 231 + 6 x 2-t)4N-7, x(Householder) < 40 x N x 2-31 + 20 x 2-t)2N,
(II-5a) (II-5b)
where t is the binary bit size of the numbers used by the computer. For the Chalk River G-20 computer, 2-t M lo-l2 for double precision arithmetic. Assuming N = lo4 (a very large matrix) then x(Givens) < 1.2 x 10-5, x(Householder) < 4 x 10-7.
(II-6a) (II-6b)
Such a matrix would take us more than 5 x lo4 hours just to reduce to tridiagonal form! Just as a comparison, the six-sweep Jacobi method has the bound (IO) x(Jacobi) < 108 x 2-t x N312(1 + 9 x 2-t)6(4N-7) = 1.08 x 10-4 (11-7) We have not made a detailed error analysis of the Lanczos process with reorthogonalizations but, since approximately the same kinds and number of steps are used as by the Givens process, the error should correspond closely to x(Givens). Thus, as far as numerical errors involved in calculating the eigenvalues are concerned, any shell-model matrix which could be fed into the memory of almost any electronic computer today could be solved for low-lying eigenvalues with high accuracy. (In fact statistical cancelling of rounding errors will improve the bounds in Eqs. (11-6) and (11-7) by several orders of magnitude.)
APPENDIX
III.
BEHAVIOR
OF THE LOWER BOUND
Yi
The method we have proposed has the property, common to most iteration processes, that the lower the eigenvalue the faster its approximation converges. We shall therefore assume that we have achieved an I-dimensional truncation space wherein the first k - 1 solutions have converged to within a preassigned accuracy and the kth approximate eigenvalue &) is fairly close to the exact value hk . Under these assumptions the kth approximate eigenvector 0:‘) may be expanded very accurately in terms of exact eigenvectors +;“, with eigenvalues hl > hl, , i.e. (III-l)
NUCLEAR
SHELL
121
MODEL
(111-2)
and the quantity
eka in Eq. (12.b) is expressed as (111-3)
If we assume that the curve of Yk (=rii) - l 3 plotted against Zdefines a reasonably smooth function Y,(Z), then Eqs. (111-1-111-3) may be used to obtain the tangent to this curve as follows:
dYk ___
(III-4a)
dZ
(III-4b)
t qy ,=f&[F-h,? + 2(g) + ‘,
We now make the assumption that all the h’s are negative and that we are in the region where A, and A,,, are well separated, i.e. (111-5) Then, obviously [-A,a
+ 2(7&f’ + Ek) A,] < [-A;,,
Replacing all the Al in Eq. (III-4b) 1 --dY, dZ > 2Ek
+ 2(4? + 4 &+,I, by A,+1 it follows
[-X,2
1 > k + 1. (111-6)
that
+ 2(?$’ + Ek) x,1 y
+ [ -h;+l +
37:’
+ Ek)~,+,IL=t+lqgq
2
(111-7)
since we are assuming that
d(Q) -x-
for
1 > k.
(111-8)
122
SEBE AND NACHAMKIN
Then (111-9) The right-hand side of Eq. (111-9) will now be positive if 72’ + Ek < (A, + 4+,)/2,
(III-1Oa)
$j’ + Ek > (A, + &+JP.
(III-lob)
and will be negative if
In other words the quantity Yk may, on the whole, be expected to be a decreasing function of Z roughly until the point is reached where Eq. (III-1Oa) is true. After that Yk increases monotonically with increasing I.
APPENDIX
IV.
A COROLLARY OF GERSHGORIN'S FIRST THEOREM
In Appendix 3 mention was made of the Gershgorin disk theorem wherein at least one eigenvalue of a matrix H may be found in a disk centered at Hii and of radius C:‘i 1Hij I. This theorem may be improved upon slightly in that the radius may be taken to be Gj = cc’,+” 1Hij 12)1/2.
To show this, consider the ith first row and column so that the tively while the eigenvalues of Householder procedure produces
row. Exchange the ith row and column with the elements Hij and H,, become Hij and Hj; respecthe matrix are unchanged. The first step in a a matrix H” where (IV-2a)
Hii = %I , H,“,
(IV-l)
= ztCZj>~ I HI, W2
(IV-2b)
and H& = 0,
(n > 2).
(IV-2c)
Applying Gershgorin’s theorem to H”, which has the same eigenvalues as H, immediately yields an eigenvalue in the disk centered at H,, with radius given by Eq. (IV-l). It should be noted that these new disks do not necessarily satisfy the second Gershgorin theorem. RECEIVED:
June 27, 1968
NUCLEAR
SHELL
123
MODEL
REFERENCES
1. For example: D. R. INGLIS, Rec. Mod. Whys. 25, 390 (1953); D. KURATH, Phys. Reu. 101, 216 (1956); HILL, AND
S. COHEN AND D. KURATH, Nucl. Whys. 73, 1 (1965); J. NACHAMKIN, Nucl. Phus. A106, 62 (1968).
P. GOLDHAMMER,
J. R.
Proc. Roy. Sot. A245, 128, 562 (1958); J. P. ELLIOTT AND M. HARVEY, Proc. Roy. Sot. A272, 551 (1963). 3. For example: M. BARNERJEE AND C. LEVINSON, Phys. REV. 130, 1036 (1963); T. INOUE, T. SEBE, H. HAGIWARA AND A. ARIMA, Nlrcl. Whys. 59, 1 (1964); M. HARVEY, The nuclear SU, model. Advances Nucl. Phys. 1, 67 (1968). 4. For example: P. W. M. GLAUDEMANS, G. WUCHERS, AND P. J. BRUSSARD, N~cl. Phys. 56, 2. J. P. ELLIOTT,
529, 548 (1964). 5. E. C. HALBERT,
6.
J. B. MCGRORY,
T. T. S. Kuo, N~cl. Phys. A103,
AND
B. H.
WILDENTHAL,
AND G. E. BROWN, Nlrcl. Phys. T. SEBE AND M. HARVEY, Chalk River Report C. LANCZOS, J. Res. Nat. Bur. Stand. 45, 255 J. H. WILKINSON, “The Algebraic Eigenvalue T. KATO, J. Phys. Sac. Japan 4, 434 ( 1949).
7. T. T. S. Kuo
8. 9. 10. II. 12. 13. 14. 15.
F. AIZENBERG-SELOVE R. THIEBERGER, Y. AKIYAMA, E. P. WIGNER,
AND
Phys. Rev. Letf. 20,
1112
(1968).
71 (1967).
T. LAURITZEN,
Nuc/.
85, 40 (1966). AECL-3007 (1968). (1950). Problem.” Clarendon Press, Oxford, Phys. 11,
1965.
1 (1959).
Nuc/. Phys. 2, 533 (1957). A. ARIMA, AND T. SEBE, to
be published. Phys. Rev. 62, 438 (1942). 16. T. D. NEWTON, Chalk River Report, AECL-2223 (1965). 17. S. G. NIUSON, Dan. Mat. Fys. Medd. 29, 75 (1956). 18. M. HARVEY AND T. SEBE, to be published. 19. For example: W. GIVENS, Oak Ridge National Laboratory Report ORNL-1574; HOUSEHOLDER, J. Ass. Camp. Mach. 5, 335, 339 (1958); P. A. Wm AND R. R. Mathematics of Computation 18, 457 (1964).
A. S. BROWN,
OF PHYSICS:
51,
lm-123
Variational
(1969)
Buildup
of Nuclear
Shell Model
Bases
T. SEBE* AND J. NACHAMKIN+ Theoretical Physics Branch, Chalk River Nuclear Laboratories, Chalk
River, Ontario, Canada
A method is proposed for accurately calculating a few of the low-lying states belonging to each nuclear spin within the shell model. A properly chosen initial vector is operated upon iteratively, generating a basis in which the Hamiltonian matrix is tridiagonal. The process may be halted when low-lying eigenvalues, up to a preset number, have been determined to a required accuracy, implying a truncation of the complete vector space. Simple formulae are given for the errors induced by this truncation, both for eigenvalues and eigenfunctions. Illustrative calculations are shown for the Of (T = 0), 2+ (T = 0), 4+ (T = 0), 2+ (T = 1) and 3f (T = 1) states of 20Ne. By taking the leading SUs states as initial vectors, quite accurate low-lying eigenvalues and eigenfunctions are obtained with a small-sized truncation space. The T = 0 states mentioned allow a three dimensional truncation to yield lowest eigenvalues with an error of less than 0.05 MeV as well as corresponding eigenfunctions having an overlap of greater than 99.5 ‘A with the exact ones. Special eigenvalue-problem techniques are discussed in the appendix.
1. INTRODUCTION
Although many complete shell-model calculations have been successfully performed for lp-shell nuclei (I), only limited successhas been met by employing the same methods in the 2s, Id shell. This is mainly due to explosive increases in size of the shell-model basis. Where, especially, more than four particles or four holes are involved most calculations employ a basiswhich has been truncated to include only what is thought to be the most important states of somesubmodel, such as the j-coupling or SU, model (2). Towards the beginning of the shell the SU, model seemsto provide a workable approximation (3) while near the other end, where holes are under consideration, the &coupling model is more often used (4). There seems to be no satisfactory truncation which works well in both regions. This last fact points out the need for exact calculations. * Permanent address: Hosei University, Koganei, Tokyo, Japan. + NRC Post-Doctoral fellow, address after June 1968: Los Alamos Scientific Laboratories, Los Alamos, New Mexico. 100
NUCLEAR
SHELL
MODEL
101
Recently Halbert et al (5) have performed exact shell-model calculations for as many as six particles in the nuclear 2s-Id shell, using forces given by Kuo (6) and Kuo & Brown (7). The dimensions of the shell-model spaces become very large reaching 537 for the 3+ (T = 1) states of 22Ne and 22Na. Although the accomplishment of this feat is encouraging it is still strongly manifested that a refinement of shell-model calculating techniques is required. Reasons for this are abundant and clear as soon as one starts to do large shell-model calculations. Unlike the atomic case, the nuclear residual force has not been very well parameterized so that a great deal of data fitting by repeated diagonalization of the Hamiltonian is inevitable. At best this process is quite time consuming. Then too, in the middle of the 2s-Id shell one encounters bases with definite spin and isospin containing more than 6000 functions (8). In our present state of computer technology this problem cannot be solved efficiently. In each nucleus the number of observed levels having the same isospin T and spin J is small so that only a few of the many shell-model states with the same T and J can be compared with experiment. The method proposed here permits us to investigate the properties of these low-lying states to any degree of accuracy. Our method is based on the Lanczos process of “minimized iterations” (9). This algorithm starts from an initial function or vector $1 and converts the original eigenvalue problem into one of finding the eigensolutions of a tridiagonal matrix. Lanczos applied the method to differential and integral operators, showing that the eigenvalue problem associated with the lowest eigensolutions of these operators may be solved very accurately by diagonalizing small matrices, if the initial function is chosen properly. The main aim of this work is to show that the nuclear Hamiltonian may be reduced to a small tridiagonal matrix with very little resulting error so long as we confine ourselves to studying low-lying states. The basic formulation of the method will be given in Section 2 and an error analysis in Section 3. Numerical calculations were applied to the Of (T = 0), 2+ (T = 0), 4+ (T = 0), 2+ (T = 1) and 3+ (T = 1) states of 20Ne in order to know how effective the method is. The results will be discussed in Section 4.
2. MATHEMATICAL
FORMULATION
The first step in our procedure is to operate with the Hamiltonian H on some unit vector +I (henceforth called the initial vector) to generate another unit vector, c#J~, orthogonal to r& . A new unit vector, &, orthogonal to @I and C& is derived next by the action of H on +2 , and so forth. Once the initial vector is chosen the sequence of operations can usually proceed until the complete space of mutually independent vectors is exhausted. The operations needed to generate
102
SEBE
AND
NACHAMKIN
the sequence of vectors +i in an N-dimensional the following equations H+,
Hilbert
space are described by
= h,+, + K,+, . . . . . . . . .
H9, = Kz--lb-l + h&v + K&c+1 . . . . . . . . . . . . . . . . H%v = &--14+-I
(1)
+ h&v,
where hi = (+i I H I 44,
(24
Ki = <9i+1 I H I +i)-
(2b)
and
It should be noted here that if +1 is orthogonal to any eigenvector of H, then the series of operations described by Eqs. (1) and (2) terminates after m steps, i.e. for some m < N, Km = 0. An important case for which this may happen is when H has some degenerate eigenvalues. In this event a new unit vector may be chosen orthogonal to the +i’s already obtained, and the operations indicated by Eqs. (1) and (2) are continued and/or repeated. We shall discuss this problem in more detail in the next section so, for the time being, it will be assumed that none of the Ki vanishes. We are now interested in solving the matrix equation HJ,=X+
(3)
for a few of the lower values of h. Here, the representation of the Hamiltonian, in the basis +r T& ,..., +N , is the N x N trigiagonal matrix
H=
h, KI 0 . . .
Kl 0 0 hs KS 0 KS h, KS . . . . . . . . .
: : . . . .
: : . . .
: : . . .
* 6-1
H,
: : . &?-I h,
We may now define the truncation of the basis of N vectors to a basis of n vectors by letting the new basis consist of the first n vectors &, &, ,..., +, . The Hamiltonian in our truncated basis is obviously represented by Hen), the n x n matrix in the upper-left corner of H. Henceforth truncation of a basis will have the above meaning.
NUCLEAR
SHELL
MODEL
103
The determinant det / h * 1 - H 1, which yields the characteristic polynomial of H, may be generated in an especially simple manner. Consider the sequence of functions {F’,(h)}, where F_,(X) = 0, F,(h) = 1, and F,(h) = (A - 111)F,-,(h)
- K~~lF&z(x).
(5)
Then it is easy to see by sequentially expanding with respect to the last row of the truncated Hamiltonian Hfn) that det / h . 1 - Htn) 1 = F,(X).
(6)
The functions Fl(A) are discussed again in Appendix I but we may sum up their properties by saying that the sequence {F&X)} is a Sturm Sequence (IO). Two important properties of (F&i)} are: (A)
The nth member of {F&l)) is the characteristic polynomial belonging to H(“), the Hamiltonian in the n-dimensional truncated space.
(B)
Between any two successive of F,(A).
zeros of
F,+,(X) there is exactly one zero
Property (B) follows because FK(X) is a Sturm sequence. Its proof may be found in standard texts, such as that by Wilkinson (10). Properties (A) and (B) are very important because they guarantee that the eigenvalues of the Hamiltonian in our truncated space all have the variational property. That is, we may arbitrarily cut off the sequence of operations defined by Eqs. (1) and (2) after obtaining n vectors +$ , generating Hen). At this point we are assured that the eigenvalues of H (“1, taken in ascending order, are greater than the corresponding exact eigenvalues of the Hamiltonian in the complete space. We may sharpen this result further by denoting the eigenvalues of Hfm), in ascending order, by @), q:“‘,..., v:,“‘. It then follows from (A) and (B) that qt.“’ z > 7)y+l’ > ... > Ai where Ai is the ith eigenvalue of the Hamiltonian, taken in ascending for each i < m the sequence of values Q(m) decreases monotonically and is bounded below by hi .l The possibility definitely exists, several of the lowest r]$ will have converged to a high degree of before the complete Hilbert space is exhausted. 1 This is merely a restatement of the well-known separation theorem (IO).
(7) order. That is, as m increases therefore, that accuracy long
104
SEBE
AND
NACHAMKIN
As a simple example of this phenomenon consider a Hamiltonian H with negative eigenvalues, the lowest of which, h, , is distinct from the others. Choosing an initial vector, +1 , which is not deficient in the ground state of this Hamiltonian we may expand cpl in terms of exact eigenvectors, +z”, and then operate on & with increasing powers of H obtaining
Eventually the term a,,&,“+~ will overwhelm all the others; the value of r for which the normalized right hand side of Eq. (8) is considered a good approximation to (Lr will depend only weakly on the size of the Hilbert space and will mainly depend on the separation of h, from the other eigenvalues. This value of r would, in fact, be simply related to the dimension of the truncation space needed to describe the ground state of H. To see this, it should merely be noticed from Eqs. (1) and (2) that
where the b’s are constants. Thus from Rayleigh’s principle, if an Z-dimensional truncation is required for the ground state then Z < r. By plotting the approximate desired eigenvalues versus the degree of the truncated space the calculation may be halted when, say, the greatest observed root changes by less than some arbitrary small number. Once the iterations have been stopped there still remains the question of just how accurate the approximate solutions are. Because of the simple form of the tridiagonal matrix it is very easy to apply the formulae derived for finding lower bounds on the eigenvalues in the Ritz variational procedure (11). These formulae as well as an analysis of the purity of the truncated eigenfunctions will be given in the next section. The accuracy of our approximation depends somewhat on the choice of an initial vector. We believe that some submodel is sure to give us a good initial vector but we leave further discussion of this point for the following sections. The fact that we may, in principle, solve the energy eigenvalue problem very accurately without a full knowledge of the Hamiltonian matrix is one of the distinct merits of our method. Sometimes, however, it might be more convenient to set up the entire Hamiltonian matrix in some well-known shell-model basis and then perform a transformation to tridiagonal form. In such cases there are two other familiar algorithms which will bring about this transformation, those of Givens and of Householder (10). A comparison of the Lanczos, Givens and Householder methods is made in Appendix II.
NUCLEAR
SHELL
3. ACCURACY
OF
MODEL
THE
105
METHOD
One of the major defects of our truncation procedure is that there is the possibility of missing a low-lying level. Although the probability of this occurring in a realistic problem is very small, it is easy to see how it may come about. Consider the initial vector, +I , and expand it in terms of the exact (unknown) eigenvectors of H, (Ly, where H+z” = hi+:“:
Repeated operation of H on +I yields WT) = Hr+,
= g /l,ra,+;x i=l
so that no new components appear in each successive @MT).Since the new basis functions, the +‘s, are linear combinations of the @fr)‘s, if +1 were orthogonal to a low-lying solution +,8” then our method would not yield a state which corresponds to $2;“. If a degeneracy exists, then for some i and j Xi = Xi and the combination (9) will always appear in each Cptr). Thus +rij will act as though it were a single vector and another linear combination of +,i”” and +j”” will be required. These defects can easily be overcome, however, by repeating the calculation using another initial vector and comparing the results with the original. Since the initial vectors are chosen from physical considerations it is unlikely that both initial vectors would be deficient in some low-lying state. In the case of a degeneracy the new initial vector will, in general, contain a different linear combination of the degenerate-state vectors, +ij, = bi +F + bj Jr,‘“. Solving with the new initial vector will then yield another independent linear combination of the degenerate vectors. If higher degeneracies are expected this process can be repeated, as in the example given in the next section. It should be realized that the diagonalization of the Hamiltonian in a truncated space defines a variational calculation. Formulae have been established for such calculations which yield upper and lower bounds on the eigenvalues as well as an estimate of the accuracy of the eigenvectors. As we have mentioned before, when there are degeneracies our truncated basis is divided into subspaces containing no degeneracies. The Hamiltonian may, in fact, be diagonalized separately in each of these subspaces. We will now assume
106
SEBE AND NACHAMKIN
that none of the eigenvalues we seek is degenerate.2 It is convenient to define what turns out to be an important concept. Denote the exact eigenvalues in the full space of dimension N by Xi and approximate ones in the truncated space of dimension I by Q, where
Now if
for i < A4 and M < Z then the sequence of eigenvalues h, , A2 ,..., & will be called separated with respect to the truncation to the space of dimension I. For separated eigenvalues we can make use of the analysis given by T. Kato (II), which reduces to a simple form in our case. Let the Z(N) approximate (exact) eigenvectors in the truncated (complete) space be denoted by ai( where (oJ+:“) belongs to an approximate (exact) eigenvalue Q(&). The quantity Q is an upper bound for Xi as mentioned in Section 2. A lower bound for hi is given in terms of the quantity Q (30) where ci2 = (4$ 1(H - qJ2 I mi)
(124
= (Qi 1H2 1Qi) - vi2.
(12b)
To obtain ci2 we first expand Qi in terms of the vectors & defined in Eq. (1). (13) k=l
The Hamiltonian can be expressed as the sum of a part acting only in the truncated space and a part which takes only the last vector, +r, out of this space. More explicitly, the first part is diagonal in the truncated space and yields viQIi when H is applied to mi while the second part transforms +r into K,+,+, , where the quantity K1 is the off-diagonal element defined in Eq. (2b). * Wilkinson (10) considers tridiagonal matrices where some of the eigenvalues are pathologically close to each other even though none of the off-diagonal matrix elements are small with respect to the diagonal elements. This situation, when it obtains, is always detectable because the Sturm sequence property (see Appendix II) allows us to count the roots lying between two numbers, no matter how close the roots are to each other. A simple remedy, should this happen, is to choose a new initial vector and re-tridiagonalize H. Almost invariably a small off-diagonal matrix element will appear, which may be set equal to zero with small resulting error.
NUCLEAR
SHELL
107
MODEL
Thus
and therefore (Qi 1H” 1 *J
= 7; + 1uIiKI j2.
(15)
From Eqs. (12) and (15) we obtain ei = j ariKr I.
(16)
If Eq. (11) holds for vi we may introduce the new quantity relation -qi < @h< hit,.
(Y~which satisfies the (17)
Kato shows now that (II) (18) The most precise (or greatest) lower bound is obtained by replacing oli by its maximum value, h,+l . When ci is very small it becomes permissible to take (~l~= qi+l instead of X,+l . In this case
Ai 3 xi = .qi -
Q2 . %+1- qi
(19)
If it happens that
then Ai will be called well separated with respect to the truncation hi 3 Yi = 7)i -
Ei s
and (21)
The quantity Y, gives a rather crude lower bound since Yi < Xi when ci is small. It should be noted that if a required eigenvalue is not well separated with respect to the truncation the quantity Yi as well as Xi may not yield a true lower bound for it. This may be seen in Figures l-5 where the lower bounds Yi are plotted as dotted curves. For small truncation spaces some of the apparent lower bounds for the higher eigenvalues are at first greater than the exact eigenvalues but after the eigenvalues become well separated the calculated lower bounds become true lower bounds and start to increase. The quantity qi is always an upper bound for Ai . Variations of Yi as the degree of the truncation space increases are discussed more precisely in Appendix III.
108
SEBE AND
NACHAMKIN
A measure of the accuracy of the eigenvectors in the truncated space is given by the quantity cio, where
The formula for the upper limit of l io as given by Kato (11), when the errors in approximate eigenvalues are sufficiently small, reduces to
where xi is the smaller of the two quantities (~+r - Q) and (Q - qiml). In conclusion, for a truncation to the space containing I vectors a complete variational error analysis may be performed with the knowledge of KI .
4. APPLICATION
TO
THE
NUCLEAR
SHELL
MODEL
Although our main interest lies in developing a method for attacking largerdimensional shell-model problems, it is perhaps best at this point to illustrate the main features of the present method with some simple examples. For these we have taken the 0+ (isobaric spin T = 0), 2+ (T = 0), 4+ (T = 0), 2f (T = 1) and 3+ (T = 1) states in 20Ne. For the first three sets the low-lying level spacings are fairly large while for the last two they are rather small. The dimensions of the untruncated 2~Id shell-model spaces for the above states are 21, 56, 44, 66 and 69 respectively. The nuclear Hamiltonian is taken to be (24)
where the first term is the single-particle spin-orbit interaction, the second expresses the splitting of the single-particle 2s and Id levels, and the last term is the residual two-body interaction. Values for [ and 77 are chosen so as to reproduce the observed energies in 17F and 170 (22), e.g. 5 = -2.03 MeV and 7 = - 1.27 MeV. The interaction Vij is assumed to be composed only of central forces having a gaussian radial form with range parameter equal to 0.7 (23), multiplied by an exchange mixture determined by p
zz.z0 ,
V13 = -60
MeV 9
V31 = -45 MeV ,
VW = 22.5 MeV.
(25)
This Hamiltonian was used by Akiyama et al (14) and yields results which are in reasonable agreement with observations in the 2s-ld shell.
NUCLEAR SHELL MODEL
109
Figs. 1,2 and 3 show the variations of the approximate eigenvaluesvi (solid lines) and the lower bounds, Yi , calculated from Eq. (21) (dashedlines), asthe dimensions of the truncated spaces increase. The graphs correspond to the four lowest 0+ (T = 0), 2+ (T = 0) and 4+ (T = 0) levels. For the sake of comparison the results of a full shell-model calculation and two other methods of truncation, to be explained later, are labelled (A), (B) and (C), respectively, and included in each figure. Using the leading SU, states 1[4] (80) J = L) as initial vectors, surprisingly rapid convergence is noted for the lowest level (first state) belonging to each spin; three-dimensional spaces allow us to obtain eigenvalues with an error of less than 0.1 MeV, and eigenfunctions whose overlaps with the exact ones
2 4 6 6 IO lZ;l4 I6 I8 m, n (DIMENSION OF TRUNCATION
SPACE)
FIG. 1. Approximate eigenvalues vi (solid lines) and calculated bounds Yj (dotted lines) for the first four O+ (T = 0) levels in 20Ne, plotted against increasing dimension of the truncation space. Levels labeled (A) are the exact eigenvalues hi . The sets (B) and (C) correspond to truncations arising from the supermultiplet and &coupling models, as explained in the text.
2
4
6
8
IO I2 I4 16 I8 20 n
FIG. 2. Same as Fig. 1 for 2+ (T = 0) levels.
110
SEBE
AND
NACHAMKIN
(T=O)
4+
711 5 kY,
I
,
f
I
2
4
6
8
FIG.
;-
+--J
3.
I
I
,
I
I
,
(A)
(8)
6)
IO 12 14 16 18 20
n
Same as Fig.
1 for 4+ (T = 0) levels.
exceed 99.5 %. For the second states the above accuracy was reached at n = 7 for O+, II = 11 for 2+ and n = 8 for 4+. Two more basis functions are necessary in each case to determine the third states to the same accuracy. As far as the first three levels with these spins are concerned, convergence does not seem to be strongly dependent upon the size of the complete shell-model space. The fourth 2+ level, however, converges more slowly than the fourth O+ or 4+ states, probably because the spacing between the fourth and fifth 2+ levels is smaller. In general the values for Yi do not start out as true lower bounds but show a decreasing trend as n increases. Eventually Yi becomes a true lower bound for the ith level and begins to increase as n increases, finally coming very close to its corresponding eigenvalue ~7~. In practise, then, we determine convergence by investigating the Yi’s as well as the 7;s. In order to display and clarify some of the characteristics of our method of solution, the values of ci in Eq. (16) are given in Table I for the lowest two Of (T = 0), 2+ (T = 0) and 4+ (T = 0) states for 12= 5, 10, 15 and 20. TABLE
I
VALUES OF THE QUANTITY ci, IN MEV, AS GIVEN BY Eqs. (12a) AND (16), FOR T = 0 STATES IN 20Ne. THE TRUNCATION-SPACE DIMBNSION IS DENOTED BY n
n
5 10 15 20
01+ 2.0 3.5 3.4 9.3
x x x x
10-l lo-’ 10-a 10-l’
a+
?A+ 2.6 3.1 x 10-Z 1.4 x 10-b 1.5 x 10-10
1.5 8.5 1.7 3.5
x x x x
10-l 10-a 10-e 10-e
A+ 3.0 4.5 x 10-l 8.9 x 10-a 2.3 x 1O-4
41+ 1.5 6.1 8.7 3.8
x x x x
10-l lo-* lo-’ lo-lo
4s+ 1.9 1.5 x 10-l 3.4 x 10-a 1.5 x 10-S
NUCLEAR
SHELL
MODEL
111
As already mentioned, the quantities ei give us measures of the accuracy of the approximate solutions. For the ten-dimensional truncation, that is n = 10, values of l i for the second 2+ and 4+ states have fallen to 0.45 MeV and 0.15 MeV, respectively. Using these quantities in Eqs. (18) and (23) one finds that for the second 2+ (T = 0) state: 9.619 MeV < X, < 9.730 MeV, (CD21i/y)
> 0.9725
(264 (26b)
and for the second 4+ (T = 0) state: 11.963 MeV < h, < 11.972 MeV, @j2 I t&> > 0.9976.
(274 W’b)
These values should be compared respectively with the exact values; 9.715 MeV, 0.9834 for 2+ (T = 0) and 11.970 MeV, 0.9999 for 4’ (T = 0). Thus we may conclude that the more accurate formula, Eq. (18), is useful for sufficiently small values of ei . (In the practical examples above, .45 MeV and .15 MeV proved to be sufficiently small.) Let us now return to the other methods of truncation referred to in Figs. 1, 2, and 3. Levels are labelled (B) and (C) corresponding to Wigner supermultiplet theory (1.5) and jj coupling, respectively. If only the basis functions belonging to the maximum symmetry [4] states were used in supermultiplet theory we would incur fairly large errors in the wave functions, exceeding 6 y0 even for the ground state of 20Ne. For this reason we follow the procedure of Inoue et al (3) and include also basisfunctions belonging to [3, I] symmetry. Functions with this symmetry interact with those of [4] symmetry through the spin-orbit interaction. Now the dimensions of the Of (T = 0), 2+ (T = 0) and 4+ (T = 0) spacesare 8, 22 and 18, respectively. This improvement of the supermultiplet model yields very good results for the lowest two states of each nuclear spin but, as seenin Figs. 1, 2 and 3, it begins to break down for the third states. Comparing with the exact wavefunctions the third O+, 2f and 4+ levels have overlaps of 0.992, 0.987 and 0.759, in that order. The fourth levels, however, cannot be trusted at all in this approximation. The truncated jj-coupling model we considered embraces only configurations of the form (I&$ (2s& (l&)7, where r < 1. In other words no configuration has more than one particle in the Id,,, orbit. Dimensions are now 10for 0+ (T = 0), 29 for 2+ (T = 0) and 25 for 4+ (T = 0). From the figures it is apparent that this approximation is rather a poor one compared with the supermultiplet theory;
-~__.
112
SJ3BE AND
NACHAMKIN
2+(T=
2-4682
II
14 16 18 20
n
FIG. 4. Same as Fig. 1 for 2+ (7’ = 1) levels. Greek letters AI and AZare omitted for clarity. 3+ (T= I 1
2
4
6
8
1
IO I2 I4 I6 I8 20 ”
FIG. 5. Same as Fig. 1 for 3+ (T = 1) levels. Most of the dotted Y, curve has been omitted for clarity. Also, for reasons of clarity, greek letters ha and h, have been omitted.
e.g. the overlaps of the lowest eigenfunctions with the exact ones are only 0.878 for O+,0.911 for 2+ and 0.918 for 4+. The eigenvalues rli and lower bounds Yi for 2f (T = 1) and 3+ (T = 1) states in z”Ne are given in Figs. 4 and 5, with the same notation as in Figs. l-3. An interesting phenomenon is displayed in Figs. 4 and 5. The lowest 3+ level converges faster than the lowest 2+ level even though the complete space for 3+ functions is larger than the complete 2+ space. Initial vectors were chosen here to correspond to the states 1[31] S = 0 J = L), simply for convenience in computing. The convergencies are rather slow as compared to the T = 0 cases.This probably arises, in part, becausethe spacings between adjacent levels are smaller in the T = 1 case. The lack of a proper initial vector is, however, also a cause of slow convergence.
NUCLEAR
SHELL
113
MODEL
As an example of the latter cause consider the four states belonging to the representation [3,1] (6,1), i.e. / S = 0, L = J), I S = I, L = J - I), 1S = I, L = J), j S = 1, L = J + 1). The expectation value of the Hamiltonian is lowest in the state with S = 1, L = J - 1. Taking an initial vector corresponding to this state, instead of the one actually used in Figs. 4 and 5, allows an accurate truncation with as few as seven vectors to describe both the 2f and 3+ lowest states. If all four lowest levels are required for both spins then it turns out that there is no remarkable difference between the two different initial vectors. Truncating at n = 20 yields overlaps between the exact and approximate eigenfunctions of greater than 99.5 ‘A for the levels shown in Figs. 4 and 5. The values of l i for T = 1 are tabulated in Table II and show slower decreasesagainst increasing n than do those for T = 0. The supermultiplet and jj models are again labelled (B) and (C), respectively, in Figs. 4 and 5. The only difference here is that only states with symmetry [3, I] are included in the supermultiplet scheme.We did this becausestatesof symmetries [2,2] and [2, 1, l] may couple directly with [3, I] through spin-orbit and two-body forces, making it difficult to obtain a truncation of proper size for a comparison with our calculations. In this model both the 2+ and 3+ spacesare 24-dimensional. Employing this supermultiplet model the lowest 2f state has an overlap of 0.951 with the second lowest of the exact calculation while the second lowest has an overlap of 0.921 with the lowest exact state. This means that the approximation predicted an inverted order for the lowest two levels. For 3+, the second lowest eigenfunction has only a 0.776 overlap with the exact one. Taking the samejj configurations as for the T = 0 casethe dimensions are 36 for the 2+ spaceand 38 for the 3+. In spite of the large spacesin this approximation, it does not yield very good results; overlaps with the exact eigenfunctions do not exceed 0.97 for any of the 2f or 3’ levels shown. TABLE SAME
2,+
n
5
AS TABLE
I FOR
2,+
II T =
29+
1 STATES
3,+
IN
20Ne
3,+
3,+
7.9 x 10-l
1.9
2.6
10
5.8
1.5
1.2
1.2 7.4 x 10-l
2.9 8.3 x 10-l
4.0 3.9 8.5 x 10-l
x IO-’
15
2.3 x IO-’
2.6 x 10-l
1.1
7.9 x 10-z
2.6 x 10-l
20
1.1
x
1.7 x 10-S
1.7 x 10-l
3.7 x 10-a
4.0 x 10-z
2.2 x IO-1
25
1.1
x 10-a
2.1 x 10-a
3.9 x 10-a
4.1 x 10-z
9.9 x 10-S
8.4 x 1O-s
10-Z
114
SEBE AND
NACHAMKIN
We would again like to point out that, by using our method, the size of the truncated space needed to obtain a good approximation depends mainly on the level structure of the states of interest and only weakly on the full size of the shell-model basis, while other methods of truncation yield bases which rapidly increase in size as particle number increases. As a further example we consider a case wherein a degeneracy arises, by taking our Hamiltonian to consist merely of the single-particle spin-orbit force. Four different initial vectors were chosen for 2+ (7’ = 0) of 20Ne, namely 1 [4] (8, 0) S = 0 L = 2), 1 [4] (0, 4) s = 0 L = 2), j [3, l] (4, 2) s = 1 L = 3), and I [3, l] (6, 1) S = 1 L = 3). Lanczos iterations were stopped in each case when an off-diagonal matrix element became very small, and the resulting tridiagonal matrix was solved for its lowest eigenstate. In all cases this state turned out to have the same energy as the configuration (&j4. Calling the four states so obtained $1 , I&, I,& and z,& the 4 x 4 symmetric matrix with elements Mij = (J,& j 4,) was constructed: 0.8723 1 -0.7658 Quoting
0.1531 -0.0640 1 i -0.3451
I.
(28)
results to four decimal places detM} = 0.0000, det(M(4 14)) = 0.0788 # 0,
(294 (29b)
where M(4 x 4) is the same as M but with the fourth row and fourth column deleted. Eq. (29a) is the criterion that +I ,..., $4 are linearly dependent while Eq. (29b) tells us that I+&& and z,& are linearly independent. Deleting other row-column combinations also gives a non-vanishing determinant. This simple but important example shows how a degeneracy, in this case of third order, may be resolved by assuming new initial vectors. (Even though it is highly improbable that a degeneracy would not be resolved in the way described above we feel that the method of choosing new initial vectors could well become an art based on intuition.) 5. CONCLUDING
REMARKS
As has already been shown our method has several characteristics which make it suitable for solving for the lower-lying states in large shell-model calculations.
NUCLEAR
SHELL
115
MODEL
Restated, they are: (1)
The special iterational character, which suggests that the size of the truncated space will be almost independent of the size of the original shell-model space;
(2)
A very simple set of error formulae, safely accurate truncation;
(3)
The time consuming eliminated.
diagonalization
which
allows
an economical
of large matrices
yet
is virtually
Another question is what is the best method of constructing the tridiagonal Hamiltonian matrix directly. As mentioned before, it is not really necessary to generate the entire matrix in the complete shell-model basis. In principle the I x I tridiagonal matrix of H may be expressed in terms of the expectation values of H, HZ,..., Hz’ evaluated in the initial state. T.D. Newton attacked the problem of generating the tridiagonal matrix directly, but since he was attempting full calculations using a Slater-determinant basis, the large problems were out of his reach (26). In all the examples we have given, the entire Hamiltonian was set up and then tridiagonalized. We feel that this may not be the best way to advance to larger-dimensional problems. At present work is proceeding on problems in the 2s-ld shell which involve more than four particles or holes. Finally but perhaps most importantly, we have seen that the properties of many low-lying states may be grasped by using very small truncation spaces. Since convergence in these spaces depends, to a great extent, on the submodel from which the initial vectors are chosen much can be added to the physical interpretation of these states. Choosing an initial vector from, say, the SU, model, the jj-coupling model, the Nilsson model (17) or the model proposed by Harvey & Sebe (18), and watching for the rapid convergence of the resulting spaces is a fairly stringent test of these models’ efficacy. In short our method supplies a testing ground for the various nuclear models.
The authors would like to express their thanks to Dr. M. Harvey for many stimulating discussions and to Prof. T. D. Newton for some incisive remarks. We are also indebted to Mr. H. H. Clayton for his kind interest in our work. One of us (T. S.) wishes to thank Mr. G. C. Hanna for the hospitality extended to him at the Chalk River Nuclear Laboratories.
116
SEBE AND NACHAMKIN
APPENDIX
I.
EIGENSOLUTIONS
OF TRIDIAGONAL
MATRICES
In this appendix we consider the problem of numerically computing the eigenvalues, hi, and the eigenvectors, Gi, of a real, symmetric, N x N tridiagonal matrix. We will be making use of the following property of Sturm Sequences (see Section 2 for properties (A) and (B)): (C) For any pO the number of times the Sturm Sequence {F&J) changes sign is equal to the number of eigenvalues of H which are greater than p,, . Starting from initial upper and lower bounds each hi is derived by the bisection method (10). We feel that the method of bisections is of enough interest to be explained in detail. The following steps define an algorithm which will isolate the eigenvalues of a tridiagonal matrix quickly and with high accuracy.
(1) Set initial upper and lower bounds on the eigenvalue desired, defining an interval. (2) Compute the Sturm sequence, using Eq. (5), at the midpoint of this interval, counting the number of times it changes sign. (3) From the results of step (2) the eigenvalue is determined to be in either the upper or the lower half interval, so that a new interval is defined which is half the size of the current interval. (4) If the new interval is smaller than some predetermined size then the sought-after eigenvalue may be set equal to the number corresponding to the middle of this interval otherwise go back to step (2). Step (1) involves finding upper and lower bounds on given eigenvalues without solving for them. In many physical problems this might be fairly easy but we make use of a simple infallible method to obtain reasonable upper and lower bounds for each eigenvalue of a tridiagonal matrix. To do this we make use of Gershgorin’s theorems which state that (IO): I) For a matrix H with elements Hii , each interval (or disk in the complex plane, if H has complex eigenvalues) defined by I Hii - P I < Gi ,
U-1)
where Gi = c I HijI, j=i contains at least one eigenvalue of H, and II) Each eigenvalue of H lies in at least one Gershgorin disk.
U-2)
NUCLEAR
SHELL
117
MODEL
If H is tridiagonal Gi never involves summing more than two terms. For ease of computation in the real, symmetric tridiagonal case we choose the off-diagonal matrix element of largest absolute size, calling it bm,, . Next, the diagonal elements are reordered according to algebraic size so that the nth one, a, , determines the interval an - 2bmax < p < a, + 2bmax,
(1-3)
and Gershgorin’s theorem guarantees our finding the nth eigenvalue in this interval. As the calculation for the eigenvalues proceeds (Steps (l)-(4) above), new bounds on subsequent eigenvalues may be obtained if the so-defined intervals overlap. Once the desired eigenvalues of H are evaluated we next seek the corresponding eigenvectors. If we know some eigenvalue, Xi, exactly and compute the Sturm sequence {F&J} exactly then it is easy to verify that the vector Gi, with components 01~, 01~,..., 01~, where =
FN-l(xi),
OIN-1
=
-HN.N--1FN--2@i)v
OIiV-7:
=
t-1”
OIN
(1-4)
HN,N-lHN-l.N-2
*‘*
HN-k+l,N-kFN-k-l(xi),
satisfies the matrix equation (Ai + 1 - H) +,i = 0
(1-5)
In the general case we cannot obtain the exact value Xi but must be satisfied with an approximation qi , for which there is no nontrivial solution of
We may, however,
(Q * 1 - H) tj+ = 0.
(I-5a)
(pi . 1 - H) +i = Xi
V-6)
solve
where xi is nonnull, obtaining the approximate
eigenvector
& = (Q . 1 - H)-l
xi
& , where (I-7)
If we now take xi = & and solve Eq. (I-6) again the new +i will be an even better approximation to the exact eigenvector +,i, if Q is close to hi. To see this, expand xi in terms of the exact eigenvectors, Xi =
C a&
W)
118
SEBE AND NACHAMKIN
Then (Q * 1 - H)-7 xc = c a&/(~
* 1 - Ai)
(I-9)
and Gs will dominate the series if, for i f j, I qi - hi I Q I r)e - h I.
(I-10)
Repeating the process purifies the eigenvector very quickly if Eq. (I-10) holds and xi is not deficient in +i . Following Wilkinson fairly closely the purification process described above is performed four times (Wilkinson recommends three), solving the linear equations by elimination and back substitution (IO). For approximate eigenvalues correct to four or five significant digits this process yielded eigenvectors, for all the physical cases tested, which were so accurate that the eigenvalues could be recalculated to ten or twelve significant digits. APPENDIX
II.
RELATIVE MERITS OF TRIDIAGONALIZATION
PROCESSES
Since the Hamiltonian is quite conveniently expressed when in tridiagonal form we feel that more should be said about reducing it to this form. The three methods usually employed for the transformation are the Givens (29), Householder (19), and Lanczos algorithms (9). The first two require the Hamiltonian to be in matrix form from the outset, that is, a “parent” matrix must be known. Theoretically, the last method is capable of generating the tridiagonal “daughter” matrix directly by knowing the expectation value of powers of the Hamiltonian in an initial state. Most frequently the parent matrix is generated. In realistic shell-model problems this is usually a large sparse matrix whose low-lying eigensolutions are required. There seems to be no doubt that if the parent matrix is specified, the Householder method is superior to both the Givens and Lanczos methods in achieving total tridiagonalization, for reasons of storage and speed. Accuracy of all three is high enough to be of minor consideration (10). The Givens method proceeds via two-dimensional rotations, the (i, j)-th element being eliminated by a rotation in the (i + 1, j)-th plane. An additional matrix, the same size as the parent matrix, is needed just to store the details of the transformation. Thus, even if the N x N parent matrix is real and symmetric, a total of N(N - 1) + 1 storage locations is needed to completely tridiagonalize it and retain all the information concerning the transformation. At each point of a Lanczos matrix iteration a new N-dimensional vector must be stored so that a total of N(N + 2) - 1 storage locations must be allocated for total tridiagonalization. Since Householder’s tridiagonalization requires only $N(N + 3) locations, it can accommodate much larger matrices.
NUCLEAR
SHELL
119
MODEL
As far as speed is concerned, the Householder & Givens processes require respectively approximately $V3 and 4N3 multiplications as well as N and $(N - l)(N - 2) square roots for tridiagonalization, while the Lanczos method requires a reorthogonalization process at each step and thus requires a total of approximately as many multiplications and square roots as the Givens method. After the tridiagonal matrix has been solved for eigenvectors the final eigenvectors of the parent matrix are found by the Givens method at only half the speed of either the Householder or Lanczos method. Even when an exotic unit vector + is to be used as an initial vector in the Lanczos algorithm, it is far more economical to first transform H by the unitary orthogonal matrix P where P = 1 - (+ - eJ(+ - elY/[l
- (e,, +)I,
(II-l)
where XT is the transpose of X, and H’ = P-lHP
= PHP.
(11-2)
Each step of the Householder method applied to II’ is equivalent to the corresponding Lanczos iteration. It should be noted that P was chosen in that form because P+ = e, (11-3) and a Householder routine may easily be adapted to perform this transformation. Only N extra storage locations are required if this recommendation is followed. If the tridiagonal representation of the Hamiltonian is to be generated directly by means of Lanczos iterations the exciting possibility exists that, given an appropriate initial vector, the number of iterations needed for very accurate determination of the low-lying states may be very small. In the future if a variational truncation based on the Lanczos procedure is not used, there may be no alternative to diagonalizing unreasonably large matrices in shell-model calculations. We have already mentioned that the actual numerical errors encountered by tridiagonalizing a matrix by the methods mentioned were very small. To get an idea of the errors introduced we will consider the percentage mean-square deviation, x, (11-4) X = Cithi l)J2/Cihi2~ where hi and qi are the exact eigenvalues of the parent matrix and the corresponding exact eigenvalues of the daughter tridhgonal matrix. Since the computer performing the algorithm does not usually work with infinite accuracy, x f 0.
120
SEBE AND NACHAMKIN
Wilkinson gives what turn out to be overly large upper bounds for x in the Givens & Householder reductions (10). They are, respectively, x(Givens) < 12 x Naj2 x 231 + 6 x 2-t)4N-7, x(Householder) < 40 x N x 2-31 + 20 x 2-t)2N,
(II-5a) (II-5b)
where t is the binary bit size of the numbers used by the computer. For the Chalk River G-20 computer, 2-t M lo-l2 for double precision arithmetic. Assuming N = lo4 (a very large matrix) then x(Givens) < 1.2 x 10-5, x(Householder) < 4 x 10-7.
(II-6a) (II-6b)
Such a matrix would take us more than 5 x lo4 hours just to reduce to tridiagonal form! Just as a comparison, the six-sweep Jacobi method has the bound (IO) x(Jacobi) < 108 x 2-t x N312(1 + 9 x 2-t)6(4N-7) = 1.08 x 10-4 (11-7) We have not made a detailed error analysis of the Lanczos process with reorthogonalizations but, since approximately the same kinds and number of steps are used as by the Givens process, the error should correspond closely to x(Givens). Thus, as far as numerical errors involved in calculating the eigenvalues are concerned, any shell-model matrix which could be fed into the memory of almost any electronic computer today could be solved for low-lying eigenvalues with high accuracy. (In fact statistical cancelling of rounding errors will improve the bounds in Eqs. (11-6) and (11-7) by several orders of magnitude.)
APPENDIX
III.
BEHAVIOR
OF THE LOWER BOUND
Yi
The method we have proposed has the property, common to most iteration processes, that the lower the eigenvalue the faster its approximation converges. We shall therefore assume that we have achieved an I-dimensional truncation space wherein the first k - 1 solutions have converged to within a preassigned accuracy and the kth approximate eigenvalue &) is fairly close to the exact value hk . Under these assumptions the kth approximate eigenvector 0:‘) may be expanded very accurately in terms of exact eigenvectors +;“, with eigenvalues hl > hl, , i.e. (III-l)
NUCLEAR
SHELL
121
MODEL
(111-2)
and the quantity
eka in Eq. (12.b) is expressed as (111-3)
If we assume that the curve of Yk (=rii) - l 3 plotted against Zdefines a reasonably smooth function Y,(Z), then Eqs. (111-1-111-3) may be used to obtain the tangent to this curve as follows:
dYk ___
(III-4a)
dZ
(III-4b)
t qy ,=f&[F-h,? + 2(g) + ‘,
We now make the assumption that all the h’s are negative and that we are in the region where A, and A,,, are well separated, i.e. (111-5) Then, obviously [-A,a
+ 2(7&f’ + Ek) A,] < [-A;,,
Replacing all the Al in Eq. (III-4b) 1 --dY, dZ > 2Ek
+ 2(4? + 4 &+,I, by A,+1 it follows
[-X,2
1 > k + 1. (111-6)
that
+ 2(?$’ + Ek) x,1 y
+ [ -h;+l +
37:’
+ Ek)~,+,IL=t+lqgq
2
(111-7)
since we are assuming that
d(Q) -x-
for
1 > k.
(111-8)
122
SEBE AND NACHAMKIN
Then (111-9) The right-hand side of Eq. (111-9) will now be positive if 72’ + Ek < (A, + 4+,)/2,
(III-1Oa)
$j’ + Ek > (A, + &+JP.
(III-lob)
and will be negative if
In other words the quantity Yk may, on the whole, be expected to be a decreasing function of Z roughly until the point is reached where Eq. (III-1Oa) is true. After that Yk increases monotonically with increasing I.
APPENDIX
IV.
A COROLLARY OF GERSHGORIN'S FIRST THEOREM
In Appendix 3 mention was made of the Gershgorin disk theorem wherein at least one eigenvalue of a matrix H may be found in a disk centered at Hii and of radius C:‘i 1Hij I. This theorem may be improved upon slightly in that the radius may be taken to be Gj = cc’,+” 1Hij 12)1/2.
To show this, consider the ith first row and column so that the tively while the eigenvalues of Householder procedure produces
row. Exchange the ith row and column with the elements Hij and H,, become Hij and Hj; respecthe matrix are unchanged. The first step in a a matrix H” where (IV-2a)
Hii = %I , H,“,
(IV-l)
= ztCZj>~ I HI, W2
(IV-2b)
and H& = 0,
(n > 2).
(IV-2c)
Applying Gershgorin’s theorem to H”, which has the same eigenvalues as H, immediately yields an eigenvalue in the disk centered at H,, with radius given by Eq. (IV-l). It should be noted that these new disks do not necessarily satisfy the second Gershgorin theorem. RECEIVED:
June 27, 1968
NUCLEAR
SHELL
123
MODEL
REFERENCES
1. For example: D. R. INGLIS, Rec. Mod. Whys. 25, 390 (1953); D. KURATH, Phys. Reu. 101, 216 (1956); HILL, AND
S. COHEN AND D. KURATH, Nucl. Whys. 73, 1 (1965); J. NACHAMKIN, Nucl. Phus. A106, 62 (1968).
P. GOLDHAMMER,
J. R.
Proc. Roy. Sot. A245, 128, 562 (1958); J. P. ELLIOTT AND M. HARVEY, Proc. Roy. Sot. A272, 551 (1963). 3. For example: M. BARNERJEE AND C. LEVINSON, Phys. REV. 130, 1036 (1963); T. INOUE, T. SEBE, H. HAGIWARA AND A. ARIMA, Nlrcl. Whys. 59, 1 (1964); M. HARVEY, The nuclear SU, model. Advances Nucl. Phys. 1, 67 (1968). 4. For example: P. W. M. GLAUDEMANS, G. WUCHERS, AND P. J. BRUSSARD, N~cl. Phys. 56, 2. J. P. ELLIOTT,
529, 548 (1964). 5. E. C. HALBERT,
6.
J. B. MCGRORY,
T. T. S. Kuo, N~cl. Phys. A103,
AND
B. H.
WILDENTHAL,
AND G. E. BROWN, Nlrcl. Phys. T. SEBE AND M. HARVEY, Chalk River Report C. LANCZOS, J. Res. Nat. Bur. Stand. 45, 255 J. H. WILKINSON, “The Algebraic Eigenvalue T. KATO, J. Phys. Sac. Japan 4, 434 ( 1949).
7. T. T. S. Kuo
8. 9. 10. II. 12. 13. 14. 15.
F. AIZENBERG-SELOVE R. THIEBERGER, Y. AKIYAMA, E. P. WIGNER,
AND
Phys. Rev. Letf. 20,
1112
(1968).
71 (1967).
T. LAURITZEN,
Nuc/.
85, 40 (1966). AECL-3007 (1968). (1950). Problem.” Clarendon Press, Oxford, Phys. 11,
1965.
1 (1959).
Nuc/. Phys. 2, 533 (1957). A. ARIMA, AND T. SEBE, to
be published. Phys. Rev. 62, 438 (1942). 16. T. D. NEWTON, Chalk River Report, AECL-2223 (1965). 17. S. G. NIUSON, Dan. Mat. Fys. Medd. 29, 75 (1956). 18. M. HARVEY AND T. SEBE, to be published. 19. For example: W. GIVENS, Oak Ridge National Laboratory Report ORNL-1574; HOUSEHOLDER, J. Ass. Camp. Mach. 5, 335, 339 (1958); P. A. Wm AND R. R. Mathematics of Computation 18, 457 (1964).
A. S. BROWN,