Application of TDHF Methods to Nuclear Physics Citation Flanders, Bradley A. (1981) Application of TDHF Methods to Nuclear Physics. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/vakw-ws15. https://resolver.caltech.edu/CaltechTHESIS:03022018-134310889 Abstract Part I presents results for TDHF calculations of realistic heavy-ion reactions. Results for a separable approximation, neglecting motion normal to the scattering plane, agree very well with the results of the full 3-D calculation. Results for the fusion cross section of two systems leading to the compound nucleus 56 Ni are compared with experimental data. The 16 O + 40 Ca results agree quite well with the experimental data, and the 28 Si + 28 Si results agree quite well with experimental data for the similar 32 S + 27 Al system. Results of calculations with the separable approximation for 86 Kr + 139 La are compared with both axially symmetric calculations and experimental results. All three show substantial agreement. Part II presents a stability criterion for the validity of TDHF solutions based on a time-dependent generalization of RPA theory. Results are tested in an exactly soluble model, the SU(3) generalization of the Lipkin model. Unfortunately, the exact solution could not be computed for a large enough number of particles to permit quantitative testing of this criterion. However, it does seem to correctly indicate the stability or instability of the TDHF path. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Physics) Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Physics Thesis Availability: Public (worldwide access) Research Advisor(s): Koonin, Steven E. Thesis Committee: Koonin, Steven E. (chair) Tombrello, Thomas A. Fox, Geoffrey C. Devi, K. R. Sandhya Defense Date: 29 May 1981 Funders: Funding Agency Grant Number NSF PHY76-83685 NSF PHY77-21602 NSF ERDA #W-7405-ENG-48 NSF ERDA #W-7405-ENG-26 NSF NSF 77-12879 Record Number: CaltechTHESIS:03022018-134310889 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:03022018-134310889 DOI: 10.7907/vakw-ws15 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 10747 Collection: CaltechTHESIS Deposited By: Mel Ray Deposited On: 02 Mar 2018 23:40 Last Modified: 10 Feb 2025 20:27 Thesis Files Preview PDF - Final Version See Usage Policy. 66MB Repository Staff Only: item control page

APPLICATION OF TDHF METHODS TO NUCLEAR PHYSICS Thesis by Bradley A. Flanders In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1981 (Submitted May 29, 1981) ii ACKNOWLEDGMENTS Many thanks to my advisors, Steve Koonin and Tom Tombrello, for their help, inspiration, and encouragement throughout my study at Caltech, and, in particular, Steve, for his guidance with my thesis project. Their example has shown me the fruits of intelligence, creativity, dedication, and hard work. I am also grateful to my collaborators, especially Mort Weiss and Paul Bonche, for the ideas and knowledge they have shared with me throughout our association. My parents, Ted and Jeanne, and my sister, Karen, deserve special thanks for the love and encouragement they gave me while I was growing up and which they continue to give me. Special thanks also go to Pauline Hayashigawa for her love and understanding, for taking care of the mundane details of my life while I wrote my thesis, and for her help in typing it up. I'd like to thank my many friends, too, for making the pleasant times pleasant and the unpleasant times bearable. I can't include all their names here, but I'm especially grateful to Barb Keenan, Steve Slutz, Ken Iverson, and Steve Trentalange. And "domo arigato" to Mr. Ohshima and the Caltech Karate Club for what I have learned about the human spirit. I would also like to express my appreciation to Caltech for its financial support. This included both teaching and research assistant- ships and the following grants: Eng.-48 and 26, and NSF 77-12879. PHY76-83685, PHY77-21602, ERDA #W-7405- iii ABSTRACT Part I presents results for TDHF calculations of realistic heavy-ion reactions. Results for a separable approximation, neglecting motion normal to the scattering plane, agree very well with the results of the full 3-D calculation. Results for the fusion cross section of two systems leading to the compound nucleus 56Ni are compared with experimental data. The 160 + 40ca results agree quite well with the experimental data, and the 28si + 28si results agree quite well with experimental data for the similar 32s + 27Al system. Results of calculations with the separable approximation for 86Kr + 139La are compared with both axially syn:nnetric calculations and experimental results. All three show substantial agreement. Part II presents a stability criterion for the validity of TDHF solutions based on a time-dependent generalization of RPA theory. Results are tested in an exactly soluble model, the SU(3) generalization of the Lipkin model. Unfortunately, the exact solution could not be computed for a large enough number of particles to permit quantitative testing of this criterion. However, it does seem to correctly indicate the stability or instability of the TDHF path. iv TABLE OF CONTENTS ACKNOWLEDGMENTS ii ABSTRACT iii PART I: TIME-DEPENDENT HARTREE-FOCK CALCULATIONS OF NUCLEAR COLLISIONS 1 Section 1: Time-Dependent Hartree-Fock Theory 2 Section 2: The Effective Interaction 9 Section 3: Brief Review of Heavy-Ion Reactions 19 Section 4: Application of TDHF to Heavy-Ion Reactions 24 Section 5: Brief Review of Earlier Calculations 33 Section 6: The Separable Approximation 37 "A Separable Approximation to Time-Dependent Hartree-Fock Calculations" Phys. Lett. 77B, 13 (1978) 41 The Filling Approximation 56 "TDHI Calculations of Fusion Cross Sections For 60 + 40ca and 28si + 28si" Phys. Rev. C 20, 641 (1979) 60 Section 8: Results for 86Kr + 139La 88 Section 9: Sumnary 95 Section 7: PART II: A STABILITY CRITERION FOR TDHF CALCULATIONS 96 Section l: Introduction 97 Section 2: RPA Theory 99 Section 3: Time-Dependent Generalization of the RPA Theory 101 Section 4: TDHF for the SU(3) Model 107 Section 5: TDRPA for the SU(3) Model 112 Section 6: The Exact Solution for the SU(3) Model 117 v Section 7: How to Compare Exact and TDHF Calculations 123 Section 8: TDRPA Results 127 Section 9: Sumnary 132 APPENDICES Al. A2. A3. A4. 133 Derivation of the TDHF Equations 134 I. Variational Derivation 134 II. Derivation by Truncation 137 III. Hartree Potential and Fock Potential 139 IV. Conservation Laws 140 (H) for the Skyrme Potential 143 I. Introduction and Kinetic Energy Term 143 II. Two-Body Terms 145 III. The Three-Body Term and Density-Dependent Two-Body Term 149 Spatial Discretization of the Energy Functional 151 I. Introduction and Zero-Range Terms 151 II. Kinetic Energy Term 152 III. Nonlocal Terms 153 IV. Coulomb and Yukawa Terms 157 Finite Difference Form for the TDHF Equations in Space 158 I. Introduction and Zero-Range Terms 158 II. Kinetic Energy Term 159 III. Nonlocal Terms 159 IV. Coulomb and Yukawa Terms 162 vi A5. A6. Implementation of the Separable Approximation 163 I. General Features 163 II. Incorporation into the Energy Functional 164 III. Incorporation into the TDHF Hamiltonian 167 Discretization of the TDHF Equations in Time 172 I. Introduction and Statement of the Problem 172 II. Comparison of Taylor Expansions for the Exact and Discrete Solutions 174 Derivation of the Discrete Evolution Operator from Energy Conservation 175 Evolution Algorithm 178 III. IV. A7. Binary Partition of the Nuclear Density 181 Bl. Expectation Values Needed for TDRPA 186 B2. Invariance of the TDRPA Eigenvalues 190 B3. Expectation Values for the SU(3) Model 193 B4. Energy and Evolution Equation as Functions of the Four Rotation Angles 196 REFERENCES 199 FIGURE CAPTIONS 203 FIGURES 207 -1- PART I TIME DEPENDENT HARTREE-FOCK CALCULATIONS OF NUCLEAR COLLISIONS -2- Section 1: Time-Dependent Hartree-Fock Theory A recent major effort in theoretical nuclear physics has been the devek:Jmenl of microscopic dynamical theories to complemenl the pl1enome11ulogi:.. d , macroscopic treatments . These macroscopic treatments, compr ising V:·~-t:1c:~ l y the entire history of nuclear physics, have been very valuable as a fr .:.m ev,' '.:'rk in which to visualize various phenomena; showing what degrees c·f freedom are needed and what laws govern the motion of these degrees of freedom . However, because these models are constructed to fit experiments on a particular phenomenon, they cannot give any underlying, unified explanation for diverse phenomena, nor show fully the relationships between these phenomena. Macroscopic theories also have a limited power for predicting new phenomena. To resolve these questions, we must turn to a microscopic many-bcdy theory. The ultimate theory of nuclear physics, involving nucleons, pions, etc . is, of course , not feasible; but approximations can be made which should be valid for certain types of problems. Time-Dependent Hartree-Fock theory (TDHF) (DI 30) is one such approxL-nalion. It consists of independent particles, in this case nucleons , moving in a selfconsistent mean field generated from the internucleon interaction. This independent particle picture may seem unreasonable, since nucleons at normal nuclear densities are strongly interacting. However, since they are fermions, the Pauli principle is very effective at blocking their interaction at low energies . A simple estimate of the mean free path of a particle p, with energy El A +er, (where er is the fermi energy) may be obtained by considering binary collisions, {p)+(q) -.. (p ')+(q') , in infinite nuclear matter. The collisior.. rate including Pauli blocking is (KO 79b,PI 66) (1.1) where vP is the velocity ~ ~~ is the free-nucleon cross section, and f is the -3- phase-space occupation. These can be calculated assuming an h;c.tropic differential cross S8Ction ( ~~ - :; ), a Fermi-Dirac distribution at low t empera2 ture for the phase-space occupation, a..'1.d a free-nucleon spectrum eP = ;!-.. F or f ~ rri high energies, EI A»ep, the Pauli blocking is ineffective and the collision rate, I'p has its cla,sical value pvP a 0 where p is the density of nuclear rnatt ~-r . Wi: h typical numbers p=.17frn- 3 and a 0 =250mb, this leads to a short mean free path v >.p = p.23fm. However, for low energies, El A«eF, Pauli blocking severely p limits the number of accessible final states+ (KO 79b) so that ( 1.2) and ~= P;o ~( ;:A) 2 ( 1. 3) With eF=41 MeV this leads to a mean free path (1.4) and a lifetime, _1_= 5.9x10-21 sec, (El A) 2 ( 1.5) rp where El A is measured in MeV. For energies up to a few MeV per particle, the 2 is reasonable . For a Fermi +There is a simple way to show that this blocking factor, ~~ 4 ., .a.. ges, the number of states per unit volume with energies less than E is proportional to E 2 . The density of states near the Fermi surface : is proportional to The PauJi princip le will require 11' and q ' to have energies above 'F and conservation of energy will require that all three particles q, p ' , q' have energies within E/A of the Fermi surface for E/A<•F · For EI A small, the density of states will be appro:X:mately constant throughout this energy range . Since conservation of energy will specify the energy of one of these particles, the nu.'Il.ber of accessible states will be limited to approximately (A:~ JfE/ A) 2 . Since the number of stat es ieJf· .a.. without Pauli blocking would have been (A:e}J2, Pauli blocking will cause a factor approximately ~~ 2 . Since conservation of energy act ually requires all three particles to shar e 4 'F this deviation E/A from the Fermi surface rather than merely being bound to within E/ A of the surface, it isn't surprising that the actual blocking factor is three times smaller than above. -4 - mean free path is much larger than the nuclear size (~10 fm) an d t h e life1-L-:1e ·s com.parable to the collision time (several 10-21 sec.). Therefore , in r.luclear •Xdisions , this picture is expected to work fer energies up to a fev- J..: e'.0 p er Ll.:'.:: ! ~·..; ::-i above the Coulomb barrier . Indeed, this independent particle picture has teen quit e suc cessful in explaining static properties of nuclei, from the early days of the she ll model to the recent sophisticated, self-consistent, mean-field calculations, such as Density-Dependent Hartree-Fock (NE 70) . The first extension of thes e ideas t o / time-dependent problems was the use of the RPA approximation (TH 6 1), a sma ll amplitude approximation to TDHF, to calculate the near-equilibrium dynamics of low-lying vibrational states in nuclei (BER 75) . This success prompted the investigation of TDHF in situations far from equilibrium, though still at relatively low excitation energies . Thus , although the derivation of the TDHF meth:::>d date s back to Dirac in 1930, the application of these calculations to the collision of heavy ions was not attempted until the mid 1970's (BO 76a) . These TDHF calculations have several attractive features . First, they are the natural extensions of successful, near-equilibrium descriptions . Since they use the same internucleon interaction (see Section 2) , they have no free pararneters. This unified approach links these different problems and completely specifies the behavior of the TDHF solutions. Second, TDHF calculations treat implicitly both single-particle and collective degrees of freedom . By not imposing collective degrees of freedom, they provide a flexible treatment of the nuclear shape and surface. The mean field is the most obvious way of comrnunieating collective information. Third, TDHF is the microscopic one-body theory to compare with macroscopic treatments of the long mean free path reg ime (KO 79b) . To discuss in more detail the nature of the approximation and the type of -5- results we will obtain, we begin with a derivation of the TDHF formula. Details can be found in Appendix Al . The exact Schroedinger equation can be derived by applying a variational principle to the action (KE 76), I = J cU; <+(t) jih:t -HI +(t )> . (1.6) where H is the Hamiltonian. Looking for stationary points of I with respect to variation of the conjugate wavefunction, +•. gives ih:t + =H '¥. The transition to a derivation of TDHF is straightforward (KE 76) . Since the total wavefunction is assumed to be composed of independent particles, it can be written as a pro- duct of single-particle wavefunctions, 1/1;; or in the case of identical fermions, where'¥ must be antisymmetrized, as a Slater determinant of the 1/J; . We then derive the TDHF formulas by looking for stationary action with respect to variation of a particular single-particle conjugate wavefunction, 1/lt This leads to an evolution formula for each 1/1;, 6 411..!!_1/IJ = - -<+IH I+> Bt 6l';• (1.7) Although it isn't obvious that this leads to a mean field formulation, (1.8) with h a one-body operator; this is shown in Appendix Al.Part I. Thus, we see that the TDHF wavefunction is the best independent-particle approximation to the exact Schroedinger wavefunction However, the independent-particle assumpti.on greatly limits the usable information we can obtain from the TDHF wavefunction, since any quantity strongly dependent on correlations between two or more particles (other than antisymmetry) will not be correct. Therefore, the theory only attempts to approximate one-body operators and will give useful results for at most few-body operators. This can be seen explicitly in an alternative derivation of the TDHF formulas detailed in Appendix Al.Part JI. In this derivation, the exact equation for the -6- one-body density matrix, p~(t) =<ir(t) Ia/a" Iir(t )> . (1.9) in terms of the two-body density matrix, (1.10) and higher-body density matrices is truncated by retaining only many- body correlations required by the Pauli principle and expressing these higher-body densities in terms of one-body density matrices (KO 79a) . The form of the TDHF equations is important for several reasons . First, they lead directly to several conservation laws which are not only desirable for a physically plausible theory, but will serve as checks for numerical computations. Since TDHF is a mean-field theory with each single-particle wavefunction evolving in time through the same hermitian Hamiltonian, h ; the overlap matrix of the 'f/;'s, <1/1; l"\bt>. is time independent (normally 6ij ) . (1.11) Therefore, the norm of the TDHF wavefunction and the expectation value of the number operator are also time independent . The variational formulation of TDHF is useful for demonstrating other conservation laws. For any operator, 0, with no intrinsic time dependence, the time evolution of its expectation value in the TDHF wavefunction. ir, may be written as (1.12) where <0> is viewed as a functional of ~1/1;~ and ~t;l Use of the evolution equation,(1.8), results in iti.L<O> = ~r 6<0> 6<H> at /o. °'Y'; (ex) 6'f//( a) 6<0> 6<H> ~1/1/(a) 61/I; (a) (1.13) the analogue of the Poisson bracket. Using 0 = H immediately results in conservation of energy. It can also be shown th.at 111.:t <0> =0 for any one-body -7 - operaLor, if the commutator. [H, O] = 0 (see Appendix Al, Fart rY) . Tb.is, fo~·· uny reasonable choice of H, that is a Gaililean-invariant interaction, the tot21 momentum and angL.ar m:.:men~U:.":l. ¥:El be com:crved by the TDIJF sobti::21. The Form of the TDHi' equations will also ca.use problems i::i foe inte!F'etc.i.tion of the TDHF results (BO 76a). Since the mean field depends on the co:r.plste set of single-particle wavefunctions, ~1/li l; the TDHF equations will be coupied, nonlinear equations. Thus, if we tried to imagine two outgoing channels, they would be coupled even asymptotically. This arises because of the demand that the wavefunction be a single Slater determinant and not a sum over various possibilities. Thus, projecting the final state onto any final eigenstates of the system gives amplitudes which are time dependent (indicating transitions) even asymptotically. These "spurious cross-channel correlations" are the result of trying to extract information from the wavefunction which was already assumed not to be in the wavefunction. Since we can't calculate amplitudes for various final states, but merely some "average" final state, TDHF is a semiclassical theory . An alternative way of arriving at the same conclusion is to notice that a Slater determinant is a wavepacket (FL 78) . Thus, for example, a static HartreeFock (HF) solution with its localized center of mass represents a wavepacket of momentum eigenstates. When one of these Slater determinants is evolved with the TDHF equations, the coupling between the different eigenstates doesn't allow the wavepacket to spread. Thus, for the example above, the TDHF evolution of the static HF solution will merely involve multiplication of these solutions by a phase exp(-i2,;ei t In} . The wave packet of momentum eigenstates can't spread j because we aren't free to arbitrarily change the center of mass wavefunction . The TDHF (or HF) equations will simultaneously give the wavefunctions for the internal degrees of freedom and for the center of mass. This non-spreading of the initial wavepacket is a sure signal that TDHF is a semiclassical theory. -8- The moral of this discussion is that TDHF theory is a deterministic theory, calculating some kind of average evolution of our initial wavepacket: and that this wavepacket is only useful for calculating expectation values of few-body operators. Thus, our primary information will be simply the density as a function of time, and our principal results will be the number of fragments and their sizes, masses, charges, and motion. As we shall see later, the rich collision phenomenology exhibited by real nuclei will make the comparison of such gross, inclusive measurements with the theoetical results quite interesting . -9Section 2: 1be Effective Interaction For the internucleon nuclear interaction we use a phenomenological force familiar from static Hartree-Fock and RPA calculations, the Skyrme potential (SK 56), with a slight modification described later. Although in principle we might like to use a potential or G matrix more directly connected to nucleonnucleon scattering (DAV 74), the computations would be far too cumbersome . In fact, even the most direct methods require some minor phenomenological adjustment to agree with the data (NE 72) . This isn't too suprising since the choice of a potential to fit the scattering data is ambiguous; the calculation of the many-body effects in a finite nucleus to determine the mean field is very difficult for even a two-body potential, requiring some approximations; and the effect of higher-body correlations is completely undetermined. It was found that an effective two- and three-body interaction could give very accurate results for both binding energies and densities of nuclei (BEI 75) . The Skyrrne force is used because it is particularly simple to calculate. The potential energy is taken to be (YAU 72, EN 75) (2.1) This three-body force gives a density dependence, which aids saturation of the nuclear force. The Skyrme form parametrizes the potential by an expansion in momentum space. & I y<2l lk'> = t(J(1+xoPa)+*t1(k 2+k'2 )+t2k ·k'+iW0 (a1+u2)·kxk' , (2.2) where t 0 ,t 1 .t 2 , Wll and Xe are constants, the Skyrme paramaters; Pa is the spinexchange operator, Pu= *(1+u 1 ·«f2); a1 and a2 are the spin operators; and k and k' are the relative momenta of the two particles. Notice that in general & I v< 2> lk'> ¢. v (k-k'), so the transformation to coordinate space does not give a local potential, <rl v<2>ir'> = 6(r-l')v(r).The transformation to coordinate space, - 10 - gives <r I v<2 ) If'>= to(1+x0Pa)o(r)o(r ')-72t1H\72o(r)]o(f ')+[v2 o(f')]o (f )l +t 2 [ 9o(r)] · [Vo(r ') ]+i W0(cr 1+a 2) · [Vo (f) ]x[~o (r ') J where r , r' are the relative coordinates of the two nucleons . (2.4) Thus we can anticipate that the simplicity of the Skyrme form will be dw.e to the expansion of the range and nonlocality of the potential as differential operators . The t 0 term can be rewritten t 0 (1 +xoPa)o(r--T')o(r), so it is a local, zerorange force . Since this term acts only ins waves ; it has strength t 0 ( 1-xc) for like nucleons (singlet spin states) , while it has strength t 0 (1+x0 ) for unlike nuc leorn ( triplet spin states) . The t 1 term acts in s-d states and the t 2 term ac ts in p states; W0 is the strength of the spin-orbit force . In all of our applications we -will be investigating bulk nuclear dynamics . Since these spin-orbit forces c ause single-particle effects which are unlikely to significantly alter the bulk dync.mics, we assume a spin-saturated system and neglect the spin-orbit force . Therefore , each orbital is occupied by two protons or two neutrons , one spin up a...11d the other spin down. The three-body interaction contains only the first term of an expansion like the above . In coordinate space, it is (2.5) a zero-range. three-body interaction. Since the potential between two nucle ons depends on the presence of another nucleon. this term introduces a local density dependence. In Appendix A2. the expectation values of these potentials are derivc;,d for a spin-saturated Slater determinant, showing that the three-body potential above is equivalent to a two-body potential (YAU 72) (2.6) - 11 - The result, from Appendix A2, for the nuclear energy density is (EN 75) H.(r) = :n + + or+~o[*(1-xo)(pi+p~)+(2+xo)PJ,Pn]+~p+Pn)PpPn (t1+3t2) B 3(t2-t1) 32 -t2 72 (pp 'Tp-Jp +Pn 'Tn -Jn]+ [ppV2pp+PnV2Pn]+ (t1+t2) 4 (t2-3t1) 16 1' (2.7) 7 (pp Tn +pn Tp-2}p ' }n] [ppV2Pn+PnV2Pp] . The particle, kinetic energy and current densities for each isospin species q ( p for protons, n for neutrons ) are defined as: Pcz(r) = l; l-¥i,{r,q,) 12 • (2.8) '"'"" fli=q 'Tq(r} = l; IV-¥i,(r,q,)1 2 • ,,=q '"'"" and hfr) = l; Im[vt(r,q,)Vl',(r,q,)J , '"fi=q 1UJt where the Vt. are the occupied single-particle orbitals and the sums are over all these orbitals with isospin q. Densities without isospin subscripts refer to total densities; p=pp+Pn· T=Tp+Tn , and J=Jp+'fn . ( For notational convenience the ti.me dependence of the above orbitals and densities has not been explicitly included.) For light systems, where the Coulomb potential is weak, which have equal numbers of protons and neutrons an assumption of isospin degeneracy seems a reasonable approximation Since this assumption considerably simplifies calculations, it was used in many of the early calculations (BO 76a, BO 76b, BO 78, DAV 78a, DAV 78b, FL 78, KO 76, KO 77) and we use it for some of our calculations. With this assumption, each spatial orbital is effectively occupied by an alpha particle. Setting each of the proton and neutron densities to one-half the total density gives R (r)=l!_r+~o,02 +~&0 9 + (3t1+5t2) (p,--'P}+ (5t2-9t1) pv2p ( 2 .9 ) 64 The p'T-j..,_ term gives rise to a non-locality ( effective mass ) in the mean field. • 2m B 16 16 - 12 - The pV2p term, which is sometimes rewritten after integration by parts as -I Vp 12 • can be identified with the surface energy since the density varies rapidly only at the surface. The determination of the five Skyrme parameters t 0 .x0 .t 1 .t 2 ,andt 3 is somewhat ambiguous and has led to a variety of Skyrme potentials. Of course, a saturating internucleon force demands t 0 negative and ts positive. Four of the conditions used to determine the coefficients are the nuclear volume, surface and symmetry energies (as determined from the semi-empirical mass formula) and the saturation density. One linear combination is then undetermined, essentially the balance between the effective mass and the density dependence in producing saturation. A more refined approach uses these values as esti- mates and adjusts the parameters slightly to give better fits to the binding energies and charge radii of magic nuclei when used in static HF calculations (BEI 75). All of these potentials give accurate fits to both the binding energies and charge radii of magic nuclei and reasonable densities when used in these static HF calculations. However, with very different effective masses the different potentials give different single-particle level densities near the Fermi surface (BEI 75) . The level density decreases with decreasing m• ( increasing t 3 ). This determines the last parameter. Presently the two favored forces are Skyrme II and III with effective masses in nuclear matter of .58 and .76 respectively. Another force, the BKN force (BO 76a), is sometimes used for isospin degenerate systems because of its simplicity. For these systems, there is only one nonlocal term, so with a choice of 3t 1 +5t 2 =0 one can select a purely local potential. This expression is approximately true for Skyrme VI with t 1=271.67 and t 2 =-13B.33 so a slight readjustment of these parameters yields a local mean-field potential. Therefore, use of this force, the BKN force, greatly simplifies calculations. However, for the nondegenerate case only the trivial choice t 1 =t 2 =0 can eliminate - 13 - both of the nonlocal terms . This choice is completely unreasonable since it also eliminates the surface energy terms t p'l2p ) . The modification of i.he 3kyrme form menl1oned earlier con.s1::;t::; u! replei.cing the~e surface energy terms by direct Yukawa interactions . (BO 76a, J~ E 78) There are two reasons for this change. First, the finite range of the nucle3.r fo r ce is more accurately approximated (HOO 77) . The expansion of the exchru."'lge term is the crucial simplification which makes the Skryme form so useful, but there is no technical reason to eliminate the convolution of the density with a finite range interaction occurring in the direct term. Second, when the calculations are done in coordinate space (as ours are) the evolution with V2 terms is unstable (BO 76a) . That is , the wavefunctions develop a periodicity equal to the number of points in the discretization formula for V2 . The connection between the surface energy terms and a Yukawa energy can be seen by a Taylor series expansion of the Yukawa energy. E n =4rra 5 ·=ffcirdr... exp(IT~· 11 a) (r) ·( ...) lr-r'j/a Pq Pq r (2.10) lfdrpq (r)p,,l'(r)+1fla f dr[pq (r)V'2pq·(r)+pq ·(r)V (r)]+ .. . ) 2 2 p9 Therefore, the surface energy terms are replaced by (2 .11) VL and Vr.r are the strengths of the interaction between like and unlike nucleons respectively and a is the range of the force . Since we choose a to approximate the range of the G matrix in nuclear matter (BO 76a), the selection of a set of Skyrme constants determines Vi, and Vu . The presence of the first term above requires adjusting the parameters t 0 and Xo of the zero-range component of the potential. Therefore, the parameters t 0 and Xo are replace d by t"'o and Xo satisfying 14· 1 ...., ...., 1 ttc(l-xo)= tto(1-xt;)--4m1 3 Vi .... ) (.-? •... 1::. 2 and ...., ...., 3 ~c(2+xo)=*to(2+xo)-4rra Vu + For the isospin degenerate case we need only one Yukawa interaction, Ey= :1 drdr'exp (l~:;·,~:)p(T)p(r') (2.13) The surf ace energy terms can be replaced by choosing (2 .14) and replacing t 0 by t~ satisfying 3"' ::::;....to-47Ta33 v :::;...to= 8 8 2 (2 .1 5) For the Coulomb energy we will use only the direct Coulomb contribution, (2 .16) TDHF studies of light ion systems have shown the inclusion of the exchange energy has negligible effect on experimental observables (KRl 78) . Studies of heavy ion systems have shown that the exchange energy is approximately 30 times smaller than the direct and that it fluctuates little throughout the collision (DAV ).++ Therefore, in summary, for the isospin non-degenerate systems we will use the energy functional, + Although the contribution of higher-order terms causes small changes in the results of static HF calculations (HOO 77, BO 76a), no ane has attempted to readjust the parameters to maintain the same accuracy of fit to energies and radii for each set of Skyrme paramete::-s. This adjustment would require a whole new set of static HF calculations, and the cha.IJ8es in the parameters would be very small compared to the differences for different Skyrme forces . ++ All these TDHF studies e.pproximate the exchange energy as the first term (Slater term) of a density-matrix expansion (NE 72), Hc(ezchange )=- ~ J1 ~e2JrG[pP(~lt-. since an exact treatment would destroy the mmplicity of the Skyrme form. Previous static HF calculations (TI 74) showed that this is accurate to better than 103 and the discrepancy varies slowly with nucleon number. These studies also showed that the excha.IJ8e contribut :c::1 is less than 103 of the direct contribution. - 15 - +f ar'Oy( jr-T' l J)[ ~[pp (r)pp (r')+pn(r)pn(r')]+ V!!Pp (r)pn(r ') +*fdr'Oc( lr-T' l)pp(r)pp(r')] . _ exp(-jr-T'lla) _ e2 .... I and Oc- .. .. I . The coefficients are given by 1r-T' I a 1r-T' where Oy- 1 VL=bt2-t1)--5 1 V!l=-A-<t2-3t1)--5 32 2rra 16 with parameters chosen from Skyrme III . 2na For the isospin degenerate systems we will use the energy functional, <H>=f dr[! ,-+ ~~2 + 116 tsa +:f ctr'Oy( lr-T' l)p(r)p(r') 3 + . (2.18) ~ f ctr' Oc( 1-r--r· l)p(r)p(r')] . 'Ib.e coefficients are given by v-( 5t 2-9t 1 ) 1fa 1 2 64 4 "' 5 8 and t 0 =t 0 -'34na 3 v . 2 with parameters of the BKN force. The values of both sets of parameters are shown in Table 1. - 16 - Table 1 Skyrme Parameters Force to t1 MeV-fm3 n Ill BKN MeV-fm 586.6 395.0 250.0 -1169.9 -1128.B -1089.0 t~ t2 5 MeV-fm -27.1 -95 .0 -150.0 5 Yn MeV-fm6 9331.1 14000.0 17270.0 0.34 0.45 - Nuclear Matter Results EIA kr K MeV -16.00 -15.87 -15.77 fm- 1 MeV 342 356 368 Force II Ill BKN 1.30 1.29 1.29 m" .... 0.58 0.76 1.0 Derived with r~e a=0 .45979 fn Force '.Vn MeV-fm3 II III BKN -104.49 -334.52 -497.73 4.01 1.74 - Vr1 V,, MeV MeV -868.53 -444.65 -619.60 -355.79 -363.04 We derive the mean field, hv, to use in the evolution of the single-particle wavefunctions, ih:t '/J;.q =hv'l/l;.q, by the variation of the energy functional, o~>) ="v1/11q. As before q labels the isospin and i the orbital. In the isospin o'l{liq r nondegenerate case this mean field will, of course, be difierent for protons and neutrons. Terms in <H> involving the variation of the density will be very sim- (2.19) Variation of the kinetic energy or current densities will require integration by parts. Th.us, with,.,= IV'l/lill 12 , ch·, =~[vo1.; · v'l/I" +V'l/J~·V6'1/Jiq] (2 .20) ( Therefore, integration by parts gives (2.21) - 17 - .. where the operator V acts on all functions to its right. Similarly for (2.22) Therefore, J .. , . . , d;q.(r') 1 ... ... .. ... ... ... ... ... d.T F(r)- 6lt~(r) c:6qq·i.f£F(r)-V'¢'ilz(r)+V·F(r)'¢'iiz(r)] . (2.23) Combining these formulas and defining q'rtq gives the result (EN 75, NE 78) hq(r)=- ~ V2 +ToLPq+Tor.JPg·+Ts(p:·+2pqPg·) (2.24) +CL ( T q-Vpq V+i rlq ·V+VJq]) +Cv(Tq·-V·pq•V+i[]q .. v+v.;q.]) + °VL UYq + Vr;r Uyq·+ Uc6qp . where the Yukawa potential is uYq= f dr'Or( ir-r· l)Pq(r') • (2 .25) and the Coulomb potential is (2.26) and the differential operators V act on all functions to the right including the wavefunction. We can see now that the nonlocality of the mean field, contained in the terms CL and Cv. has been expressed as a second-order differential operator just like the kinetic energy. By rearranging the terms slightly, we can define an effective mass, ~ ~ ~(r) - 2m I CLpg+CvPq• (2.27) which will appear in the mean field as -V· ~(r) V. For the isospin degenerate case, with the BKN force, the mean field is quite simple and can be written immediately (BO 76a). - 18 - n 2 3..... 3 h(r)=-~ +~()0+~ 31> 2 +VUy+*Uc m 4 1 6 , (2.28) where f dr'Oy( ir-r· J)p(r') (2 .29) f ar'Oc(jr-r·1)~ . (2 .30) Uy= and Uc= Therefore, this BKN force results in a local mean-field potential. Although a local two-body interaction with finite range will give rise to a nonlocal Hartree-Fock potential (as shown in Apprendix Al, Part III), we see that we can choose a nonlocal two-body interaction which will result in a local Hartree-Fock potential for spin-saturated systems. - 19 - Section 3: Brief Review of Heavy-Ion Reactions (For a recent review see (LE 78) . ) The early heavy-ion beams of about 20 years ago were confined primarily to 7Li, i 2 c, 14N, and 160 at energies less than 10 MeV /nucli::on. Scattering sL.:.die s were done to probe the optical potential between nuclei, and resonances revealed the occurrence of transient nuclear molecules (BR 60). Compou_"ld nucleus formation and decay was used to study the behavior of highly excited nuclei with large angular momentum. Then, about 10 years ago, with the use of heavier beams, a new phenomenon called deep-inelastic reactions was discovered . Two fragments similar to the target and projecWe are emitted, but with pronounced damping of the initial kinetic energy (LE 73) . In the most dramatic cases, there can even be complete energy relaxation; the two fragments are emitted with energies corresponding to Coulomb repulsion, just like a fissioning compound nucleus . However, the masses of the fragments are completely different from the expectations of a fissioning compound nucleus . The lack of a fissioning compound nucleus was convincingly demonstrated when reactions leading to nearly the same composite system had different final fragments (TA 75) . Thus, the energy damping and mass exchange show that a composite system is formed, but the entrance channel dependence shows that the system splits again before complete equilibration can be reached. These reactions, therefore, allow the study of interaction times and equilibration times of dit!erent degrees of freedom for a tremendous variety of targets, projectiles, and energies. This opens a new field of collective phenomenology correlating the scattering angle, energy loss, mass of the products, and the charge of the products . We can then develop models using parameters such as interaction potentials, viscosity coefficients. friction tensors, and diffusion coefficients (NO 76) . - 20Or, in our case, we try to describe all these collective phenomena from a microscopic approximation using a previously determined nucleon-nucleon interaction. The application of classical concepts, like that of an ion-ion trajectory, to the problem is justified since the wavelength associated with the relative motion of the nuclei is small. This wavelength is (NO 76 p.3) (3.1) where v is the relative velocity, µ is the reduced mass, ~ is in fm, Ecm is the center of mass kinetic energy in MeV, and A 1, A2 are the atomic mass numbers of the target, projectile. Thus, for some of the cases we will be describing later, this wavelength can range from .2 fm, for 160+ 160 at Ecm =55 MeV, to .004 fm, for 86 Kr+ 139La at Ecm =312 MeV. Also, the Sommerfeld parameter, 1), is large for these collisions (NO 76 p.19). (3.2) Here dmm is the distance of closest approach on a Coulomb trajectory; and Z 1 , Z 2 are the charges of the two nuclei. This parameter ranges from about 5, for our 0+0 case, to about 400, for our Kr+La case. This means the scattering can be described with a classical defiection function, and that classical Rutherford scattering is expected when the minimum distance of approach lies significantly outside the sum of the radii of the two nuclei. Therefore, at large impact parameters Rutherford scattering is expected. At small.er impact parameters, the attractive nuclear interaction bends these trajectories forward in angle (WI 73). For still smaller impact parameters, when the nuclei touch, the damping of kinetic energy begirul along with mass, charge, and angular momentum transfer. Also, the temporary sticking and rotation of the system bends the trajectories further forward in angle. At still smaller - 21 - impact parameters, with significant overlap of the two nuclei, a composit '=' system is formed; in some cases, a compound nucleus . For the formauon of !:I. compound nucleus the target C1.nd pwjeetile mu::;L completely lose their individltul characteristics, and the system must eY h~ ·..::it equipartition in the occupation of all accessible degrees of freedom. Tr~::n, ':.!'le decay process will be independent of the formation process, except for s:ich restrictions as the conservation of energy and angular momentum. The system is still considered fused into a compound nucleus if a few pre-equilibrium nucleons or alpha particles are emitted. For light nuclei fusion is quite clear because the composite system is stable against fission. Therefore, we observe an evaporation residue together with light particles: mainly protons and neutrons, '\\ith alpha particles if the system has high angular momentum. Thus, for light systems, with small Z 1 Z 2 , the complete-fusion cross section is a large part of the total reaction cross section; and only a small fraction of the reaction goes into the inelastic, grazing collisions. For heavier systems the situation is much more complex. These systems can be unstable and fission - either because they are simply too large, or because they have too much .excitation energy and angular momentum - so the complete-fusion cross section goes to zero. Evidence of fusion-fission behavior is the complete relaxation of the kinetic energy of the system, and mass and charge distributions of the fragments which are independent of the initial fragments . However, these same factors, which would cause the compound nucleus to fission, will inhibit its formation; causing the system to break apart before becoming completely equilibrated. Thus, for heavier systems even the fusionfission cross section decreases, and the deep-inelastic reactions become more important. For later sections, it is useful to discuss in general terms some of the final- - 22 fragment correlatioru; we will be looking at. and the type of information we can learn. For instance, as mentioned earlier, the correlation between energy arid mass led to the disc0vcry of deep-inelastic reactions. Thu _·; , con t.Jurs of lll·2 nwnber of events versus the total kilJ.etic energy and mass of the fragm en '. s show two schematic possibilities. (Fig . la and lb). both of which have been observed (TA 75) . Both of these reactions show quasi-elastic peaks at the initial masses and energy. Fig . la shows evidence of fusion-fission since the energyrelaxed products are centered about the average mass. Fig . lb shows evidence of the deep-inelastic reaction since the energy-relaxed products still maintain nearly their original identities . Thus, using the previously discussed classical concepts, we would say that for the system in Fig . la a certain range of impact parameters have interaction times longer than the equilibration times for both mass and energy. For the system in Fig . lb the energy equilibration is much more rapid than the mass equilibration, and a certain range of impact parameters have interaction times which are between these two equilibration times . No impact parameters have interaction times as long as the equilibration time for the mass . By looking in more detail at this correlation, we can relate the amount of energy damping, mass transfer, and the spread of these two distributions (HU 76). Thus, for example, a certain increase in energy damping corresponds to a certain increase in mass transfer; but the time scale, so far, is completely undetermined. Looking for correlations with scattering angle proves very useful because this allows the estimation of an interaction time (SC 77); and thus, sets a definite time scale for the equilibration processes. For a given angular momentum (or impact paramet er) the deviation from Coulomb scattering is a measure of the rotation of the composite system. Therefore. a knowledge of the angular momentum and moment-of-inertia of the composite system gives an estimate of - 23 - :~he interaction time, TtTtt = .~ /).'() ; where 1 is the moment-of-inertia, L the & · IS '-liar momentum. and 6,1') the -tev:iation from Coulomb scattering . Typical int c: ~·ac­ tion times range from 10-22 to several times 10- 21 sec . (SC 77) . Thus, for instance, a gradual energy loss as interaction time increases would lead to a classical trajectory in E vs . 19. like that in Fig . 2a (WI 73); exc e pl , of course, that we can't differentiate between positive and negative angle scattering, so the trajectory folds back at 19.=0° . Since each point in the trajectory corresponds to a particular angular momentum, we can also plot E vs . l and 1' vs. l. From the plot of 19. vs. l, we can derive an interaction time, Ttnt vs . L. These are also shown in Figs . 2b, 2c, and 2d for the same system as in Fig . 2a. In a realistic case, showing spread in both E and f), the trajectory ir: E vs . 19 becomes a ridge in contours of the cross section. The example in Fig . 2 is fai:ly typical of moderately light systems (GA 76), but the variation of both interaction times and equilibration rates causes the shape of these contours to vary for different systems . For heavier systems, the trajectory drops more steeply in energy (SC 77). In summary, we can make a few general comments about these heavy-ion reactions. First, the widths of the mass and energy distributions increase as the mass transfer or energy damping increases (HU 76). This shows that these deep-inelastic reactions are relaxation phenomena. Second, the kinetic energy equilibrates much more rapidly than the masses of the fragments (LE 73). Third, looking at different isotopes shows that the neutron excess, or ratio N /Z, equ:'.librates even more rapidly than energy (TN/ z~10-22 sec.)(GA 75, KRA 77) . And fourth, the large variation in all these rates and interaction times for different systems leads to quite difierent collective behavior, and presents a formidable challenge to a unified, microscopic theory. - 24 Section 4: Application of TDHF to Heavy-Ion Reactions The application of th~ TDHF method to heavy-ion collisions requires large scale numerical calculations. We write the TDHF equations in coordinal~ space and use finite difference methods to follow the evolution of the single-particle orbitals. These orbitals are described by their values on a uniformly spaced Cartesian mesh within a rectangular box. Vanishing boundary conditions are imposed outside this box. We use a coordinate space representation rather than an expansion in a set of basis functions as in static calculations because it more uniformly describes the wide variety of shapes which occur and because the matrices are sparse. Our discrete space-time mesh has a typical spacing, 6.x, of approximately 1 frn and a time step, 6.t, of approximately 5 x10-24 sec . Therefore, two separated nuclei require a spatial box approximately 30 x 30 x 16 (the smaller dimension is perpendicular to the scattering plane) and the evolution to the final state requires on the order of 100 time steps . For our heaviest nuclei , with approximately 100 wave functions, the computational effort is substantial. Storage requirements exceed that available on a CDC 7600, requiring tapes and substantial input/output; the computing time can be as high as 5 hours for a single impact parameter. F-ortunately, the Cray we used most recently can meet the storage requirements and also decrease the computing time to about one hour per run In our calculations, we impose a reflection symmetry with respect to the reaction plane, thus cutting our numerical work in half. For collisions of identical nuclei, we also impose a point reflection symmetry through the center of mass of the total system, cutting our numerical work by another factor of t wo. Since the TDHF equations, ih:t '/ij =hV;; .. are first order in time, they can be solved formally by using a time evolution operator, (4.1) -25where U is the unitary operator U(t,t 0 )=Texp[ .' -;;.J dTh(T)] (4.2) lo and Tis the time-ordering operator. Therefore, we only need to know the set of single-particle waverunctions at one time to solve the problem. Our discrete approximation to this operator, U(t +At ,t ), is discussed in Appendix A6. The calculation or h'l/I;, required for the evolution, is considerably simplified because the only nonlocal terms are differential operators. Therefore, the discretization, as discussed in Appendices A3 and A4, results in a sparse matrix for h with nonzero elements only near the diagonal. Our initial condition is an impending collision between two nuclei. At a large separation. when nuclear forces between nuclei and Coulomb induced distort.ion are negligible, the system is completely described by the nuclei. their energies (kinetic plus Coulomb) and their impact parameter (or angular momentum) . 'We choose two static Hartree-Fock solutions at positions which result in a large separation without either nucleus getting too close to the walls We can generate a translating HF solution, '1', by of our spatial box. +=exp[ ik A]t. where t is the static HF solution, composed of single-particle wave.ruD.ction.s <fJ; and R= ~ i 1 is the center of mass coordinate of the nucelus. (A is the mass number of the nucleus .) The nucleus then bas a total momentum .. The existence of these translating solutions is intimately connected with NC the conservation or momentum. The time-dependent single-particle wavelfunctions of i1 are -;1 where where iJ (4.3) =:... .. c: ~ and m is the nucleon mass. We do our calculations in the ~ .. . center-of-mass frame, so we choose A 1=-K2 =K and our initial state becomes - 26 - +=exp[ikR 1]1f> 1exp[-ikR2]lf>2 . (4.4) two translating static HF solutions. The total kinetic energy is given by TKE= (M) 2 2µ (4.5) I A1A2 with reduced mass µ=m (Ai+A ). The total energy is 2 Z1Z2e 2 E=TKE+ I.. .. I Ri-R2 and the total angular momentum is given by I (4.6) l=11Kx(R1-R2) (4.7) (or equivalently the impact parameter is IR 1 -R2 Isin1' where 1' is the angle .. . . between Kand R 1-R 2 ) .+ The static HF solutions are calculated using the same energy functional as the TDHF solutions . (Otherwise they would not be the ground state of the nuclei in the collision) These solutions are obtained by an imaginary time step method (FL 78). We use an operator similar to the actual evolution operator, U(t ,t 0 ), but with t replaced by -i;. This results in an exponential decrease, exp( -~-r ~. of higher states relative to the ground state of the static solution. Since this operator is nonunitary, we e,1.so renormalize the wavefunctions after each time step. This procedure enables us to use a code almost identical to our TDHF code and results in very rapid convergence to the ground state static HF solution. For the exact solution, we can easily show that the imaginary time step method will result in the lowest energy state which has any overlap with the initial trial state. Any+ can be expanded +(O)=l:°'i~i where the ti are the complete set of " eigenstates Ht,=£'-&~..:. The exact evolution of+ is given by +(t}=l;°'Lexp[ ~t ]~, ' + Difierent diacretize.tion errors for the different wavefunctiom result in slight]y different boosts. Since an we.vefunctions are coupled by the mee.n fteld, this results in 6 small initial trmisrer of collective ldnetic energy to internal e:zcitation. (4.B) - 27 - where wi =E-~ . Replacing t by -i; gives '1:,( -i: ) =exp( -WoT n:Cli exp [ - (Ult -Wo)T N'i (4 .S) ' Renormalizing+ as a function of time gives (4.1 0) and (4.11) For the HF solution the method is very similar, 'T +(-ir)=N(1)Texp[-~ { d1'h(i1')]+(0) , ( 4.12) but the dependence of h on the wavefunction +means that a tiny overlap of the intial guess 'lf(O) with the ground state, IP 0 , is no longer a sufficient condition for convergence to the ground state. If the initial guess is bad enough so the mean field h (0) looks nothing like the mean field in the ground state, h 0 , it is possible that the self-consistency will result in an excited state. similar to '11(0) , 1~ith a mean field similar to h(O) . In practice, this is not a problem because we always know the approximate ground state. The resulting mean field is always close enough to h 0 so that+ evolves closer to IP 0 and the convergence is very rapid. This initial state is evolved by our discrete approximation to the evolution operator (in both space and time) until some final time. The complete set of wavefunctions ~'¢-;(t)J contains far too much information, as we are primarily interested in the proton and neutron densities . We observe one of two possibilities, either the density separates into two fragments again rather quickly or stays as one fragment long enough to be considered fused . This fusion criterion can be improved by considering the mean radius of the system (or a collective coordinate describing fragment separation as discussed later) although in some cases the determination can be ambiguous. - 28 - As an example, we will jump ahead to sorr;e data which are cliscusseJ more <r;;!> fully later. The root mean square radius of the system, ( ~)*. for four .:::.ngular momenta is plotted in Fig. 3. All show a rapid decrease initially, as the two nuclei approach one another, followed by a slower increase. For the pt:ripi.1eral angular momentum, L= 165, this increase is very steady with th~ d~nsit:," separating into two fragments at about r=9.16 fm and the fragments moving apart at a speed somewhat slower than the initial speed. The curve for L = 80 is a typical deep-inelastic event. The rms radius can oscillate several times before beginning a steady increase leading to separation into two fragment~ at about r = 8. 9 frn.. The lowest curve for, L = 50, is an unambiguous fusion event . The rms radius shows oscillations, indicative of the excitation of the system, but no net increase . The curve for 1=5 is one of the most ambiguous cases since tr. . e rms radius shows a slow but definite increase . Extrapolation of this increase would result in the system reaching a typical separation radius of 9 fm at a time approximately 160 x 10-22 sec, far beyond our calculations . However, in any sytem interacting this long, the two-body collisions will become important . The excitation of these additional internal modes will damp the energy in collective modes and result in fusion. Therefore, events such as this are considered fused and we can use plots such as these to determine whether or not to continue our calculation a little farther to look for separation. For the angular momenta which lead to fusion, we can learn little except that the system fuses . The behavior of the system after fusion is better described by traditional statistical models (NO 76 p.224). We can calculate a fusion cross section, aF, by a sharp cutoff formula C1F= 1r2 ~)2l+1) k , (4.13) '1 where hk is the relative momentum of the nuclei and the sum is over all l which - 29- fuse, lF. When fusion occurs, it occurs for a range of l from l< to l» Then we can do the sum, (4.14) The classical correspondence is most obvious from the first formula. Since the total angular momentum squared L 2 =1f!L (l +1)=(hkb ) 2 , where b is the classical impact parameter, substituting above gives uF=7r(b~-bn. the area between these impact parameters . For lighter nuclei at low energies, fusion dominates and the fusion cross section is the main result. For heavier nuclei, scattering becomes more important. For the angular momenta which separate, we get much more information, primarily from the densities. By integrating the densities over each fragment separately, we can calculate their mass numbers (A 1, A2). charge numbers (Z,, Z 2) and centers of mass (R 1 • R2 ) at each time. (After separation the masses and charges are very nearly constant.) The motion of the centers of mass determines the important dynamic quantities in the center of mass frame: A 1A2 -+ -+ 2 TKE=~(Ai+A )(R1-R2) . 2 (4.15) Z1Z2e 2 E=TKE+ I.. .. I ,and R1-R2 L=mA1R1x(R1-R2)=m.A2A2x(R2-R1) These quantities are actually estimated at the half time step, t -~t, by .. R,(t)-il,.(t-M) Rt=----At (4.16) and i4= ~(t)+ACt-At) 2 We also calculate the orientation of the principal axis of inertia and various moments of the density distribution for both the total system and the fragments . In addition. we define fragments even when they aren't separated, as - 30 - discussed in Appendix A7. Although this procedure is somewhat arbitrary when the fragments overlap greatly. it determines possibly useful collective coordinates throughout the evolution. To calculate the scattering angle, we need to translate our initial nuclei back to infinity and our final nuclei forward to infinity along Rutherford trajectories (GO 50). Therefore, the change in angle before the calculation begins is fl1'0 =cos- 1 [ .::-4 _1___,,..1)]-cos- - l ~ • ei ei 1[ ~ri (4 .17) and the change after the calculation ends is where (4.18) Z 1• Z 2 • A 1 , and A2 are all different. (Of course A 1+A2 =Ar and Z 1+Z2 =Zr and the total mass and charge are not changed.) Therefore, the scattering angle is given by (4.19) The ambiguity in the sign of fl1' due to the multiple values of cos- 1z is resolved by choosing cos- 1z between 00 and 1800. This results in fl1'>0 as we would expect for a repulsive potential. For scattered events. we also define two moments, contact and separation, from the pictures of the density projected onto the scattering plane. Thus, we consider the fragments to be separated when a line in the scattering plan dividing the density distribution has a maximum projected density, Ppr= J d.zp, less than Pc:r=.lfm-2 . (1bis would correspond to a neck width of .6 fm at - 31 - P•at =.17fm-5 .) Defining these events allows us to approximate the interaction Ume by 7'tnt=t (separation)-t (contact} and allows comparison of other quantities, like the rms radius of the system, at contact and separation. The last results we obtain require more information from the wavefunctions than the total proton and neutron density. These are the single-particle ener- gies and the mass and charge widths of the fragments. The single-particle energies, (4.20) e;=<'t/I; lh l't/I;> are obtained quite simply since h 11;> is calculated for the evolution of 1/1; (see Appendix A6, Part N). The mass and charge widths can be obtained by defining a number operator for particles in a particular fragment (BO 76a) (4.21) where ir+, ir are creation, annihilation operators for a particle of position r and r' is a Heaviside function (MA 70) so the integral extends only over one fragment (as defined by the split in Appendix A7). <NR>= J drtJ(r)~'jt,(r)'Jlt(r) (4.22) ( =~w.a i The dispersion of this number operator is given by (4.23) The first term gives <.FJR'>=J cfrtJ(r)++(r)ir(r}+ J drdr'tJ(r)~{r')++(r)lf!+(r')ir(r')'l'(r) (4.24) =~Wg; +~(UJ.:'WJ;~wji) ( 'i Since (4.25) (4.26) Thus for a Gaussian distribution, we can define a full-width half-maximum for - 33 - charge or mass by (4 .27) ,,·here the sums over ij include only prntons ft:.r the cb.c.,,rge w1::lth anC. bc':.h is2spins for the mass widths . We note that this width is inherently limited by the one-body nature of TDHF (DAS 79). (Since this r is a two-body operator, the TDHF results will be less accurate.) Performing a unitary transformation to diagonalize W\i gives O'~R= ~Wu (1-wu) (4.28) i so with ~wii =AR i (4 .29) and the absolute maximum is 2 ( )-A aRRma.x-4. (4.30) - 33 Section 5: Brief Review of Earlier Calculations The earliest investigation of TDHF for nuclear collisions was a onedimensional calculation of the dynamics of slabs of nuclear matter, infinite and homogeneous in two directions but finite in the other direction (BO 76a). These authors introduced the BKN force and ignored Coulomb forces . These slabs showed nearly constant saturation density in the interior with surf aces about 2 fm thick. Single slab dynamics were studied, both breathing mode oscillations and the interaction with an external barrier. However, for our purposes their calculations of slab collisions are more important. These revealed a wide range of phenomena - fusion, resonances, deep inelastic scattering and fragmentation as the bombarding energy was increased. At the lowest energies, up to 1 MeV per nucleon. the nuclei fused. The coherent translational energy went into internal excitation, that is, oscillations of the system. At higher energies, these oscillations became more vigorous and could lead to scission. In this energy range, from 1 - 2 MeV I A, a delicate interplay between single particle effects (orbitals bouncing between the walls) and collective effects (scission of the walls) caused resonances, sometimes fusion, sometimes inelastic scattering . From 2 - 15 MeV I A, there was deep-inelastic scattering where the fragments separated with reduced relative kinetic energies (sometimes only 10% of the initial) and internal excitation. At still higher energies fragmentation was ob served with 3 or 4 outgoing fragments . These results encouraged the first realistic calculations: two dimensions with axial symmetry . The se calculations were first done for head-on collisions (BO 76b, KO 76) but were very soon extended to any impact parameter by the use of a r otating frame approximation (KO 77) . 'Ibis approximation reduced 'I'DHF to two dimensions (r, z) by requiring axial symmetry about the rotating line between mass centers. It then added a rotational energy to the Hamiltonian - 34- £2 H'=H+ 2 /(p) (5.1) where L was the conserved relative orbital momentum and I was the moment of inertia, a prescribed functional of the density. This added a centrifugal term to the HF potential W'= W-*'->2 OJ 6p - ~ d.1' - T L where "'- (5.2) At large separation I was taken to be that of two point masses, and, when the nuclei coalesced, it was taken to be that of a rigid rotor. The wavefunctions were evolved with W' while simultaneously e.> was integrated to give the orientation of the symmetry axis in space, 1'. This approximation should be best for peripheral or head-on collisions and at lower energies. These first calculations for symmetric collisions of 160 or 4tlca showed inelastic scattering and fusion similar to the experimental data but with too little dissipation and too little fusion. Since these deficiencies could be due to either TDHF itself, the effective interaction, or the two dimensional approximation, it was important to do full three-dimensional calculations. Calculations were done for the same systems, yielding fusion cross sections· at various energies (BO 78) and detailed scattering results at one energy (FL 78) . The results at large and small impact parameters agreed with the axially symmetric results but significant differences appeared at intermediate impact parameters, where nonaxial deformations were large. In some cases, the axial asymmetry, <z2-y2> , was as large as 0.2. The 3D calcula<z 2 +y 2> tions provided a microscopic realization of tangential dissipation which decreased both the angular momentum and energy in collective modes, particularly at higher energies, and led to increased fusion. The scattered solutions for 180 at one energy. showed a behavior typical of our later calculations. The highest angular momenta showed no energy loss and - 35Coulomb scattering . Smaller angular momenta led to decreasing scattering angle (becoming negative when nuclear attraction dominated Coulomb repulsion) and increasing energy loss until the system fused . This behavior is to be expected (see section 3) although the details can be compared with experiment. 'Ib.e unexpected result was that at the lowest angular momentum there were also scattered solutions with nearly constant energy loss (about 35 MeV of the original 52.5 MeV) and little loss of angular momentum. The fusion calculation for 160 showed that this low-l "cutoff" appeared at Elllb =54 MeV and persisted at higher energies. This lower angular momentum limit occured because too much energy was transferred to a single collective mode, in this case a vibrational mode . For TDHF calculations this mode is strongly excited because the long mean free path of the nucleons results in an "almost transparent" nucleus, but it has been shown to be a more general phenomenon (BR 78, TS 74, SI 76) occurring even in drops of water (AD 68) . Thus, the upper and lower angular momentum limits can be thought of as rotational and vibrational limits. Both the upper- and lower-l cutoffs were approximately proportional to ..JE-ET, though with difierent proportionality constants and energy thresholds, ET. This can be interpreted as a limitation of radial kinetic energy since E,.u=Er will result in l=[ ~ (E-Er)]*. For the upper-l cutoff this means the nuclei can climb in over the fusion barrier, for the low-l cutofi this means that a more compact configuration (smaller I) can climb out over the fission barrier. These 160 fusion cross sections fit the data very well up to Eia1> = 120 MeV, in contrast to optical model predictions which greatly overestimate the fusion cross section above about 60 MeV. (The potential fit elastic scattering for EllllJ:s:70 MeV.) The higher cutoff was in excellent agreement with the optical model up to 80 Mev and showed only a slight lowering above this energy due to - 36- the onset of peripheral deep inelastic processes. This discrepancy between the optical model and experiment was then resolved in TDHF by the existence of a low-l cutofi'. The fusion calculations for .40 Ca showed the same general behavior of the upper limit but they showed no lower cutoff. The energies were all below the threshold since a low-l cutoff was seen in the higher energy 2D calculations which should be exact for l =O. These cross sections were less accurate and showed some sensitivity to the interaction parameters. However, all of the forces were local forces so these calculations may give us more information about the effective interaction. The experimental search for this low-l cutoff is as yet unresolved and work is continuing. More will be said about this later for our specific calculations. The need to extend these calculations to more realistic forces (even if still the Skymre type) and heavier nuclei led to a renewed investigation of 2D approximations (DAV 78a, 79,80) since the 3D calculations were already saturating computer capabilities. Heavier nuclei are important both because the increase in deep inelastic events will yield information about equilibration processes in TDHF and because the classical concepts we are using will be more appropriate. Refinements to the axially symmetric approximation were made by changes in the moment of inertia (DAV 78b) . Several prescriptions were tried which permitted angular momentum to be transferred to internal excitation, thus allowing tangential dissipation. These results were a better approximation to 3D results than the previous 2D results and agreed quite well at lower energies. However, at higher energies, above the threshold of the low-l cutofi', there is no escaping the importance of nonaxial deformations for intermediate impact parameters. - 37Section 6: Separable Appro:ximation The deficiencies of axially symmetric approximations and the intractability of full 3D calculations prompt the investigation of other approximations which can limit the number of variables treated explicitly but still retain those degrees of freedom necessary for an accurate solution. Once such approximation, the separable approximation (DE 78a,78b), is tested in the first published paper included here (KO 78). The separable approximation reduces the calculation to two dimensions by forbidding motion perpendicular to the scattering plane . However, it imposes no restrictions on motion in the scattering plane and will allow the important features of nonaxial deformations and tangential dissipation to be calculated self-consistently. The freezing of the wavefunctions in this transverse direction is motivated by the mean-field dynamics which govern their evolution in the TDHF approximalion. The saturation of nuclear forces results in a nearly constant density and thus a nearly constant mean field in the nuclear interior. Therefore, the orbitals which are initially moving in the scattering plane will continue to move predominantly in the scattering plane during the collision. The separable approximation constrains the single-particle wavefunctions to the factorized form. (6.1) where z is the coordinate normal to the reaction plane and x ,y are coordinates in this plane. The functions X• are real and time independent so only the 2D complex wavefunctions rp, need to be evolved in time. The time evolution is obtained from the TDHF equation, ih:t '/Ji =h'l/I;. Since XJ is time independent, the equation becomes (6.2) - 38- The functions XJ are normalized so we can write the evolution equation as 'ill.:t 'I'; (z ,y)=(h ';'P; )(.:z: ,y) (6.3) (in general h and h' will be nonlocal) with the projected Hamiltonian given by h '(z ,:i:',y ,y')= J d.zd.z 'XJ(z)h(x ,z',y ,y',z ,z ')X,; (z') . (6.4) This equation can also be derived from a variational principle with. ~<H> 6 •( ) '/>; z .y ' (h J'Pi )(z ,y ) (6.5) (Variation of only rp;• would leave a XJ on the left-hand side of the Hamiltonian and give the same result as Eq. (6.4)) Notice that this projected Hamiltonian is, in general, different for different orbitals. Since the mean field (that is, the same field for all wavefunctions) is importa.nt for maintaining the orthogonality of the single-particle wavefunctions 1f;1 , we must check the orthogonality in the separable approximation. At the initial time <rp; Irp"><x; lx">=6i.1 can be satisfied with either <Xj lxt>=oiJ or <rp1 \rpi>=oii permitting a natural grouping of the 1/IJ · The orthogonality of the 1/1; is trivially preserved for wavefunctions with <x; IXi >=0 since these functions are independent of time. However, maintaining the orthogonality of the 'P; for the other wavefunctlons with <x.; lx,>~0 is not trivial since ~'Pi jrp,>=<rp; lht-h; I'Pt>. Orthogonality requires ~ =h; and Xi =XJ if <XJ lx,>iJfO. Therefore, we need a set of orthogonal functions in the z direction, Ba(z ). and there is a group of wavefunctions 'PJ for each a. For a symmetric collision. like those considered in this first paper, this is the obvious thing to do. However, in asymmetric collisions, like those considered later in Sections 7 and 8, this requires that the same set of orthogonal functions be used for the two different nuclei. We choose these wavefunctions, Ba(z), to be one-dimensional harmonic oscillator functions . For the 16<> and 40Ca nuclei considered in this first paper, oscillator states are a fairly accurate approximation to the ground state - 39- wavefunctions and for oscillator states the factorization, 1/J; (r)=r,o;(x,y)x;(z) is exact. For heavier nuclei we might consider other possibilities for the transverse wavefunctions but the 3D HF wavefunctions can't generally be factored as above . Fortunately, as shown in the paper, factoring and using oscillator states doesn't cause much change from the 3D results . The static HF calculations begin with 3D oscillator states, so we simply freeze these wavefunctions in one direction and use the imaginary time step method with the separable approximation to find the ground state wavefunctions SD;(x ,y) . We then choose the time independent oscillator frequency to minimize the energy of the static solution. The actual incorporation of this approximation into the computer code and the savings obtained are discussed in Appendix A5 for the full nonlocal Skyrme force (protons and neutrons treated separately). For this paper, we use the much simpler BKN force, isospin degenerate and without the (p-r-j 2 ) term, but otherwise the discussion is applicable. The major computational effort, both in time and storage, is the evolution of the wavefunctions. This evolution is reduced to a 2D problem resulting in considerable savings. However , the entire problem is not two dimensional and the calculation or the energy and potential will still require 3 dimensions. The results we obtain with this approximation for symmetric collisions of 18 0 or 40Ca and the comparisons to the 3D results are described in detail in the paper. The conclusion is that this approximation is not only better than the axially symmetric approximation, but is surprisingly accurate. Therefore, we will use this approximation for more TDHF calculations with heavier nuclei and more realistic nonlocal Skyrme forces . The insensitivity of the results to the choice of oscillator parameter is also important. For both 0 and Ca, runs done at the most sensitive angular - 40 - momenta, those at the fusion limits, showed that a 10% change in the oscillator parameter would lead to, at most, one unit (h) of angular momentum difference in the fusion cross section or in comparable trajectories. This is undoubtedly related to the success of the approximation, since even small motion in the transverse direction could be important if the results were sensitive to the shape in this direction. It is also fortunate since the transverse wavefunctions for heavier nuclei are more complicated than harmonic oscillator functions, and the factorized form limits the accuracy of the descriptions of these heavy nuclei in the z direction. Also, this insensitivity is almost a necessity for the calcula- tion of asymmetric systems since the requirement of one set of wavefunctions Bo.(z) in the transverse direction limits the accuracy of the shapes of the two different nuclei. -41- A SEPARABLE APPROXIMATION TO TIME-IEPENDENT HARTREE-FOCK CALCULATIONS S. Koonin, B. Flanders, H. Flocard, and M. Weiss Phys. Lett. 77B, 13 (1978) -42- A SEPARABLE APPROXIMATION TO TIME-DEPENDENT HARTREE-FOCK CALCULATIONS* S. KOON!Nt and B. FLANDERS W. K. Kellogg Radiation Laboratory California Institute of Technology, Pasadena, California 91125 USA H. FLOCARD Division de Physique Theorique Institut de Physique Nucleairc, B. P. N.<?. . . 1, 9ll+OG Orsay, France and M. WEISS Lawrence Livermore Laboratory University of California, Livermore, California 9!~550 USA Abstract We discuss an approximation to three-dimensional (3-D) TDHF calculations for heavy-ion reactions which neglects motion normal to the scattering plane. Quantitative comparisons with 3-D results indicate the accuracy of this method and an orderof-magnitude reduction in computation time. * Supported in part by the National Science Foundation [PHY7683685, PHY77-21602] and U.S. Energy Research end Development Contract No. W- 71~05-Eng-48.) tAlfred P. Sloan Foundation Fellow. -43Numerical applications of the time-dependent Hartree-Fock method (l] to heavy-ion reactions (2-5] and fission [6] provide a microscopic description of large-amplitude nuclear dynamics. Several calculations have approximated these systems in two non-trivial spatial dimensions by constraining the TDHF determinant to be axially symmetric (2 1 3 1 6]. However, recent fully three- dimensional (3-D) solutions for heavy-ion reactions have shown that tangential degrees of freedom are important for realizing fully the dissipation inherent in a mean-field approach [4,s]. There are also indications that tri-axiality is essential for an accurate TDHF description of fission dynamics [6]. Since "exact" 3-D solutions presently involve a great deal of computational effort for large systems, it is desirable to seek approximations which limit the number of variables to be treated explicitly, but which retain those degrees of freedom necessary for an accurate solution. This Letter explores one such approximation for heavy-ion reactions, in which motion normal to the scattering plane is forbidden (7). with 3-D calculations for 16 0 + 16 0 and 40 By direct comparison ca + 40 ca reactions we show that the resulting solutions are surprisingly accurate. The associated order- of-magnitude reduction in computation time therefore opens the possibility for meaningful calculations of larger systems using more sophisticated effective interactions. The motivation and specifics of a TDHF approach to nuclear dynamics have been discussed elsewhere [1,2,4]. The fundamental feature new to this work is that the complex time-dependent single-particle wavefunctions are constrained to be of the factorized form proposed by Devi and Strayer [7) (1) The coordinate normal to the reaction plane is z, while (x,y) are the -44coordinates in this plane. The real functions X are time independent, so that only the complex functions ~ involving two spatial variables must be evolved in time. The ansatz (1) is thus comparable to the axially symmetric factorization used previously [2,3 1 6], although here the coordinate z is "frozen" rather than the azimuthal angle about the symmetry axis. Note that non-axial configurations are still permitted by the form of (1). The motivation for the above factorization can be found in features of mean-field dynamics familiar from previous calculations [2,4]. In a TDHF approach, the evolution of the single-particle wavefunctions is governed only by their common, time-dependent mean field. Initially, the determinant is constructed to represent approaching, self-consistent static HF solutions for the colliding ions. During the collision, as orbitals propagate freely in the nearly constant mean-field of the nuclear interior (l], their motion is predominantly in the scattering plane. This can readily be seen in fig. 2 .. of ref. [2 ] , which represents a 4o Ca + 4oc a co 11 1s1on. There, despite the dramatic dynamics along the symmetry axis, the transverse extent of the system changes little during the collision. There is consequently little tendency to deform the system normal to the plane, and an approximte self consistency in this direction is maintained. The equations of motion for the functions ~ follow from the TDHF equa- For an effective interaction which results in a local HF potential w, tions. they are of the form J 02 ) . -it= o~. [- 'fl.2 ( 02 · 2m 2 +2 ~· ox oy + (Tz).J + W.(x,y) J J l.'fl. Here, W. is the projected HF potential J (2) -452 /x, y) =[co dz lx.(z)l W(x,y,z), J (3) W -oo and (T ). is the kinetic energy of XJ.: z J (4) Our calculations employ the energy functional of refs. [2,4,5]. This includes zero-range two- and three-body forces, a finite range Yukawa interaction, and the direct Coulomb energy. tion of the ground-state properties of local HF potential. It is known to give a fair reproduc16 0 and 40 ca [1,5] and results in a The functions X are chosen to be one-dimensional hannonic- oscillator functions and are assigned to the wavefunctions ~ at t = 0 according to the three-dimensional harmonic-oscillator classification of the static HF The time-independent oscillator frequency solutions. is chosen to minimize the energy of the individual static HF solutions before collision. Finite difference approximations to the equations (2) are derived by variation of a discretized "5-point" energy functional defined on a uniform 1 fm cartesian mesh, as described in refs. [2,4]. The Poisson and Helmholtz equations for the one-body Coulomb and Yukawa potentials are solved by the conjugate-gradient method discussed in ref. [4], and the functions~ are evolved in time using the exponential method of that reference. Time steps of 6t = 0.002-0. oo4 x 10-21 sec are used and energy and norm conservation for the entire system are equal to or better than the values quoted in ref. [4] for the 3-D calculations. One evidently necessary condition for the accuracy of the separable approximation (1) is that it reproduce adequately the exact static HF -46- solutions, at least in the reaction plane. The properties of these solutions for several oscillator closed-shell nuclei (N = Z = 8, 20, 40, 70) are shown in Table 1. Both the separable and 3-D solutions have been obtained with the imaginary time-step technique of ref. [4]. For each separable solution, the oscillator frequency for the functions X, ~w, has been chosen to maximize the binding energy per nucleon, BE/A. The binding energies of the separable solu- tions are consistently somewhat smaller than the 3-D values since they correspond to a more restricted class of variational wavefunction. The discrepancy increases with mass as oscillator states become a poorer approximation to the exact HF wavefunctions, although even for A = 140 it amounts to only 0.3 MeV per nucleon. 2 The root-mean-square radii, (r )i, of the two solutions agree to within 1% for all values of A. 2 The values of Q/A(r ), where Q is the mass quadrupole moment, show that even though the separable approximation treats one cartesian coordinate differently from the other two, it results in only slightly oblate solutions. Note that the deformations listed for the 3-D solutions indicate the accuracy of our finite-difference techniques, since these should be zero in the exact HF solution. In fig. 1 we show the coordinate-space densities of the separable and 3-D static solutions. In the reaction plane (z = o), the separable results agree well with the 3-D calculation, particularly in the surface region. In the direction normal to the reaction plane (x = y = O), where the separable approximation is expected to be worst, the agreement is poor. However, this appears to be of small significance for dynamical TDHF calculations, which depend strongly onl y on those properties of the ions within the reaction plane. In fig. 2 we compare separable and 3-D time-dependent trajectories for the 4o Ca+ 4o Ca system at Ecm = 139 MeV. used are defined in ref. . . The fragment separation and angular coordinates [4] in terms of the over-all mass-quadrupole tensor -47- of the system and the centers-of-mass of the two fragments. The results as a function of impact parameter follow a pattern familiar from previous 3-D calculations [4,5): highly inelastic events at small impact parameters, fusion at intermediate angular momenta, and inelastic scattering for peripheral collisions. The overall qualitative behavior of the separable approximation agrees excellently with the 3-D calculations. The relatively minor discrepancies in the L = 20 ~ (highly inelastic) and L = 30 ~ (fusion) trajectories indicate that the nuclei are somewhat too "stiff" in the separable approximation, and scission is therefore more difficult. This trend is confirmed in the L = 80 ~ trajectory (inelastic scattering), although the discrepancy here is actually smaller than it appears since the deflection function is changing rapidly with L in this region. and for 16 0 + 16 Trajectory plots for 40 Ca + 40 Ca at other bombarding energies 0 collisions show similar agreement between the separable and 3-D calculations. Final-state results for some are quantitatively compared in Table 2. 16 0 + 16 0 and 40 ca + 40 ca collisions Separable values for the final orbital angular momentum, final fragment kinetic energies, and center-of-mass scattering angle are generally in very good agreement with the 3-D calculations. particularly that in every case save one ( 40 Ca + 4o Ca, E. = 96 MeV, L. 1 1 Note = 10 ~) the fusion region is reproduced by the separable approximation, in contrast to the axially symmetric approximation where no fusion is found at all in these collisions [2]. Note also that the apparently large discrepancies in Ef for 4o Ca + 4o Ca at E. = 400 MeV amount to only some 20 MeV out of over 300 MeV 1 dissipated in the collisions. The largest differences between the separable and 3-D results generally occur where the final state is very sensitive to the initial conditions. L. 1 Thus, the separable solution for 40 Ca + 40 Ca at Ei = 96 MeV, = 10 t'1 crosses the outer barrier to scission with less than 2 MeV of radial -48- collective kinetic energy, while the 3-D calculation is just above the low-L cutoff for fusion. Similarly, the slight discrepancies at E. = 96 MeV, L. = i i 60 h and Ei = 192 MeV, Li = 80 h may be attributed to the rapid change of scattering angle with impact parameter in the region just above orbiting. Since such singular regions of L-space are only a very small part of the total cross section, the separable approximation accurately reproduces both the 3-D fusion cross section and the deflection function. We have also investigated the sensitivity of several of our calculations to choices of ~ other than that which minimizes the energy of the static solutions. The separable results pre- sented in Table 2 are relatively insensitive to variations of ± 200/c in thus parameter. · In summary, a quantitative stu d y o f 16 0 + 16 0 an d 40 c a + 40 c a co 11·isions · has demonstrated the accuracy of a separable approximation to the 3-D TDHF equations. The time-independence of the wavefunction in the coordinate normal to the reaction plane results in almost an order-of-magnitude decrease in cornputational time. The techniques presented here may be extended inunediately to Skyrme-like forces having an effective mass and allowing for isospin degrees of freedom. Large-scale calculations with these methods for mass synunetric systems should then provide a decisive test of mean-field theories for heavyion reactions. Furthermore, while there appears to be no need to do so at present, the separable determinant can be made somewhat more flexible by allowing for a simple time-dependence in the transverse coordinate; e.g., a time-dependent oscillator frequency. Self-consistent equations describing such motion can readily be derived from the time-dependent variational principle [8]. The method can be simply extended to moderately mass- asymmetric systems by using the same ~ for each nucleus; the relative insensitivity of the results in Table 2 to hw favors such an approach. However, the extension -49- of this technique to strongly asymmetric systems is not at all obvious and work in this direction is in progress. We are grateful for discussions with Drs. Paul Bonche and John Negele during these calculations. -50- References (1) P. Bonche, s. Koonin (2) S. E. Koonin, K. T. R. Davies, V. Maruhn-Rezwani, H. Feldmeier, and J. W. Negele, Phys. Rev. C13 (1976) 1226. S. J. Krieger and J. W. Negele, Phys. Rev. Cl5 (1977) 1359. [3) V. Maruhn-Rezwani, K. T. R. Davies and S. E. Koonin, Phys. Lett. 67B (1977) 134; A. Kerman, to be published in Proceedings of the International Conference on Nuclear Structure, Tokyo, 1977. [4) H. Flocard, s. E. Koonin and M. s. Weiss, Phys. Rev. c, in press. [5) P. Bonche, B. Grmmnaticos and S. E. Koonin, Phys. Rev. c, in press. (6] J. w. Negele, s. E. Koonin, P. M"oller, J. R. Nix and A. J. Sierk, Phys. Rev. c, in press. [7) K. R. Sandhya Devi and M. R. Strayer, Phys. Lett. 77B, 135 (1978). (8) A. K. Kennan and S. E. Koonin, Ann. Phys. (N.Y.) 100 (1976) 332. -51- Table l Properties of the static solutions 160 4oca A = 80 A= 140 BE/A (MeV) 2 1. ( r ) 2 ( fm) 2 Q/A(r ) ~w (MeV) 7.29 2.61 -0.014 13.9 ( 7 .43) (2.60)- (-0.006) 8.11 3.38 -0.003 (8.30) (3.35) (-0.003) 7.84 4.17 -0.022 (8.10) (4 .14) (-0.001) 6.91 5.00 -0.020 ( 7 .23) (4. 96) (-0.000) Values in parentheses correspond to the 3-D solutions. 10.9 9.1 7.6 Ei and Li are the initial c.m. energy and orbital angular momentum; Lf' Ef' and Gf are the final orbital angular momentum, c.m. fragment kinetic energy, and c.rn. scattering angle. values. 3-D results are given below the separable -53- Figure Captions Fig. 1. Comparison between separable and 3-D static HF solutions for several oscillator closed-shell nuclei. Curves labeled by z = 0 (x = y = 0) are the separable solutions in (nonnal to) the reaction plane. The fully 3-D solutions are the solid lines. Fig. 2. Comparison between separable and 3-D trajectories for 40c a + 40ca co 11'isions · at a center-o f -mass energy E i = 139 Me v • In cases where the 3-D trajectory is not shown, it is too close to the separable result to be plotted. -54- 0.2 A= 80 -'E -- 0.1 --.... ... I ") .. • 00 0.2 2 4 r (fm) ·., ·.'. ' ,_ 6 0 I ' ·...... - ....' ' 8 ~- '' .•. •• \ E -- ~0.1 A= 140 :.\ •• •• .... 00 6 \ .·.. ' • - 4 r (fm} \ rt> 2 2 4 r (fm) 3-D 6 0 2 -----SEP (z =0) Fig. 1 4 6 r (fm) ............. SEP {x =y = 0) 8 -55- oo 4oca + 4oca E10 b = 278 MeV -SEP -- -12 90°1----11--_....-+___.._.,__-+-~-f-~..f.----li----+~-+~-+-~-f'lr-~+--~f----190° 12fm rII rII rII <.000~ rrrr ;::JI =I =I O-NUJ 000 II II II 11 ooo '=I =1 =I 180° Fig. 2 t - 56- Section 7: The JiUling Approximation The previous calculations were for closed shell nuclei, 160 and 40 Ca. For open shell nuclei, however, the static HF solutions will be deformed. Therefore, the previous initial conditions are ambiguous because we must also specify the initial orientation of the ions . Difierent orientations have been shown to lead to substantial differences in the fusion cross sections calculated for 12 C+ 12C [Nilsson. Dhar, Koonin unpublished], so a plausible presumption is to average over all such orientations . However, heavier nuclei are relatively less deformed than was 12 C in these calculations ({:J= 2 (a-:) a+ ~0. 4 with a,b the major and minor axes) so that such orientation effects are probably small for energies not too near the Coulomb barrier . Averaging over orientation would substantially increase the computation involved, even beyond the need to perform several runs, at ditferent orientations, for each energy and angular momentum. Lack of reflection symmetry with respect to the reaction plane would double the size of our spatial box; and, in fact, a general orientation would no longer permit the separable approximation to be made since the values of the wavefunctions out of the plane will depend on both the values in the plane and the tilt of the nucleus . 'Ib.erefore, we adopt a filling approximation in which the valence nucleons are uniformly distributed throughout the valence shell, resulting in spherical HF solutions. This represents a crude orientation averaging. In the filling approximation, the density matrix is generalized as p(r,r',t)=~n{l/IJ(r.tfil'/(r',t) where the n; . (7.1) ; are time-independent occupation numbers determined by the ground states of the colliding ions . These n; are one for filled shells and a fractional value for the partially filled shells. Since these n; are not all zero or one, the many-body wavefunction cannot be expressed as a single Slater determinant - 57- and these calculations are no longer true HF or TDHF calculations. (In effect, this creates a wavepacket of difierent orientations and binds them together by the use or a common mean field.) However, once the densities have been defined, the TDHF evolution scheme is still applicable. Since the particles are independent, moving in a mean field and conserving their normalization, the evolution should be just like TDHF as long as the contribution of each particle to the mean field is properly normalized. Therefore, we regard the filling approximation as only a slight generalization of the ordinary theory. To actually incorporate this approximation into the computer code, it is much more convenient to normlize the 'V!i to ni rather than 1. <'l/.'i 11/li >=ni . . Therefore, the calculation of the densities Pv· -r,, and }q doesn't need to be changed. Instead, relatively minor changes are made in the code: in the normalization routine, the calculation of single particle energies, and the calculation of the mass and charge dispersion. The single particle energies become e;=~'l/lilh 11/1;> ni (7.2) and the total single particle energy becomes ESPT=L;niei i (7.3) The dispersion of proton (or neutron) number becomes (7.4) where the operator Di means that the integral in r is done over only one fragment (see Eq. (4.22)) and the sums are over only protons (or neutrons). Use of the tlll.ing approximation results in static solutions with smaller binding energy. For the case considered in the next published paper (BO 79), 28 Si, this difierence is quite large. Shown in Table 2 are energies for several systems calculated with the BKN force . - 58 - Table 2 nucleus exo. u1 0 -127.6 40 -342.l Ca 28Si prolate -236.5 28Si oblate 26 Si filliru!: sep. aoorox. -116.9 -324.9 -208.4 -205 .0 -185.2 Oxygen is a closed shell nucleus, in an oscillator basis (~, ~, 7'l.z) the occupied states are (000), (100), (010) and (001). Calcium is the next closed shell nucleus, including the states (110), (101), (011). (200), {020) and (002). The first two silicon solutions used n;=l and included only 3 of these valence states; prolate about the z axis used (101), (011). and (002), while oblate used (110), (200) , and (020). (The HF solution using the separable approximation will retain the oscillator states in the z direction and the symmetry of the wavefunctions in the x-y plane .) The third silicon solution used n;=* for all six states in the valence shell. Therefore, the binding energy discrepancy mentioned in the paper is due to the filling approximation rather than the force or the separable approximati on. The second paper included here describes in detail calculations for two systems 160+ 40 Ca and 28 Si+ 28 SL These would lead to the same compound nucleus 158 Ni but use of the filling approximation for 28 Si results in a di.fierent system. Although both systems have 26 protons and neutrons, the u10+40Ca calculation uses 14 wavefunctions whereas the 28Si+ 28 Si calculation uses 20 wavefunctions. Therefore, the silicon system has more degrees of freedom to dissipate energy. but it also has more total energy and more energy to dissipate at a given born- barding energy because the nuclei e.re less tightly bound. The energy differences seem to dominate, resulting in both a smaller maximum angular momentum for fusion and the opening of the low-l window for fusion at a smaller energy. This seems reasonable since some of the degrees of freedom in the silicon case are, in a sense, fractional degrees of freedom. Regardless of their - 59 - relative size, however, these two effects obscure any entrance channel effE' cts which could be discovered by a detailed comparison of these calculations. (Just as possible entrance channel effects throw uncertainty into the comparison of the two effects above .) Both of these calculations can be compared with experiment and both give reasonable agreement, as discussed in the paper. The 160+ 40 Ca calculations are compared with 16 0+ 40 Ca experiments. The fusion cross sections show quite good agreement but the experiments don't extend to high enough energies to give evidence of the low-l fusion ·window. The 28 Si+ 28 Si calculations are compared with 52S+ 2 ? Al experiments since there were no data for 28Si +26Si at the time of publication. This can be done because single-particle effects aren't calculated properly anyway and the bulk dynamics should be very similar for these two systems . The experiment doesn't measure fusion cross sections but does give evidence of a high-l and low-l cutoff for fusion in quantitive agreement with the calculated results . The calculation also gives quantitative agreement for the energies of scattered products just above and below these fusion cutoffs . The experimental references of this paper (24,25) have since been published (VI 79, NA 79) . The major addition is that 180+40 Ca seems to have too few deep-inelastic events for the low-l cutoff to be correct, however, no numbers are quoted . There are at least two experimental groups actively looking for low-l cutoffs in systems relevant to TDHF predictions . R. Vandenbosch, et al are study- ing 16 0+ 160, and J. Barrette ,et al are studying 28 Si+ 28 Si. The problem is difficult, however, especially because or the small cross section at these low l . -60- TDHF CALCULATIONS OF FUSION CROSS SECTIONS FOR 160 + 40ca AND 28si + 28si P. Bonche, K.T.R. Davies, B. Flanders, H. Flocard, B. Grammaticos, S.E. Koonin, S.J. Krieger, and M.S. Weiss Phys. Rev. C 20, 641 (1979) -61- TDHF CALCULATIONS OF FUSION CROSS 28 28 40 16 0 + ca AND SECTIONS FOR si + si P. Bonche a) . , K.T.R. Davies b) , B. Flanders c) , H. Flocard d) , B. Grarrunaticose), S.E. Kooninf), S.J. Kriegerg), and M.S. Weissh) a) Service de Physique Theorique, CEN Saclay, BP n°2, 91190 Gif-surYvette, France b) Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37830 c) W.K. Kellogg Radiation Laboratory, California Institute of Technology , Pasadena, CA 91125 d) Division de Physique Theorique, Institut de Physique Nucleaire, BP n°1, 91406, Orsay, France e) Service de Physique Theorique, CRN de Strasbourg, BP n°20, Cro, 67037 Strasbourg, France f) Alfred P. Sloan Foundation Fellow, W. K. Kellogg Radiation Laboratory, California Institute of Technology, Pasadena, CA 91125 g) Department of Physics, University of Illinois at Chicago Circle, Chicago, IL 60680 h) Lawrence Livermore Laboratory, University of California, Livermore, CA 94550 -~2- ~NUCLEAR REACTIONS 16 oc 40ca,x) and 28si( 28si,x) in time-dependent] ~artree-Fock approximation. Fusion and strongly damped collisions. ABSTRACT Time-dependent Hartree-Fock calculations are carried out for the systems of 56 16 0 + 40 28 . 28 . Ca and Si + Si. Cross sections for the formation Ni are qualitatively similar in both cases, although there exist significant quantitative differences which reflect the importance of the entrance channel. dow for fusion. Both systems exhibit an angular momentum win- The results of the calculations are compared with currently available experimental data. -63- I. INTRODUCTION Although the time-dependent Hartree-Fock (TDHF) approximation was formulatedl) 50 years ago, it is the last 3 years in which the approximation has been seriously applied to the calculation of physical quantities in heavy ion collisions. In this brief span, we have witnessed a remarkably rapid evolution from the early one-dimensional calculations 2 ), 3), which reproduced the qualitative features of colliding heavy ions, to increasingly more ambitious two-dimensional calculations 4 )-l 4), and finally 8 to three-dimensional calculations 12 )-l ) which yield quantitatively accurate descriptions of certain experimental data. In the present work, the TDHF approach is applied to the calculation of fusion cross sections . bo t h o f which . 1 ead to the f or two systems, 16 0 + 40Ca, an d 28Si. + 28Si, compoun d nuc 1eus 56 . Ni. Although the primary purpose of this work is to compare the results of the TDHF calculations with experimental data, we shall also discuss the accuracy of the two-dimensional rotating frame and separable approximations to the fully three-dimensional TDHF calculations. In addition, we shall discuss the effect of removing the isospin symmetry restriction on the nuclear wave function, and of varying the parameters of the nuclear potential. II. METHOD OF CALCULATION Most of the calculations presented here employ the parameterization of the Skyrme potential 19 ) referred to as Force I in ref. 18). In this version of the potential, the non-local terms in the mean field proportional to the Skyrme parameters t 1 and t 2 are identically zero. At the highest energies, the full Skyrme potential, including the non-local -64- terms.has also been used to calculate the fusion cross section. The ef- f ect on the fusion cross section of varying the 5kyrme parameters is also discussed in ref. 18), and in ref. 11). The calculations have been carried out using the two-dimensional . . 12),14) , and t h e two- d.imensiona . 1 rotating . f rame separabl e approximation approximation4),S),lJ) to the TDHF equations. In addition, in order to . . veri. f y t h e accuracy o f t h e separa bl e approximation, t he 16 + 40 O Ca fusion cross section has also been calculated at selected energies using the fully three-dimensional TDHF codelJ),l 7)• It should be noted that in effecting this comparison we have employed codes, both the separable code and the three-dimensional code, which enforce isospin symmetry. That is, neutrons and protons which occupy the same orbital are constrained to be degenerate. This degeneracy is not enforced in the rotating frame approximation4),S),lJ) code, and we shall therefore be able to conunent upon the effect of lifting the isospin degeneracy. As the various numerical methods which are employed to solve the TDHF equa~ions have been discussed extensively in the references cited, we shall not conunent upon them further here. · cross sec t•ions f or Th e f usion calculated as follows. 160 + 40c a and 28 5·i + 285·i have been For each of the energies studied, the maximum value of angular momentum for which fusion takes place, R> ~' and the minimum value of angular momentum for which fusion takes place, I<. ~. are separately de termined to a precision of approximately 2~. The fusion cross section is then calculated at each energy by the sharp cut-off formula 18) , - 1T O'fusion - k2 , (2.1) -65- in which k represents the relative wave number in the entrance channel. . system 16O + 40c a, t h e sum extend s over a 11 partia . 1 For t h e asymmetric 28 28 waves for which fusion occurs. For the symmetric system si + si, the sum in eq. 2.1 is replaced by twice the sum over all even partial The values of R for which fusion does waves for which fusion occurs. occur are determined at each energy by evolving the system at various impact parameters for a sufficiently long period of time to off er convincing evidence that the system will not undergo prompt fission. A more detailed discussion of the latter Point may be found in ref. 11). Although the sums to which we refer above may easily be evaluated exactly, it is consistent with the nature of the approximation to use, for either symmetric or asymmetric systems, the result cJf . ~ usion - k tv <~> +1)2 - (Q< +1)2 J [ = 1T~ 2 2µECM [(I.> +1)2 - (~, +1)2 J (2.2) in which µ represents the reduced mass of the system under consideration. III. RESULTS A. 160 + 40Ca We have performed calculations for the 16 0 + atory bombarding energies from 40 MeV to 350 MeV. 40 ca system for labor- The calculations have . . 1 separa bl e approximation . . 12),14),30) been carried out in the two-dimensiona , , . . . 4),5),13) • and in the two-dimensional rotating frame approximation I n addition, selected Points have been calculated using the fully three.. 1 13),17) • d imensiona TDHF code In all cases which have been compared, the separable fusion cross sections have been in complete agreement with the -66- fully three-dimensional TDHF calculations. A graphic demonstration of this agreement is given in fig. 1, in which the fragment separation coordinate r, defined as in ref. 5), is plotted as a function of time for a collision at a ·1aboratory energy of 290 MeV, in which the initial angular momentum,~, is equal to 48. The results of the two-dimensional separable approximation are in very good agreement with the results of the fully three-dimensional calculation, until t ~ 9xlotime the system has fused. 22 s, by which Plots of the quadrupcle moments of the mass distribution exhibit similar behavior. . 16 40 The angular momenta for which the O+ Ca system fuses are depicted as a function of the center of mass energy in fig. 2. The values plotted are derived from calculations which employ the separable twodimensional code with degenerate neutron and proton orbitals. The max- imum angular momentum which the fused system can sustain is approximately 60 units of ~. This value is in good agreement with the prediction 58 ~ 2 of the liquid drop model 0) for the A = 56 system. The fact that there exists a minimum angular momentum,P<(E), for which the 16 40 0 + ca sys- tem fuses, represents a dramatic prediction of the TDHF calculations 27 ). Referring again to fig. 2, we note that as the center of mass energy increases above 100 Mev,9< likewise increases. This reflects the in- ability of the colliding ions to convert sufficient translational energy to internal excitation energy. The inability of the transient fused system to dissipate more radial kinetic energy than a certain amount is an essential feature of the mean field dynamics 21 ). The increase of 2< with energy causes the fusion cross section to decrease more rapidly than l/ECM' which represents the energy dependence of Eq. (2.2) in the absence of}~, for energies greater than ECM ~ 160 MeV. At the latter -67- energy, the saturation angular momentum, k N 60, has been reached, and, in the absence of J<, the energy dependence is given simply by the kinematical factor l/EcM• The actual existence of a lower fusion limit could thus be inferred from the decrease with energy of the experimental fusion cross section for energies ECM ') 160 MeV. 16 40 . The 0 + Ca fusion cross section, as given by Eq. (2.2) is compared with the experimental cross section in fig. 3. The theoretical results presented have been calculated using both the results of the two-dimensional, separable approximation with constrained neutron and proton orbitals, and the two-dimensional rotating frame approximation in which no isospin symmetry has been enforced. As we have previously noted, the separable results are in complete agreement with the results of fully three-dimensional TDHF calculations (with the same enforced isospin symmetry). The cross section computed using the results of the rotating frame approximation exceeds the cross section calculated using the results of the separable approximation by approximately 120 mb in the energy range from just above the fusion barrier to ECM N.100 MeV. This difference may be attributed entirely to the enforced neutron-proton degeneracy in the separable calculations. At center of mass energies above 100 MeV, a significant amount of the dissipated energy excites nonaxial modes, and the rotating frame approximation breaks down. Thus, although the upper limit for fusion, J~N 50, is accurately given by the rotating frame approximation at ECM N 100 MeV, the approximation incorrectly indicates that the system will not fuse at J = 15. The latter results are qualitatively the same as those found in ref, 11), in which the lifting of the isospin restriction on the nuclear wave function 40 40 resulted in an increase of approximately 200 mb in the ca + ca fusion -68- cross section, and the excitation of non-axial modes caused the rotating frame approximation to break down in the 40 40 ca + ca system at a center of mass energy on the order of 100 MeV, the approximate energy at which the angular momentum window appears in the three-dimensional calculations. Referring to fig. 3, in which the fusion cross section is plotted as a function of the energy, we note that the theoretical results are in good agreement with the experimental results at low energy. This agree- ment at low energy affirms the facts that the nuclear sizes and the inter-ion Potential are accurately given by the TDHF wave functions and the nucleon-nucleon interaction. However at energies ECM ~ 100 MeV, the theoretical cross section exceeds the measured cross section by ~ 25 i. • In discussing the discrepancy between theory and experiment it should be noted that a significant amount of fusion followed by fission is consistent with the prediction of the liquid drop model mass 56 system. 22 ) for the Moreover such events would not have been observed in the evaPoration residue experiments with which we are making our com. parison. . . . However, separate investigations 23 ) d • d icate th at t h ere o not in 0 exists sufficient fusion-fission cross section to account for the discrepancy between the theoretical and experimental results. The fact that the TDHF cross section exceeds the experimental cross section especially warrants further study, since corrections to TDHF would be expected, at least naively, to result in an increased fusion cross section. We have also performed a single calculation at ECM = 208 MeV, using . . . . . ll) , an d the separable approximation with the f u11 Sk yrrne interaction the isospin symmetry restriction removed. The fusion cross section corn- puted on the basis of this calculation is 823±121 rob. Referring to fig.3, -69- we note that the cross section at this energy, as calculated using the results of the separable calculation with the local Skyrme interaction and degenerate neutron and proton orbitals is approximately 600 mb. The difference in cross section of the calculations is thus 220±120 mb. Since approximately 120 mb of this difference can be attributed to the removal of the isospin symmetry restriction, no definitive statement can be made as to the effect of the non-local terms in the Skyrme patential. -70B. 28 si + 28 si We have performed calculations for the 28 28 si + si system for lab- oratory bombarding energies from 70 MeV to 220 MeV. The calculations have been carried out in the two-dimensional separable approximation, with the initial, static wave function for the system determined in a filling approximation, in which the occupation of each ls-Od orbital of t h e separate d 28s· . . set equa 1 to 0 • 5 • i ions is . approx1ma. In essence, this tion represents an average over the possible initial orientations of the . o f t h e separate d d e f orme d Sl ater d eterminants 28 . . Si ions. As such it is a physically reasonable (and computationally necessary) approximation to the actual scattering. It nevertheless represents a departure from, and approximation to, a true TDHF calculation. The angular momenta for which the 28 28 si + si system fuses are depic- ted as a function of the energy in fig. 4. The maximum angular momentum which the fused system can sustain can be inferred to be approximately 16 40 so ~- Since the fused 0 + Ca system has been calculated to sustain 60 ~, it might appear that the maximum angular Jlk>mentum depends upon the entrance channel, and is not a function only of the total charge and mass of the system, as given by the liquid drop Jlk>del. This conclusion, how- ever, is considerably obscured by the fact that the tion has been effected in the filling approximation. 28 28 si + si calculaWe shall discuss this point Jlk>re fully at the conclusion of this section. As in the case 28 28 16 40 of the 0 + ca system, we find a fusion window in the si + si vs. 28 28 ECM plot. For center of mass energies greater than 54 MeV, the si + si system exhibits no fusion when the angular Jlk>mentum of the system is less than a cut-off value, ~JE), which increases as a function of bombarding energy. An unexpected result of the 28 28 si + si calculations is the -71- partial closing of the fusion window in the energy range above 100 MeV. This effect is manifest in fig. 4, in which the fusion region in the . p1otte d f or the 28Si. + 28Si. system. angu1 ar rnomen t um-energy p1ane is The results of the calculations effected in this energy range have been explicitly plotted. Thus, the crosses represent calculations in which • . . . t h e 28s.i ions scatter, and the dots represent ca 1culations in which the 28 si ions fuse. An island of non-fusion events is clearly visible in this energy range. While such dramatic variation, associated with only minor changes in the initial conditions did occur in the very first 2 one-dimensional calculations ), it is more probable that the appearance of the island signals a breakdown of the two-dimensional separable approxirnation. For a fully three-dimensional calculation effected in the cen- ter of the island of non-fusion events indeed leads to fusion. It is thus possible that the two-dimensional separable approximation breaks down at high energy, where the fusion cross section goes to zero. Accordingly, we shall assume that the island of non-fusion events is spurious, and shall ignore it in the computation of the fusion cross sections. . . f or 28Si• + 28Si. is . p1otte d as a f unction . The fusion cross section of energy in fig. 5. similar to that of 16 The behavior of the cross section is qualitatively 0 + 40 ca, although the latter system exhibits a. more pronounced plateau at intermediate energies, and drops more precipitously at high energy. the faster decrease of the A comparison of fig. 2 and fig. 4 shows that 16 0 + 40 ca cross section is due to the more rapid closing of the fusion window of the latter reaction. As the energy . 16 40 increases,9< increases more rapidly for the 0 + Ca system than for the 28 si + 28si system. It should be noted, however, that this high -72energy comparison may be complicated by the unexpected behavior of the 28 28 si + si system at energies above ECM = 100 MeV, as given by the separable approximation. In fig. 6, the center of mass energy loss and scattering angle, as calculated in the separable approximation, are plotted as functions of angular momentum for 28 28 si + si at ECM = 111 MeV •. The curves are qual- itatively similar to results obtained5 ) for 160 + 160 and 40ca + 40ca. However, an important distinction is the fact that the systems and energies investigated in ref. 5) did not lead to fusion. As discussed above, the scattering events observed in reactions in which the angular momentum of the system is in the range 27~ J. ~ 39 are quite possibly spurious. Fully three-dimensional calculations would likely l~ad to fusion in this interval. Although experimental data are not available for the system, measurements have recently been rnade 32s + 27Al system. 25 28 28 si + si ) on the simiiar In their investigation of the latter system at a laboratory energy of 175 MeV, Natowitz ex a1 25 ) can explain the observed spectrum by assuming the existence of two strongly damped components in the energy spectrum of the reaction products. One component is inferred to correspond to an angular momentum ~ ~ 48, not inconsistent with the liquid drop expectation that the A = 59 system will not support angular momentum ,.,, ~58. In contrast, the second component is inferred to corres- pond to very low values of angular momentum, ~ ~ 15. As a result of their analysis, Natowitz et a1 25 ) conclude that the observation of scissioning nuclei with such small angular momenta can be explained if the lowest partial waves in the entrance channel lead to strongly damped collisions, rather than to fusion. This situation is precisely that found in the -73285.l system. At ECM = 80 MeV, the same center of mass enerp,y as 32 27 . 28 in the 5 + A1 experiment, we find that the si + 28si fusion window is open for angular momenta 17 ~ j ~ 48. Further comparison can be made with the fragment energies for collisions at angular momenta just outside the window. Natowitz, et a1 25 ) find a center of mass fragment kinetic energy of 16.6 MeV for the R~ 15 component, and a fragment kinetic kinetic energy of 26.4 MeV for the~ ~ 48 component. of the TDHF calculations 28 28 si + 5i system at ECM = 80 MeV show that the fragments in collisions with ~ = 16 carry 16 MeV kinetic energy, and fragments in collisions with .l = 49 carry 28 MeV kinetic energy. The results of the 28Si. + 28s i. ca 1 cu 1ations . . . h t h.is interpretation . . are t h us consistent wit of the 32 5 + 27 Al data. We mentioned above that any conclusion as to the existence of entrance channel effects is significantly obscured by the fact that the filling approximation has been employed in the calculations for the 28si + 28si system. One such consideration is the fact that the Q value of the reaction is .considerably affected by the filling approximation. Thus, while the total binding energy of the separated 16 0 and 40ca ions, as computed using the local Skyrme interaction, is 442 HeV, 28 MeV less than the experimental value, the total binding energy of the separated 28si ions is but 370 MeV, 103 HeV less than the experimental value. The 75 HeV discrepancy in binding energy of the two systems can be used to explain, at least in a qualitative sense, the reason why the low angular momentum cut-off occurs at a much lower . t he 28s.i + 28s.i system t han in . t he 160 + 40c a system. energy in effects must be considered as well. i (16 + 40) Thus, for the 16 0 + Ot h er 40 Ca system, = 14 orbitals are involved in the calculation, while for the -741 (40 40) 28 . + 28 . " . . 51 Si system, b ecause o f t h e f"i 11 ing approximation, + = 20 4 orbitals are involved in the calculation. The intent of this discussion is to emphasize that no clear-cut statement as to the importance of entrance channel effects is possible on the basis of the present calculations. IV. SUMMARY We have calculated fusion cross sections for the systems and 28 si + 28 si, both of which lead to the compound system systems exhibit an angular momentum window for fusion. 16 40 0 + ca 56 Ni. Both 16 40 The 0 + ca 24 ) . reasona bl e agreement wit . h experiment . resu 1 ts are in , and the 28si + 28si results are consistent with the experimental data 25 ) on the . · 1 ar 32S + 27Al system. s1m1 . . 1 y s1m1 . . 1 ar f or the 160 + 40Ca The f us1on resu 1 ts are qua l"itat1ve and 28 ences si + 26) • 28 si systems, but do display s.ignificant quantitative differ- Thus, the 16 0 + 40 Ca system can fuse with an angular momentum l = 62, while the maximum · angular momentum with which the 28si + 28si system is calculated to fuse is e = 50. For comparison, the ·limiting angular momentum which the A = 56 system can sustain is given by the liquid drop model as R. = 58. The results calculated using the two-dimensional separable approximation are in very good agreement with the results of the fully threedimensional calculations, except in the case of th.e most energetic of . th e 2s5 i. + 2s5 i. reactions. The rotating frame approximation, as in past calculation~ is accurate just above the Coulomb barrier, but breaks down for energies greater than that at which a fusion window develops. The effect of removing the neutron-proton degeneracy increases the -7516 O+ 40 ca f usion · · b y approximate · l y 120 rob over t h e energy cross section range studied. The effect of the non-local terms in the Skyrme interac- tion is not clear at the single energy at which comparison was attempted. The difference in fusion cross section, as computed using the full Skyrme and local Skyrme interactions is not inconsistent with the difference which could be attributed to the fact that the full Skyrme calculation does not enforce isospin symmetry. Perhaps the most striking result of the TDHF calculations remains the prediction of an angular momentum window for fusion. While both the 16O + 4oc a and 285 i. + 285 i. resu 1 ts are consistent . . h experiment, . . wit it would be exciting if a definitive experiment could determine whether such a window indeed does exist. -76- ACKNOWLEDGEMENTS This research was supported in part by the National Science Foundation [PHY77-21602, PHY76-83685, NSF77-12879], and by the Division of Physical Research, U.S. Department of Energy, under contract W-7405-eng26 with the Union Carbide Corpbration. -77- REFERENCES 1. P.A.M. Dirac, Proc. Cam. Phil. Soc. 26, 376 (1930). 2. P. Bonche, S. Koonin, and J.W. Negele, Phys. Rev. Cl3, 1226 (1976). 3. P. Bonche, J. de Physique 4. s. Koonin, Phys. Lett. 5. S.E. Koonin, K.T.R. Davies, V. Maruhn-Rezwani, H. Feldmeier, S.J. 37 Collogue CS, 213 (1976). 61B, 227 (1976). Krieger, and J.W. Negele, Phys. Rev. ClS, 1359 (1977). 6. V. Maruhn-Rezwani, K.T.R. Davies, and S.E. Koonin, Phys. Lett. 67B, 134 (1977). 7. R.Y. Cusson and J. Maruhn, Phys. Lett. 62B, 134 (1976). 8. J.A. Maruhn and R.Y. Cusson, Nucl. Phys. A270, 437 (1976). 9. J.W. Negele, S.E. Koonin, P. Moller, J.R. Nix, and A.J. Sierk, Phys. Rev. Cl7, 1098 (1978). 10. K.T.R. Davies, V. Maruhn-Rezwani, S.E. Koonin, and J.W. Negele, Phys. Rev. Letts. ~, 632 (1978). 11. S.J. Krieger and K.T.R. Davies, Phys. Rev. Cl8, 2567 (1978). 12. K.R. Sandhya Devi and M.R. Strayer, J. Phys.~, L97 (1978); Phys. Letts. 77B, 135 (1978). 13. K.T.R. Davies, H.T. Feldmeier, H. Flocard, and M.S. Weiss, Phys. Rev. Cl8, 2631 (197S). 14. S.E. Koonin, B. Flanders, H. Flocard, and M.S. Weiss, Phys. Letts. 77B, 13 (1978). 15. R.Y. Cusson, R.K. Smith, and J. Maruhn, Phys. Rev. Lett.~' 1166 (1976). 16. R.Y. Cusson, .J.A •. Maruhn, and H.W. Meldner,. Phys. Rev.~, ?589 (197 8). 17. H. Flocard, S.E. Koonin, and M.S. Weiss, Phys. Rev. Cl7, 1682 (1978). 18. P. Bonche, B. Grammaticos, and S.E. Koonin, Phys. Rev • .£!Z_, 1700 (1978). ..;78- 19. T.H.R. Skyrme, Phil. Mag •• .!' 1043 (1956); Nucl. Phys. 2, 615 (1959). 20. S. Cohen, F. Plasil, and W.J. Swiatecki, Annals of Phys. ~, 557 (1974). 21. While the actual value of J." does depend upon the symmetries imposed upon the TDHF wave function, it is likely that the limit would remain in a totally unrestricted TDHF calculation. Although it is true that the removal of enforced symmetries, such as axial symmetry, and neutron-proton degeneracy, can enlarge the fusion window, o:up ·experience suggests that the limit will remain in a purely mean field calculation. 22. Referring to fig. 15 of ref. 20), we note that for values of 1~42, the fission barrier for A= 56 is calculated to be less than 8 MeV. Thus, on the basis of the liquid drop model, the A = 56 system is expected to fission for angular momenta 42 <I.< 58 • 23. D.G. Kovar, unpublished; F. Plasil, unpublished. 24. S.E. Vigdor~ D.G. Kovar, P. Sperr, J. Mahoney, A. Menchaca-Rocha, C. Olmer, and M.S. Zisman, unpublished. 25. J.B. Natowitz, G. Doukellis, B. Kolb, G. Rosner, and Th. Walcher, unpublished. 26. It should be iterated that the initial 28 si wave functions are not pure HF wave functions, but are calculated in the filling approximation. Thus, as discussed at the conclusion of section III, no clear cµt conclusion can be drawn as to the importance of entrance channel effects. -79- 27. The dynamical inhibition of fusion at low impact parameters has also been observed in Ca + Ca calculations based on the excitation of damped surf ace modes 28) , and in other macrophysics studies 29) • 28. R.A. Broglia, C.H. Dasso, G. Pollarolo, and A. Winther, Phys. Rev. Lett. 40, 707 (1978). 29. C.F. Tsang, Physica Scripta lOA, 90 (1974). K. Siwek-Wilczynska and J. Wilczynski, Nucl. Phys. A264, 115 (1976). 30. In the separable approximation, the required orthogonality of the orbitals (as a function of the coordinate perpendicular to the reaction plane) demands that the same harmonic oscillator parameter be used for both the 16 0 ion and the 40 . Ca ion. In the present cal- culation, the parameter employed is the mean of that appropriate to the individual ions. Since the latter parameters differ by but 107., the effect upon the calculation is negligible. This point will be fully discussed in the forthcoming publication by B. Flanders, P. Bonche, S.E. Koonin, and M.S. Weiss. -soFIGURE CAPTIONS Fig. 1 Comparison, as a fUnCt U)Il Of time I of the fragment separation . between coordinate in a co 11"ision 16 . 0 and 40c a ions at a labor- atory energy of 290 MeV with J. = 48. The solid line gives the result of the two-dimensional separable approximation, and the dashed line gives the result of the fully three-dimensional calcul ation. In both calculations, the neutron and proton orbitals are degenerate. Fig. 2 Fusion region in the angular momentum-energy plane for the system 40 16 0 + ca. The calculations have been effected with the separable code using the local Skyrme interaction, and with neutron and proton orbitals degenerate. Fig. 3 Comparison of the calculated fusion cross section with the re. sults of experiment 24) for 16 0 + 40 Ca. Plotted are the results of calculations using the separable code with the local 5kyrme interaction and degenerate neutron and proton orbitals, the results of calculations using the rotating frame approximation with the local 5kyrme interaction and unrestricted neutron and proton orbitals, and the result of a single calculation at 208 MeV using the separable code with the full Skyrme interaction and unrestricted neutron and proton orbitals. Fig. 4 Fusion region, in the angular momentum-energy plane, for the syst em 285.i + 285.i.. The calculations have been effected with the separable code using the local Skyrme interaction, and with neutron and proton orbitals degenerate. For energies above 100 MeV, the results of the calculations have been explicitly plotted, 28 5i ions scat28 ter, and dots used to indicate events in which the 5i ions fuse. with crosses used to indicate events in which the -81· F ig. 28 28 s ca 1 cu l ate d f usion · · f or · cross section Si. The calculations 5 l· + Fig. 6 have been effected using the code as described in fig. 4. 28 28 Deflection and energy loss functions for si + si at a center of mass energy of 111 MeV. The calculations have been effec- ted using the code as described in fig. 4. -82- 160 +oca ELAs=290 MeV 12 L = 48n - - separable 10 ----3-D 8 r(fm) 6 4 2 o~--~--~~~_._-----L-----L----' 0.2 0.4 0.6 0.8 t(10- 21 s) Fig. 1 1.0 1.2 1s 0 +4oca FUSION REGION 80 60 Lin 40 I 00 w I 20 0 40 80 120 160 EcM(MeV) Fig. 2 2 00 240 16 4 0+ °Ca FUSION CROSS SECTION EcM(MeV) 2000 250 150 100 70 50 40 35 30 I SEP.,LOCAL SKYRME,DEGEN. n-p I ROT.FM.,LOCAL SKYRME, ~ .c 1600 E g 1200 ..........,. ·-en :J r· ". . "". " ••••. .•. . . . ,, NON-DEGEN.n-p I EXP. ~. SEP.,FULL SKYRME,NON-DEGEN.n-p I:.f-1----I-------- ·~t~l····,__. .·.. --- .... 't- 0 800 I ~ I • ~:--..., ' 'I, 400 0 0.01 0.02 E~M(MeV- 1 ) Fig. 3 0.03 ', ~ 2sSi +2s Si FUSION REGION 80 60 Lili 40 I 20 0 ./'\,,.' ' ' ' ' ' 'A'·'- A"" I 20 40 ' I co V1 I ~~Gl~rt"" "~ 60 80 100 EcM(MeV) Fig. 4 120 28 28 Si+ Si FUSION CROSS SECTION EcM(MeV) 1600 250 150 100 70 I I ~1200 .c E 50 40 35 I '-""" I co c: °' I ·~ 800 :3 't- 0 400 0 I I 0.02 0.01 EcM1(Mev-1) 0.03 Fig, 5 -...~100 - --+----r I 8cM (DeQ.)-----180 140~',,, 100~ ' '',,,, 60~ • -~~F l 1b ~ I 20 LOS~/ 60 E CM (MeV) ,...__ 40 ...-- ~ ·~- I 20 ~I ~----.. I ,.- 70 I 60/I -60 I -100 -140 -180 180 ENERGY I ' LFUSION~ Fig. 6 I " I I I• I 80 I 00 ....., I - BB- Section 8: Results for 86Kr + 1s9La The final calculations presented here are separable, 2-D calculations for the heavier system, 86Kr + 1391.a, using the Skyrme III force. This system is heavy enough so that deep-inelastic processes dominate and the TDHF predications for the various equilibration processes can be obtained. For this larger system, the classical approximations discussed previously are more accurate and the bulk dynamics is less infl.uenced by single-particle behavior. This system has also been studied experimentally (VAN 78, DY 80) and with axially symmetric TDHF calculations (DAV 79,80), providing ample results for comparison with our predictions. The experimental results do show substantial deep-inelastic scattering but also show some fusion-fission events. For this system, the full, nonlocal modified Skyrme force must be used (again neglecting the spin-orbit force and occupying each orbital with two nucleons of opposite spin) as in Eq. (2.17) and (2.24). Since neither of these nuclei have closed shells in an oscillator basis, we also use the filling approximation. Thus, for eei<r the 36 protons fill the lowest three shells and the fractional occupation of the fourth shell is 0.8, while the 50 neutrons fill the lowest four shells and the fractional occupation of the fifth shell is 1 /3. For 189La, the 57 protons result in a fractional occupation for the fifth shell of 17 /30 while the 82 neutrons result in a fractional occupation for the sixth shell of 6/21. (However, tor real nuclei, including the spin-orbit force, both of these nuclei would have closed neutron shells N = 50, 82.) Static HF solutions for each nucleus are obtained by the imaginary time step method. For an oscillator parameter 1r.> =4.10, which maximizes the total binding energy of the system, the binding energies and rms radii of the two nucleons are given in Table 3. - 89 - Tc.b le ~ : nucleus i BE l (M1-::V'\ I rm:s., 1 rmsn i i a mL) I (frn) i I . 88lv. 6 '1~. 4 1391,a 1032.5 I 4 .86 I .;.2 ~ I 4.41 i I 5.01 I The experimental values are BEKr = 749.2 and BELa =1159.5. (All energies are in MeV.) From Table 1 of the first paper, we expect that the separable approxL'nation would cause about 20 MeV and 40 MeV less binding, respectively, than the full 3-D calculation. Because of this difference in binding energies, our calculations will have an extra 200 MeV of total energy at a given bombarding energy (collective kinetic energy) . However, to the extent that it is the initial relative velocity of the ions which determines the dynamics , this should be unimportant. The spherical aspect of the filling approximation should not present this kind of problem, ho·wever, since both nuclei have small deformations . ({3< . 1 see page 56.) This is related to the closure of the neutron shells . For the TDHF evolution of the system, a spacing ta= 1.2 fm and a box size 35 x 34 x 16 are used . The discretization errors can be estimated from the kinetic energy in the z direction (odd multiples of 1/2 1'.>) . These range from 0.3% for the first level to almost 6% for the sixth level. Indeed, discretization · errors in the boost result in a rapid initial transfer of about 16 MeV, or 53, fr om collective kinetic energy to internal excitation. This is worse than fer lighter systems because higher orbitals with more rapid spatial variation are included. These errors, however, are not related to the conservation of the discrete total energy because the operator h was derived from a variational principle. The conservation of energy (and norm) are dependent on the accuracy of the time evolution operator. Using ~t = .08 or.09 x 10-22 sec . and 2(4) terms in the exponential in the first (second) time step, total energy is conserved to 1 part in 600 and the normalization to 1 part in 5 ,000. A partial test run with M smaller - 90 by a factor of 2 showed less Lhan 1 part in 900 difference in the physical ob E.~rv­ ables : kinetic energy and the total and fragment rms radii. A sel of c.:dlculalions Wiis performed at E 1@ = 505 MeV lo determine:: bo ~ h the fusion region and the characteristics of the scattering , and two rur...s were performed at El® = 710 MeV to check the existence of a fusion region. Looking first at the fusion regions for E 1ab = 505 MeV, our results pr edict fusion for l =5 through l =60 (in units of h) while l =O and l =BO scatter. Thus, we would predict a fusion cross section 190±50 rnb . This agrees quite well with the experimental value of 170±50 mb (Lent = 66±10) which is identified as fusionfission yield from both charge and angular distributions . In contrast, the axially symmetric (AS) calculation shows no fusion at this energy. However , the Skyrme Il force was used in these calculations and a single run at l =O with Skyrme lII did fuse . The scattering at l =O differed from that in lighter systems because the lighter nuclei bounced back from the heavier nuclei rather than passing through it. This is undoubtedly due to the greater strength of the Coulomb potential. However, since the AS calculation should be exact for l =O, this effect is probably due to the slightly increased stiffness of nuclei in the separable approximation. This is demonstrated by comparing the rms radii of the systems. The AS result tor l =O has a minimum rms radius approximately 6.4 fm whereas our results have a minimum rms radius about 6.75 frn for any 1 ~ 80 . (This corresponds to a fragment separation about 9.8 fm.) At slightly larger l, the impact of the two nuclei is somewhat more gradual, giving the nuclear shapes more time to adjust. Runs performed at El.ab = 710 MeV show fusion for both l =O and l =50. This is consistent with both the experimental determination, 77±12, and with AS calculations using Skyrme II , which show fusion for !~65 . (For these two r :.ms, our minimum rms radius is about 6 .4 frn whereas AS calculations give about 6.2 fm.) - 91 - Figure 4 shows results for the final total kinetic energy and scattering angle in the center of mass. Our results are compared with both the axial calculations and the experimental contours for the double differential cross section. (The TDHF results are labelled by the initial angular momentum.) Both calculations reproduce the general characteristics of the experimental contours but with less energy dissipation, particularly for the fully damped events . The separable calculations provide slightly more dissipation in general than the axially symmetric calculations but the principle difference seemB to be a much larger small-angle contribution from the separable calculations. Experimental values are quoted for partially damped and fully damped cross sections - fully damped events are defined to be those below and to the left of the dashed line in Figure 4. These values are a(fully damped) = 1020±200 rnb and a (partially damped) = 770±165 mb . The axial results are quoted as a.ro~975 mb (l =0-100 and 125-175) and upDr:i:; 745 mb. Our results show l =0-90 and 110-180 in the run damping region yielding <1.FD~ 1100 mb. The partially damped region is hard to estimate in TDHF calculations because of the "energy damping" due to the numerical errors, so that even at l =250, when the nuclei pass without getting closer than 14.6 fm (their combined radii are only about 11.4 fm.), there is some energy loss. Although an energy loss of almost 20 MeV is not particularly good, it must be remembered that this is only 2% of the total energy of -1380 MeV. Both predictions of the fully damped cross section are consistent with the experiment. For small angle events (~cm <42°) the difference between the two TDHF calculations is much more dramatic. Our results indicate l =110-180 in this region giving a cross section ~ 800 mb . We also expect that fusion-fission events (with an angular distribution =~ = s~1'- ), will fall in this region approximately a quarter of the time, contributing another 40 nib. For the axially symmetric calculations, the comparable estimate is only - 92 - about 600 mb. Unfortunately, the comparable experimental number was not published. Figure 5 shows the total kinetic energy in the center of mass, Ecm, the scattering angle in the center of mass, 1'., and the interaction time, Tw as functions of the initial angular momentum for l~BO . The fusion region and the run at l=O, which should have fused, have been omitted. We see that Ecm damps below the Coulomb barrier very quickly. The experimental Coulomb barrier is about Ecm.=250 MeV, in agreement with the estimate (B.1) ro(Af +A/-) The final states have energy below the barrier because they scission with a very deformed shape (not two tangent spheres). Full energy dissipation is reached in about 2 x 10-21 sec. The scattering angle, 1'., is characteristic of a heavier system with a stronger Coulomb force. Thus, it is predominantly positive, like the AS results, although in contrast to these, there is a small region of negative angle scattering . For the highest l, 1'. is very close to pure Coulomb scattering. Our interac- tion times are generally larger than the AS calculation - sometimes by as large as a factor of two. Figure 6 shows the changes in charge and mass numbers, AZ and AA, the charge to mass ratio ZIA and the widths of the number distributions, r z and r A. Since the AS results for the width are quoted as full-width hal!-maximum, r. and the experimental results as dispersion, all, both labels are used for the axis, where r=..JBln.202. The AS results are quoted for AA and r.A while experimental results are quoted for AZ and I'z. - 93 - Our calculatio·1S show more mass transfer to the larger fragment ¥.ith a typical M~-3 for events with intermediate damping, compared to ilA~-1 for the AS calculations. Tht experime!''..al valuu are gn7en for two amounts d e!"". . 2 ~ ;;y damping : Q = -123 MeV gives ~z~-1 and Q = -:75 MeY gives ~z~+l. Therefo:· e, we seem to slightly overestimate mass and charge transfer, while AS calculations underestimate them. We never show mass transfer to the light fragment as do the most highly damped experimental events. but none of our calculations show as much damping either; our Qmax ~ - 110 Mev. Our value for AIZ shows an initial rise due to a greater transfer of neutrons to the heavy fragment and never seems to equilibrate to the ZIA of the total systern, denoted (ZI A)r. The AS results show the expected behavior ·with Z I A approaching and then fluctuating about (Z I A)T . Our results for the widths are somewhat smaller than the AS results, but they use a slightly different prescription for calculating the widths in the filling approximation because of ambiguities in the implementation. The AS calculations (DAV ) use a2=~'11.\Wii-~)tinj lwii i2 ' , (8.2) ij but since this isn't zero for the initial 86Kr nucleus. it is modified to a2(tru.e)=a2-uj;.cd C=O • (8.3) In either case, the results are lower than the experimental values. These values are u}=7 for Q = -125 MeV and uz=20 for Q = -175 Mev. Thus. for intermediate damping, we are about a factor of 4 too low, while for the full damping (which we do not get in our results). we are a factor of 10 too low. For comparison, the absolute limit on u2 for an independent-particle wavefunction ((DAS 79) and Section 4) is about 20 (depending on how the filling approximation changes the widths) close to the value observed experimentally. This underestimation of the widths is a common feature of TDHF calculations and is related to the fact that - 94 - u2 is a two-body operator. - 95- Section 9: Summary Jn conclusion, a one-body microscopic theory, TDHF using the Skyrme interaction, calculates quite successfully the gross properties of many different types or nuclear behavior for collisions or nuclei up to several MeV /nucleon above the Coulomb barrier. The separable approximation yields substantial agreement with both full 3D TDHF calculations and experimental results . Fusion cross sections are generally within about 10% of the experimental values for systems from 86.Kr+ 139 La . 16 0+ 16 0 to In heavier systems, the predictions of these separable approxima- tions for the deep-inelastic scattering are also quite accurate for energy vs. scattering angle and for the mean charges and masses of the fragments, except at the most fully damped events. The most fully damped events in these calculations show too little energy dissipation. The widths of the distributions, however, are not accurately given. The calculations predict widths smaller by a fac tor of 2 to 4 (for the most fully damped events). This isn't surprising since the widths are expectation values of two-body operators. The inadequacies of these calculations are generally assumed to be due to the TDHF approximation itself. However, in the absence of full 3-D calculations including the spin-orbit force and averaging over orientations, we cannot, with certainty, know how much or the inadequacy is indeed due to the TDHF approximation itseU. These calculations will, however, have to wait for more powerful computers or more powerful computer budgets . Indeed, they may be inappropriate because of the use of an effective interaction and the inherent limitations of the TDHF description. -96- PART II A STABILITY CRITERION FOR TDHF CALCULATIONS - 97 - Section 1: Introduction Large-scale numerical calculations have shown that the TDHF approximation can successtully describe the gross features of heavy-ion collisions. The situation remains unsatisfying, however, for two related reasons: there are as yet no quantitative criteria which might be applied to signal the brecik-down of the mean-field approximation in any given physical situation, and tractable, systematic corrections to TDHF which become important when such a break-down occurs have yet to be developed. Although we can otrer no guidance in solving the second problem, this contribution is relevant to the first problem, namely the development of a necessary criterion for the validity of TDHF solutions. Our approach rests on the following general considerations . It is well known that the TDHF equations for a system of A nucleons may be derived from a variational principle (KE 76), in which the action is made stationary with respect to the functional variation or the A complex single-particle wavefunctions parametrizing the TDHF determinant. However, the variational principle makes no statement about the second variation of the action, which is connected with the stability of the TDHF solution under small perturbations. Consider, for example, the situations shown in Fig . 7. Let 4>(t) be a. solution to the TDHF equation (one which makes the action stationary). and suppose that at some time t=O, 4>(0) is perturbed to a new determinant +(O) which differs only slightly from t(O). The subsequent evolution of+ can then lead to one of two situations. In Fig. 7a, +(t) remains "near" t(t) for all t>O and so 4>(t) is stable. Fig . 7b, i.P(t) is unstable, and +(t) diverges from t(t ). Of course, the stability or instability of CJ> may change with time, and might depend upon the type of perturbation applied to generate it (conditional stability) . The standard interpretation of the TDHF wa.vefunction is that expectation values of one- or few-body operators generated from it correspond to the aver- - 98 - age evolution of the exact Schroedinger wavepacket. Consistency in this interpretation requires that this mean evolution (and hence the TDHF path) be well defined; i.e ., that~ lie in a "valley" as in Fig. 7a, rather than on a "ridge", as in Fig . 7b. We, therefore, propose that the stability of TDHF solution with respect to small perturbations is a necessary {although not sufficient) criterion for the validity of the TDHF approximation. The RPA theory (TH 61) has been used to test the stability of static HF solutions; we propose a time-dependent generalization to test the stability of TDHF solutions. Section 2: RPA Theory Ordinary RPA theory attempts to calculate the energies of collective states which are similar to the rlartree-r ock ground stale . The theory looks ior periodic solutions for state:; "close" to the HF state. Because of the condition that these states be "close" to the HF state, about which we will be specific later, we are actually calculating the curvature of the energy surface for deterrninantal wavefunctions at the HF state . This curvature is the important property to test the local stability of the HF state. If all the normal modes have periodic solutions, the HF state mu.st be a local minimum and it will be stable . However, if any mode does not have a periodic solution, because its frequency is ima- ginary, then the HF state will be unstable with respect to that mode . As in the derivation of TDHF, we minimize the action frit,<'o/lili:t -Hl'o/> (2. 1) with respect to variations in 'o/• in order to derive equations or motion for 'o/. For the state '11 we want to use a general Slater determinant. Thouiess (TH 61) showed that any N-particle Slater determinant not orthogonal to the HartreeFock determinant ii> can be written in the representation appropriate for .P as (2.2) 1J is a normalization constant and a+ and a are creation and annihilation opera- tors, respectively. The subscripts p or h refer to states which are not occupied or occupied, respectively. in the original HF state 14>>= IT a,t IO> . (2.3) h=l To show that it is a determinant, we need to use the fact that all terms in which the same creation or annihilation operator occurs more than once vanish. (2 .4) - 100 Since there are only N different occupied states h we can stop the sum over n at N and rewrite the sum as a product of N terms. lir>=[fl(l + ~c.Ph0p+a,J][IJan-t: IO>] p h. (2.5) h.' Each term in h will give 0 unless h'=h. Therefore, lir>= Il(a,,,++ ~Cph<Jp+) I 0> h. (2.6) p Thus, we can see that ir is a determinant of the wavefunctions Xh =cI>h +~cph «I>p • (2. 7) p Since unitary transformations among occupied states will not change the Slater determinant of the wavefunctions, we see that the state iris a general determinant not orthogonal to cI> . States orthogonal to cI> can be constructed by using =p one or more wave functions of the type Z11. ~ ~ «I>p. We now see that we are calculating the TDHF evolution of the state '11 where its time dependence is contained completely inc.Ph, the particle-bole occupation numbers. We write the action in terms of cl'h , c.Ph, and c~ and accomplish the variation with respect to +• by a variation of c~ . so the evolution of '11 is described by the equation for c':p1i. as a function of cp1i. and c~ . To study the local stability of the HF state, we assume cp11. is small and truncate the calculations of +and the action to second order in cp1i.. Therefore. after the variation. the equation for cl'h will be linear in cJlh and c;i.. We look for periodic solutions by assurn1ng (2.8) and then solve for the normal modes and their frequencies <.J. These must all be real for the solution~ in Eq. (2.B) to be periodic. Of course, a larger rJ will imply a sharper local minimum and a more sharply confined HF state. - 101 Section 3: Time-Dependent Generalization of the RPA Theory For our time-dependent generalization of the RPA theory, we again minim- ize the action for a general determinantal wavefunction, (3.1 ) but in this case, the state 4> and the creation and annihilation operators for perticles and holes also evolve in time . We define an operator ~=a1,+ah. which creates a particle-hole pair and write (3 .2) Since we are still concerned with local stability, we still truncate the expressions for + and the action to second order in cph: (3.3) We also evaluate TJ to second order, by setting <'11J'11>=1. (See Appendix B1 for details.) (3.4) 80 that (3.5) For the action we need <+IHI+> and <+li:t I+>. ( h=l) The first is straightforward. <+I H I+>=( 1-e •·c) <H> +<He· A.+>+ <c \AH> (3.6) 2 +<c··Aflc·A+>+*<(c··A) H>+*<H(c·A+) > 2 We assume a second-quantized Hamiltonian -"'t ,,, 1 "' H - /..J 4ao.+a~+ 4 '-' T',.S~a.a.+a.~+a.4ri.y • o.'16 (3.7) where laJ is a complete single-particle basis, t 4 =<alt I"> (3.8) are the matrix elements of the one-body kinetic energy, and (3.9) - 102 - are the antisymmetrized matrix elements of the two-body interaction. The second expression, <ir Ii :t I ir>, is more complicated because the gen- eralization must agree with the RPA theory for ii> an HF state . If we evolved each wavefunction rp11, with the standard HF operator h as irP11, =h rp11,, we would get rp11,=e 1ht 'ht.!IF· where hrph.HF=e11,rphJIF. This would result in a time dependent Iil>>=exp(-il:eht) Iil>eF> (3.10) h rather than the time independent Iil>HF> needed to agree with ordinary RPA theory. Since this time-dependent phase is unobservable, the correct operator ii must include only the physically relevant changes in <I> . The usual h in the single-particle wavefunction evolution i¢u.=ho.prtJp. with hu.fJ=tap+ l: VimfJh. (since h pap=Oo.pOo.h see Appendix Al) is a unitary transformation mixing all the various states. However, since unitary transformations among hole states or among particle states will not produce any observable changes, we want to take the bole-hole and particle-particle components out of h. To do this , we can use the one-body density matrix p=~ l<Ph><rph I. which is a projector for Iii>>. If we h remove the hole-hole components, php and the particle-particle components, (1-p)h(l-p) we are left with h =ph (1-p)+ (1-p )hp (3.11) Ji,,;,·= O=hn'" (3.12) 80 that ~ =~ and ~ =h,,p Since we now include only the physically relevant changes i > is orthogonal to its time derivative . <tJ :t l<P>=-i<il>lph(1-p)+(1-p)hpl~>=O (3.13) With this evolution operator, the static <Pw is a solution to the evolution equation - 103 - .8 '1.-1 il>HF>=h Iil'HF>=O 8t (3.14) We can now evalute the other term from the action, <it Ii :t lit>=i[ -¥lfti.c··c)+c··b' +<i! .J+>+<c··.A c .J+> (3.15) +}2( <i! .J+c ..J+>+<c ..J+c ..J+> +<il> 11-c··c +c··l+c .1+ +c··1 c .1+ +}2(c··1) 2 +}2(c ·.A+) 2 1~> J We now know the correct evolution operator to use, so i I4>=h l ii» (3 .16) ia;t=[ii.~+] ; and ~=[h,a1J ~=[Ji ..A+] Therefore, the net result is <itli :t lit>=i[-}2fti.c··c)+c··e]+(1-c··c)<h>+<ii c ·A+> (3 .17) +<c··.A ii>+<c· ·.A ii c·.A+>+}2<li(c·A+) 2 >+}2<(c··.A) 2ii:» We can notice that after the first two terms, the right-hand side looks precisely like <it l HI ir> but with H replaced by h. Now we want to minimize the action by varying c~ but remembering that the total derivatives, such as :t (c··c) do not vary since the endpoints are .fixed. (3.18) and .A-<+li !t l+>=iC.Ph-c:Jlh <ii>+~ii>+~ii c·A+> 6c:Jlh u (3.19) +*(~c··.A ii>+<c··.A ~ii>) Since the stability of the TDHF path is determined by the curvature of the energy surface normal to the path. another modification of the ordinary RPA theory is needed. We require that the initial perturbation of il> be orthogonal to - 104- the TDHF path. (An initial perturbation along the path would be equivalent to a shift in time.) The direction of the path is given by :t 4J. so that !Ji is also orthog- onal to the path and the perturbed wavefunction "-¥ is orthogonal to :t !Ji . (3.20) Only the c·A+ terms contribute (See Appendix Bl for details.) The orthogonality condition on "-¥ becomes (3.21) These conditions are most simply imposed by the method of LaGrange multipliers, adding the terms (3 .22) to the action. The final result for the variation of the action with this constraint is ii:Jlh =-cl1h ( <H>-<.h >) + <Ap.i.. (H -h)>+<Ap.i.. (H -h)c ·A+> (3.23) +*[~c··A(H-h)>+<c··A ~(H-h)>]-~ The evaluation of these expectation values in Appendix Bl gives the result (3.24) where the matrix Ap,."J''h'=hpp•6M·-6ppJtn.•11. + Vph'hp' (3.25) is hermitian and the matrix (3.26) is symmetric. We can see that for a static HF solution with '7n=O, hpp·=e,,6pp·· and ht..·1r. =s1r.6M' where the e's are the single particle energies; this equation reduces to the usual RPA equation. To test the stability (or measure the curvature of the energy surface) at a particular point along the TDHF path of~. we must imagine - 105 - suspending the evolution of q, and studying the motion of '-Ir constrained to be perpendicular to this frozen direction of motion for cl> . Thus A , B, and h depend on the TDHF wavefunction cl> and vary from point to point in parameter space. However, while testing the hypothetical motion of '-Ir in a particular direction, at a particular point these are constant. Note that we are only able to equate choosing a particular point in parameter space with freezing the time evolution of~. A, B, and~ because we have chosen the correct evolution operator h to isolate only physically relevant changes. This means that Eq. (3.21) becomes (3 .27) for all times t in the motion of '-Ir and any particular time we choose to freeze the motion of If>, tF. Therefore, (3 .28) al.so and we can solve for A. by multiplying Eq. (3 .23) by~ and summing over ph . This gives (3 .29) +~(tF)Bph;p'h'(tF)c;.h' (t )] and substituting>.. back into Eq. (3.23) gives (3 .30) where ~ ~h;2he pJllt.;p2h2=6.PP26'M2- E11.p.,,..~.h' (3.31) 'h' We can put this equation into a more elegant form by noticing that P is hermiti.an and that it is a projector off the TDHF path, i.e. (3.32) We can write the equation as - 106 - ii:,m= ~ [~;p'h'cp'h'+B,m;p 'h'c;.h . ] (3 .33) p'h' where A=PAP is hermitian and B=PBP• is symmetric, just like the RPA matrices A and B. We look for periodic solutions by setting C,m (t) =X,m e -ir.it + Y;., e Hr.it (3 .34) and assuming r..> is real. Equating terms with e -ir.it gives CJXpn =~ [~;p'h'x;,'h'+B,m;p·h·Yp~·] (3 .35) -GJY.Ph = ~ [~;p'h' ~·h.+ B.Ph'3''h'x;'h' ] (3.36) p'h' and equating terms with e +ic.it gives p'h' In the usual RPA matrix notation. these equations are written. "'rn=[_~. -~·l m (3 .37) If all the eigenvalues r..> are real, then the TDHF path for ~ lies in a valley and the path is stable. Once again larger Ir..>j imply a more sharply confined path. Although a specific evolution operator h was chosen to facilitate the derivation of these equations, we would suspect that unitary transformations among particles or among holes would not change the eigenvalues, e.>. In Appendix B2 we show that these transformations would result merely in a similarity transfermation of the matrix I~. ~.1 -B -A and thus, we prove that the eigenvalues are not changed. This is important because the time evolution can still be defined merely by the particle-hole and hole-particle elements, h,,n and Ii,,,,, as in TDHF. Therefore, we will not need to evolve any additional wavefunctions or parameters to do the stability analysis . - 107 - Section 4: TDHF for the SU(3) Model To test our ideas about stability, we apply them to an exactly soluble model. We consider the SU(3) generalization (LI 70, HOL 74) of the familiar SU( 2) (Lipkin) model (LIP 65) . Although the TDHF approximation for SU(2) has be en invcstigated in some detail (KRI 77), it is not sufficiently complex for our purposes. There is only one collective coordinate and associated velocity (the rotation angles ~. If') and hence, no collective degrees or freedom exist normal to the TDHF path. The SU(3) model consists of N distinguishable identical particles, labeled by the index n. each of which can be in three single-particle states having energy ~i . The Hamiltonian is (4. 1) where Vti is the strength of the two-body interaction between states k and l . V~ = VUc and Y.t.t: = 0 (4.2) with (4.3) The operators a+ and a are creation and annihilation operators so that each term of ~ takes a particle, n, out of state l and into state k and each ~ counts the particles in state k . We want to write a general collective determinant, i.e .. in which all particles behave identically. We define new creation, annihilation operators a+, a by (4.4) where A is a unitary matrix independent of n . Now a.,! 1 creates the occupied TDHF orbital; a.,!2 and a.,!3 create the remaining two particle orbitals . Thus, a general collective determinant is written (4.5) - 108 - and the matrix elements AtJ determine the actual occupied and unoccupied states. Since these A ' s may be complex, six parameters would be required to completely define the TDHF orbital; however, since both wavefunctions (created by a+ or ex+) are normalized and the overall phase can be ignored, we are left with four parameters needed to define the TDHF state . Since the operators a+ are in a fixed basis, the time dependence of the operators a+ and ultimately of the TDHF state is contained completely in the matrix elements of A . Therefore, the TDHF equation governing the evolution of 4> must be transformed into equations governing the evolution of the matrix elements of A and the parameters used to characterize A . The TDHF equations are .B h "'wn= rfJn . (4.6) h=~"*'~G\'.~G\'.ni (4.7) where h is a one-body operator l:;l n The matrix elements ~ are independent of n because the particles are identical. We write the TDHF equation in terms of the a+•s as (4.8) and hrpn =~~~cxri·ka'.n 'lcx,!1 IO>=~h.t1~ IO> l:;l n' (4.9) J: Therefore, the equation for the operators a+ becomes an + 'l. . +1 =""'".LrU CXJ1A: • (4.10) " We know that the time dependence of ex is contained completely in the matrix A . (4.11) Th.us, if we define a matrix (4.12) we can write Eq. (4.10) as - 109 - (4.1 3) iT it =./1t; 1 We can easily show that iT i:; herrnltia...""l., just ash is hermitian. If we igr1ore the winding phase, given by tht · pw·ely imaginary T, 1, '""e tgs.:n need just 4 parameters for the two complex equations with k = 2 or 3. For our time-dependent generalization of the RPA theory, we have shown that we can ignore hole-hole and particle-particle components of the evolution operator, T 11 , T22. TP/:3 and T33 . Therefore, throughout the calculation we need only the four parameters used for a normalized general TDHF orbital (ignoring the overall phase). Using the matrix elements of h , which are calculated in Appendix B3, the TDHF equations can be written as two equations for matrix elements of A : (4. 14) for j =2,3. We choose, for our model Hamiltonian, equal level spacing (e 1 =-e, e 2=0 , and e3 =e) and equal interaction strengths (Vkl= ~ for k~l). We choose to parameterize the matrix A in terms of four rotation angles 1'1 ."'2,1Jl 1 , and 1/12 and write A as (4.15) C0~1 A= co~sim~ 1 e "'"°'11 cos~sin1' 1 e •ti -simj2sin1' 1 e '"'2 1+cos 2~(co~ 1 -1) sin"'2cos1'2 (co~ 1 -1)e -iC'fii--12) sin~s~ 1e -i1i! sin~cos"'2( cos1' 1 -1)e •C?i-"'2> 1 +~in21'2( cos1' 1 -1) Note that this A is defined to be the inverse of the transformation matrix in (KRI 77) . So for ,, 2 =0 this choice results in the same definitions : my 1' 1 is his a.., my 1/1 1 is his 1/J. The physical meaning of the "coordinates" ~ 1 and ~2 is quite clear. Since we can easily confirm that A is unitary in Eq. (4.15), we can write the occupied TDHF orbital (4 .16) -110- Therefore, the probability of finding n in the first level is cos 2 ~2 • in the second level is cos 2~sin2~ 1 and in the third level is sin~ 2 sin2~ 1 . Since all particles behave the same in the collective state, multiplying the above probabilities by the number of particles, N, gives the mean number of particles in each level for the TDHF state ii> (in agreement with<~> as calculated in Appendix B3) . The meaning of the "associated velocities" 'f/1 1 and 'f/12 is less clear but we can see from the evolution equations following, that 1/1 1=1/1 2 =0 will imply ~ 1 =~2 =0. In Appendix B4 we transform the TDHF equations for the matrix elements of A into four coupled nonlinear first-order differential equations for the four angles . (4.17) ,,P 1=-t(l +x[ cos2'f/11(-cos~ 1 +cos~~in2~ 1 )+cos2'f/12 sin2~2 sin2~ 1 -cos2('¥'1-lt2) sin2~sin~ 1 ]) ,,P2 =-t(2+x[cos2'f' 1 cos~sin~ 1 +cos2l'2 (-cos 2~ 1 +sin~sin2~ 1 ) ~os2(l' 1 -"2) cos~sin~i]) Note that for ~=O these agree with (KRI 77) 's TDHF equations. Here xis the usual strength parameter x=- V(N-1) . (4.18) t The determinant characterized by these four angles has a total energy H(~1.""2.l'1.'Yt2)=N e(sin2""2sin~ 1 -cos2~1-xsin2~1 (4.19) [ cos~ 1 ( cos ~cos2'f/1 1 +sin2""2cos2'f/12 ) 2 +co~sin2~2 sin2~ 1 cos2(1/11-1'2)]) (see Appendix B4). The TDHF evolution of ii> must conserve energy and substitu8H lion of the TDHF equations for the four angles into 8t shows that indeed, - 111 - a; =O. We can also notice that for x positive, an attractive potential between particles in different states, the minimum energy for a given '?9 1, '?9 2 occurs at 1/;1 =1/;2 =0. Thus, H('t9 1 ,'t92 ,0,0) can be viewed as the potential energy and we can orient ourselves by following the position of the TDHF state, '?9 1 and '?92 on a contour plot of the potential energy surface . - 112 Section 5: TDRPA for the SU(3) Model To apply our ideas for a time-dependent generalization of the RPA theory, we must cast our general equations (3.30-3.33) into the specific form suitable for this SU(3) model. Thus, a normally occupied state, previously denoted by h , will now be denoted by a particle n in state 1. A normally unoccupied state, previously p, will now be a particle n in state a=2,3. Thus, our particle-hole creation operator ~=a7,+~ . now becomes Knai=a.taan 1 and its adjoint~ now becomes ~la · Our particle-hole occupation numbers, cp11,. become Cnal• or dropping the irrelevant 1, Cna· The previous sums over p and h now become sums over n and a. Thus, our perturbed state becomes Ii'>=1]ffexp(})naatiatx.n i) l 4>> (5.1) fta with the TDHF state (5.2) Of course, for jir> to also be a collective state, all particles must behave identi- cally and Cna must be independent of n, i.e. Cna= -Jica where C11 =~cna· This res- " triction on the perturbation of Iir> will be used only after the equations have been solved. 'We now need only a prescription for 11.ph in order to write down the TDRPA equations for the SU(3) model. Before h,,A=~=<A,,hh> . Now this equation becomes (5.3) Since all particles, n, are identical in our full Hamiltonian H, these hnai must also be independent of n and are written hai· The term <Knia(H-h)> vanishes again just like ~ (H -h)> and the bole-hole and particle-particle components of h must vanish. With these substitutions, the TDRPA equations for cna become - 113 - Wna= l; Cn·a.C-Onn•Oua'( <H>-<h >)+<Kn111(H-h.)Kn,·11·1) (5.4) n'a' After the evaluation of these expectation values, as discussed in Appendix B3, we get (5.5) where .Anu;n•a·=Onn·(hua'-Oaa·hu)+ Ya1111·(1-0nn·) (5.6) Bn.a~·u= Vaa'll ( 1-0nn•) and . =" """' V*Pnq ~2Vi;~A,1~1\q,, (5.7) ij hw =l;A..t Atr + l; 2 Vii (N -1)Aj1 At1 A;A: ~ =h~ ( 1J Because the n dependence of the matrices A and B is of the simple forms 1 or Onn· we can rewrite the above equation, using Ca as Wna= l;(h=·-Oaa'h u-Ycrua' )cnu+ Ycr1 la' Ca' (5.8) a' +2; (- Vcr,,.u)c:a' +2; Ycra·11 c; -Mai . ,,. rt With the equation in this form, we see that all the coupling between different n's occurs in the collective terms Cv and c;. The constraint l;~(tF)c,ph(t)=O now becomes ph Lhia(t.P) Ccr(t )=O (5.9) CJ Just as in the general case, we solve for >.by (5.10) CJ We find the expression for C,, by adding the equations for each n. Thus, we get (5.11) With -114- . Aov=hoo·-~orth 11 + "' Vol la' (N -1 ) =Aa·a (5.12) and Baa·= Yaa'll (N-l)=Ba·u The matrices A and B have the usual RPA symmetries, A=A+ and B=BTr. Using the constraint from Eq. (5.10) to find >.in Eq. (5 .11) gives (5.13) 0 Now to solve the TDRPA equations, we .first solve the collective equation for ta. (5.14) where h 2 =I;h 1uh 01 . Because there are only two values for CT, the constraint f1 l;h 1uCu can be used to solve for f1 (5.15) This fortunate circumstance occurs because the one-body HF Hamiltonian couples the hole state to one linear combination of our two particle states. Excitation of the other linear combination is the single RPA mode transverse to the TDHF path whose stability we are testing. Thus, the equation for u=2 is (5.16) where A= h31h13 h21h13 h1i!l-s1 h21h12 h2 ~- h2 A32- h2 ~s+ h2 Ass (5.17) h31h1s h21h1s h21h1s hl1hrn h 2 B22- h 2 Bs2- h 2 B2li3+ h:uh 2 Bss (5.18) (which is real) and B= (which is complex) . We now look for periodic solutions (5.19) With c.> real. When we separately equate terms With e """"'' and terms with e +ic.it , - 115- we get ~~]=[-~. -~·]~:] (5 .20) Therefore, our collective mode has a frequency(,,) given by setting IA-r.> B I j-s• -A ·-r.> =O . (5 .21) Notice that a real frequency r.> is only possible for A real, as it is in this case, so (5 .22) Of course, we could just as well have substituted for C2 in the equation for C3 . This would give (5.23) with (5.24) Thus, we would get the same value for r,;2 verifying that we have found a normal mode . This is the one collective mode normal to the TDHF path. Of course, the matrix equation above also has the trivial solution, X2 = J2=0 meaning C2=C8 =0 . This also means that A=O and the constraint ~h 10 C0 =0 fl becomes trivial. Now we can look back to our equations for the individual Cna and try noncollective perturbations for which cACJ is dependent on n and ~cn 0 =0 n with cnCJJl!O. With C"=O our equations for cna become uncoupled, i.e. singleparticle like. The equations are (5.25) 1 .... 1 . A ao'=hCJa'-6ar1h11-Vcruo'=A a'a (5.26) B' 11o'=-Vfla'l1 =B'u'a. Making the periodic assumption and separating terms as before gives modes with frequencies r.> given by - 116 - (5.27) where (5.28) and I B'I A' A'• E=det B'• (5.29) Thus, we get two values for (,)2 for the two normal modes of noncollective oscillati.ens. However, the exact Hamiltonian H, which treats all particles identically, would not mix these states with any collective state. This means that this instability will not be apparent in a comparison of TDHF with the exact evolution of a collective state. This instability will instead cause doubts about the appropriateness of a collective state description of the many-body wavefunction. - 117 - Section 6; Exact Solution for the SU(3) Model We choose a two-body Hamiltonian which is exactly soluble (Ll 70) (6.1) Yv= vji ' (6.2) V-a =O ' ~ = L: Cln~ Cln; ' n 0,: creates a particle n in state i, and f1ni annihilates a particle n from state i. Applying the usual fermion anti-commutation rules, we get (6.3) and the nine operators ~ are generators for U(3) . Remembering number conservation, L:~=N, this becomes SU(3) . ( In order to find the eigenstates and eigenvalues to solve the problem exactly, we need to find a complete set of basis states. Group theoretic methods have been published (LI 70) which derive this basis set, but for our purposes physical intuition will suffice. First, we must remember that we are only interested in collective states, that is, states with all particles, n, behaving the same. Therefore, we need not worry about each particle and its mean occupation of each level, but only the mean number of particles in each level. This is analogous to Krieger's (KRI 77) using only states with J=*N like the state with N in lowest level. In our case, H will not mix a collective state like this state with any noncollective states since H treats all particles alike . Thus, a natural basis can be written jpq >meaning p particles in the second level, q in the third level and, of course, N-p-q in the first level. To calculate the number of basis states lpq > for a given N, we can pick any number from 0 to N for p and then pick any number from 0 to ,N-p for q. The number in the first level is then determined, N-p-q . Thus, for q we have N-p + 1 choices - 118 - after p is chosen, with N + 1 choices for p. Thus, the number of states NS is NS= f (N+l-p)= (N+1)(N+2) 2 (6.4) p=O that is, the number of choices for p times the average number of choices for q, We define IOO>=TI0.:1 IO> . (6.5) n so that I00> is all N particles in the lowest level. We can write the general basis state jpq>=c(p,q)~ic;gi IOO> (6.6) with c (p ,q) a normalizing constant since ~ is the collective operator to move a particle from state j to state i . To caiculate c (p. q) or any of the other expectation values we will need (<G\i> and <H>), we need to learn how to manipulate these ~·s to arbitrary powers. (Note that since [~i.~i]=O the ~i's in Eq. (6.6) can all be grouped together as can the ~i's.) For q =O these states jpO> are identical to the basis states of (KRI 77) since ~i is his J +· We know that ~IOO>=O unless j=l and G 11 jOO>=NIOO>. Therefore, to evaluate ~; jpq >, we can commute the ~ through G2 i p times, and through G3 i q times, collecting G21 's and G:n's to the left and moving the other ~ 's toward the right. For instance, to verity the meaning o! the states above, we need to show G22 jpq >=p jpq > and Gss jpq >=q Iqp > . From the commutation relation derived earlier we get (6.7) and therefore (6.8) Similarly, ~ 119- (6.9) so Gu lpq >=(N-p -q) lpq > (6.10) For ~ 2 we get (6 .11) and (6 .12) Therefore, (6 .13) Similar relations show that [Gss.~1 ]=O (6.14) [Gss.G3iJ=qGg1 and Gss jpq >=q Jpq > (6 .15) Thus, we have chosen our basis to be eigenstates of these number operators Ga;. Since Ga: lpq > was proportional to lpq > we didn't need to evaluate c (p ,q) yet. However, for the other elements of ~· we need to know c(p,q); for instance ~ 1 lpq>= c(p,q) lp+l,q> (6 .16) c(p+l,q) Before finishing these calculations, we must evaluate all the other ~ jpq >. By a similar repeated application of the commutation rule in Eq. (6.3). we calculate the following commutators: [G11.G11 ]=-p~1 (Gu.G§1 ]=-q~1 [~.Gti]=p~l [~.G§1]=0 [Gss.G§1 ]=qGg1 [G12.G1.J=p~1 (G11-G22)-p(p-l)G11- 1 · [Gss.G11]=0 1 ( G13,G§1 ]=qG'g} 1(Gu-Gss)-q (q-1)Gg} 1 (G1s.G11 ]=-pG11 1 ~ [Gs2.G!1 ]=p~1- 1 Gs1 [G12.G!1 ]=-qG!l 1 Gs2 [~.GK1 ]=q~1G!1 1 (6.17) - 120 [~lo~i)=O [ G~n. Ggi)=O [Galo~1 ]=O [~1.GgiJ=O (~3.~1]=0 [ Gs2. Ggi)=O These relations then give net results for the action on the states lpq > of Gu IFq >=(N-p -q) jpq > G22 IFq >=p jpq > Gas jpq >=q lpq > (6.18) - c(p-l,q}p c (p q) ( N-q- (p -1 )) Ip-1,q> G12 Ipq>I G21lpq>= c ~(p+.f,q)) lp+l.q> G13lpq>= c~~:.1) q(N-(q-1)-p) jp.q-1> G jpq>= c(p,q) jp,q+1> 31 c (p ,q +1) G23 lpq>= c(p,q) qlp+l.q-1> c(p+l,q-1) G32 jpq>= c(p.q) p lp-1 q+l> c(p-1 ,q+l) I We can calculate c (p ,q) by a repeated application of these rules to evaluate <OOI Gf3~2~1Gg1 IOO>. Thus p!(N-q)! Gf2Gt1G§1 lOO>= (N-p-q}! G81 IOO> (6.19) A similar relation for Gf3 resulting in a factor (Z~\! gives N!p!q! <00 I GfaG'1'2 Gl1 Gg1 I00>= (N-p -q )! (6.20) Substituting into the previous Eqs. (6.18), the net result is G12 lpq >=-..Ip (N-q-p + 1) lp-1,q > (6.21) Gislpq >=-..lq(N-q-p +1) lp,q-1> ~i lpq >=-..l(P +l)(N-p -q) IP+ l ,q > Gs1 lpq >=-..l(q + 1)(N-p-q) IF ,q + 1> G23jpq>=-..l(p+1)q lp+l,q-1> G32jpq>=-..l(q+l)p lp-1.q+l> Now we can immediately write down the one-body density matrix elements <p'q' IPv lpq>=i;-<p'q' I G"' lpq> (6.22) - 121 for each ij by writing the above factors and appropriate delta functions . By setting q=O andp=l/2N+m, G,a 1 agrees with (KRI 77)'s J+ and G22 agrees with his L;flp+flp, although there is a misprint in his Eq. 2.12. We must look to his Eq. 2.14 p to kee that Eq. 2.12 should be ( 1 N z+ m ~Omm· . Now we can calculate expectation values of H and thus, eigenvalues and eigenstates of H. We need GJ acting on Jpq >; that is, two applications of the rules in Eq. (6.21) . Thus, with e 1 =-e, e2 =0, e 3 =+e and Vv = f(l-oii ), the result is <p'q' 1Hjpq>=e(-N+p+2q)6pp·Oqq· +~ (6 .23) where A=v'p (p-l)(N-p-q +l)(N-p-q +2)6p-2,p' Oqq' (6 .24) +...J(p + l)(p +2)(N-p-q )(N-p-q-1)6p+2.p' Oqq' +...Jq (q-l)(N-p-q + 1)(N-p -q +2)6pp· 6q- 2 ,q· +...J(q + l)(q +2)(N-p-q )(N-q -p-l)opp' 6q+ 2,q· +...J(p +l)(p +2)q (q -l)op+2 .p · oq-24 +...J(q + l)(q +2)p (p-l)op-2,p' 6q+24 H is hermitian, in fact, with the choice made for Vi.i and phase of c (p ,q), it is real and symmetric . To make this more apparent, we can rewrite some of the terms, for instance p +2=p' as p =p' -2. Now switching p with p' and q with q ' obviously leaves A unchanged. For any given number of particles N we can set up the complete basis states, write down the matrix elements of H and then diagonalize H to find its eigenvalues (real) and its eigenvectors (real and orthogonal) . For even a 4 moderately large N the dimensions of H will be very large (greater than ~ ) . However, H connects only states with Ap=0,+-2 and Aq=0,±2 which makes the problem easier. First, we group states with p ,q even; p even, q odd; p odd, q even; and p ,q odd. This means that <H> becomes block diagonal containing 4 - 122 blocks which can be diagonalized separately, each approximately 1/16 the size of the full <H>. Regardless of the size of N, there can be no more than 7 nonzero elements in any row or column - all quite close to the diagonal for many reasonable schemes of ordering the states. For this kind of sparse matrix, ' the standard routines for eigenvalues and eigenvectors of a matrix work very well and very quickly. Eigenvalues are labeled e A and eigenvectors I>..> . An arbitrary state can be expanded l'!Jt(O)>=~ l>-><>..l,P(O)> (6.25) A and its time dependence is straightforward. l,P(t )>=~I >..>exp(-ie>.t)<>.. !if(O)> A Now for any given wavefunction, we can follow its evolution exactly. (6.26) - 123 - Section 7: How to Compare Exact and TDHF Calculations Our purpose is, of course, to compare with the Hartree-Fock approximation. At this point, we can compare our ground state energy, the smallest eA, with the eqergy of the Hartree-Fock state (and perhaps higher states to periodic TDHF solutions if these can be found), but we want to compare the evolution of an arbitrary TDHF state to its exact evolution .. To do this, we must expand an arbi- trary TDHF state in terms of these eigenvectors and, secondly, we must identify the 4 parameters characterizing a TDHF state from an arbitrary state. Although the eventual goal is to expand a TDHF state 11/1> in terms of the eigenstates IX>; these are not known yet and are different for each N,V different . Therefore, we expand 11/1> in our basis states lpq >; lv>=~lpq><pq 11/1> (7.1) pq Our general TDHF state was defined as (7.2) a general determinantal state, where the states created by a are linear combinations of the three levels, an obvious set of basis states . Thus, using the defintions of a and lpq >. we get <pq 11/1>=<0ITICl.n1Gfo~2fl( f;Ai'iaJi) I O>c (p ,q) n (7.3) n' \=l We need to expand IT~Ai'ia;., . Thus, we get a general term of the structure n' i (A{1)N-n-m(A{2)n(At5)m TI (N-n-m) Cl.nt1fICI.nt2fICJ.nt3 , (n) (7.4) (m) where the number in parantheses under the product is the nwnber of n' 's in the product. The difierent choices or particles for each level give terms of identical structure. First, we need to choose from N particles N-n -m to be in the first level; this gives the standard binomial coefficient (N --n _:;! ( n +m) ! . Then - 124 - with n+m particles left, we need to choose n for the second level - (n+m)! n!m! choices. The net result is then <nq I"''>= ~ r Y' + )n(A + )m N! I n":'m (A+11 )N--n-m(A 12 13 (N-n-m)!n!m! (7.5) n+m<N where /=<OIIT1Zn1G¥sGuaP IT aJ1ITU.:.2IlaJs IO> n (N--n-m) (n) (m) (7.6) G12= ~aJ· i!Zn" 2 commutes with u.n+: 1 and acts on CJ.nt 2 . Terms with n" #-n' give n" zero since !Jn" 2 10>=0. Terms n" =n' leave aJ 1 . Since there are n factors of flnt 2 , there will be n terms left, each with one less G12 , one less CJ.n· 2 and one more aJ 1 • Multiplying all the G12 through gives dpnn!<OITI0n1G¥s TI a.J1TICJ.nts !O> n (N-m) (7.7) (m) The G15 's behave similar to the G12 's leaving (7.8) The net result for <pq 11/1> is <pq 11/1>=( (N-p ~~ )!p!q ! )ff(A;1 )N-P-q(A~1 )1'(Ai1)q (7.9) This can be checked quite simply to be sure <1/1 lpq ><pq 11/1>=1 . Note also that if we constrain q =O in Eq. (7.9), the result agrees with Krieger's results for SU(2) (KRI 77). Now that we can expand any TDHF state in terms of our basis states for the exact solution, we can evolve any TDHF state exactly. We need, however, to compare this exact solution with the TDHF evolution. Since the TDHF states form only a restricted set of many-body wavefunctions, namely determinantal wavefunctions, there will be no way to uniquely match wavefunctions in general. In particular the wavefunctions themselves can not be matched since for a TDHF state 1/1( t) - 125 - <pq 11/t(t)>rt ~ <pq IX>exp(-ie>.t)<Xlp'q'><p'q' 11/t(D)> (7.10) '>qJ'q' The TDHF state (or its expansion parameters) depends only on the 4 parameters used to specify An 1. The one-body density matrix Pii =-}<Gt;> for the TDHF state also depends only on these 4 parameters piJ=A0 A/1 . We use this one-body density matrix to compare to the exact evolution. Of course, the exact <GlJ > depends on more than 4 parameters. A general 3 x 3 complex matrix would require 18 parameters, an hermitian matrix half of this or 9 parameters, and since the Trp=l. we are left with 8 parameters for<~>. This number is reduced to 4 for the TDHF state by the condition that for determinantal wavefunctions p2=p. Thus, there is ambiguity in determining the relevant 4 parameters for comparison. We also want to know Tr(p2 -p) for the exact state since the deviation from the TDHF value of 0 gives us an idea of the deviation of the exact wavefunction from a determinantal wavefunction and an idea of the ambiguity involved in matching p to determine the 4 parameters. Although the choice of which elements to match may vary in specific cases, generally, the most physically relevant property of these states is the mean occupation of the various levels. Therefore, we want to use the diagonal elements Pv. to calculate the magnitudes of the An 1 and thus, the parameters ~ 1 and "'2· (We can also check on some of the numerics involved in the calculation by checking Trp=l.) This means that in general, the magnitude of the off- diagonal elements will not match. To compare the phases relative to the first level, we use p 12 and p 13 to calculate the phases 'lf 1 and 1f2 , but now the phase of p 23 will not match 'lf 1 --1/12 . Therefore, we define 1'1 (exa.ct )=cos-1{v<G 11 >/ N) 1 "'2(ezact )=tan- ( v<Gs:i>/ <G22>) 1' 1 (e.:m.ct )=~os- 1 (Re[ <G 12 >/ <G21 >])xSi.gn (Im[ <G 12 >) (7.11) - 126 - •1J12(exact )=*cos-1 (Re[ <G13>/ <Gs 1>])xSign (Im[ <G13>]) We also check that (7.12) Since there are many exact states which result in these same 4 pararneters, the prescription for finding <pq l'lf> can be thought of as finding the unique exact state which gives these 4 parameters and the auxilary condition p 2 -p=O . (This state must be unique since we find a unique expansion.) We also find this prescription useful to compare our lowest energy eigenvector to the Hartree-Fock state. Since H mixes states with l:!.p ,t!.q =O, ±2 and G mixes states with l:!.p ,l:!.q =0,±1 only the diagonal elements <~> will be nonzero for any eigenvector. Thus, the phases of the off-diagonal elements are zero just like static HF solutions. We compare ,,. 1 and ,,.2 by the prescription in Eq. (7 .11) and also compare Tr(p 2 -p) . Tr(p 2 -p)(exact)=p~1 -pi2 -p~3 -1 since pis diagonal. - 127 Section 8: TDRPA Results For orientation, contours of the energy surface (Eq. ( 4.19) with 1/; 1 =1/;2=0) are plotted in Figures 8-13 for values of x= .5. 2.5. 5, 10, 20 and 100 (attractive forces) . The surface is relatively fiat and featureless for x= .5 but as x increases, it develops rather deep valleys and quite a deep minimum. Static HF calculations can be performed analytically (HOL 74). yielding the position and depth of the minimum, as well as saddle points and maxima. For ~ 1. the minimum energy lies along the line 'jj 1 =O (all particles in the lowest level) and has a value E/nE=l. When x increases past 1, the strength of the attraction between different levels becomes strong enough to decrease the energy by beginning to populate the second level. The position and energy of the minimum are then given by 'jjl =~os- 1 ( ~ x 1 'jj2=0 .· (B.1) 1 E/nE=-:i"\x+~-* . At x=3. the potential is strong enough to begin populating the third level also, and the minimum moves off the ~ 2 axis. Its position and energy are then given by 1'1 =sin-1([ 2Tu~3 ]*) ; ~=Micos-1( 3 ) zx-3 (B.2) E/nE=-(t+? . The TDHF calculations can be checked by calculating the evolution near this minimum in energy. The solutions should then be approximately equal to those of the ordinary RPA calculations, which can be solved analytically. Since the time-dependent generalization agrees with ordinary RPA theory for static HF solutions. these positions can be substituted into Eq. (5.22). The solution involves only straightforward but rather messy algebra and trigonometry. We are most interested in ~3 where the energy surface has a more interesting - 128 - topology. In this case, the RPA frequencies are ( ~2=4( f-3 ±[ ~*) £ 3 3 Therefore, we expect an average period (8.3) eT- 4rre - (c.>++CJ_) • (8.4) and beats in the amplitude with a period (8.5) Table 4 shows the two periods and the beat period for several values of X· Figure 14 shows that the TDHF solutions do oscillate with these periods . Table 4 Period x=5 x=lO x=20 x=100 tT+ 0.975 0 .508 0.262 0.054 tT- 1.519 0.611 0.286 0.055 eT,, 2.721 3.035 3 .126 3.142 Figures 15-18 show results for the TDRPA frequencies for several paths to compare with our expectations from the contours of the energy surface . Beginning with a case of extreme topography, x=lOO, the first example is not plotted in this way although the details of the motion are shown in Figure 26. For this path near the minimum, r.>2 is positive and large everywhere. The path is a small elipse and ( ~ 2 varies only between +1.32 x 104 and +1.42 x 104 . Figure 15 £ shows a TDHF path moving down a valley. The frequency is initially negative because the motion is from the side of the valley down into it. However, once the path begins moving along the valley, the frequency becomes positive and large, and it stays positive and large . Figure 16 shows a path which is unstable almost everywhere. At this higher - 129. energy, Lhe path moves closer to the maxima where the peaks flatten and the curvature of the surface is negative . The tew points where the frequency is positive, mo::'.ly 0n the inne1· loop, cross the ridge closer tc.· the minimwn wher::: the curvature of the top of the ridge is positive. Thus , there can be a ridge effect (if ' the ridgP has a saddle shape) since the direction in which stability is being tested changes from across the ridge to along the ridge . The next example, Figure 17, shows a path with more nearly equal stable and unstable regions . The central right region in the vicinity of the minimum, where the curvature of V is positive, has positive frequencies and should be stable . The upper right region shows another ridge effect. The path is unstable while moving along either side of the ridge, but appears stable briefly as the top of the ridge is crossed and the direction of the stability test crosses along the top of the ridge . The lower left region shows interesting behavior where two paths cross the initial rightward path. The inner loop, lying closer to the minimum, is stable but the outer loop isn't. The initial rightward path tests stability in a different direction and is unstable. Therefore, in the vicinity of the inner loop, the surface must have a saddle shape . These features persist at other values of x with less drastic topology as seen in Figure 18 for x=5. Here, the path is unstable while moving back and forth across the valley (at either end), but is stable while moving along the valley. Just as for ordinary RPA theory, the values of r,;2 scale like x2 - on the order of 104 for x= 100 and 25 for x=5. Therefore, the value of r,;2 does, indeed, describe the local curvature of the potential energy surface in the direction normal to the TDHF path and can thereby identify relative minima, ridges, valleys and saddle shapes . In this case, we could easily draw contours of the energy surface; however, in more realistic cases there are far too many parameters to allow the construction of this kind - 130 of energy surface, so the smallest values of (,)2 could prove useful. We next compare these TDRPA calculations with exact calculations, first checking the HF results. Figure 19 shows that the energies of the exact ground states agree quite well for N~15. (Lines are to guide the eye since N is an integer.) However, comparing the positions and Tr(p2 -p) for these states gives poorer agreement. The trace, shown in Figure 20, shows that none of these states looks much like HF (independent particle) states which have Tr=O. However, the important quantity tor comparison is the position ~ 1 .~ . This comparison is made in Figure 21 for x= 5, 10, 20. For x=5, agreement is quite poor. As x increases, the agreement gets much better and also seems to agree for smaller values of N. For x=lOO, both N=10 and 15 give good agreement with the position of the ground state. This is a problem because the method of comparing the evolution of TDHF and exact states is to compare these angular positions. Indeed, for calculations with small x. the exact solution doesn't look like TDHF regardless of the value of ~. Tb.us, in Figure 22, for x=5, the direction of motion seems to be right when r.:}>O and wrong when r.>2 <0, but neither of the exact results agrees well with the TDHF results . (Notice that the TDHF path moves so slowly relative to the exact, that it is drawn for twice as long a time so that it can be seen more clearly.) The exact solutions, of course, depend on the number of particles N and TDHF is accurate only for large N. Therefore, the agreement could be poor because TDHF is a poor approximation to the exact behavior or because N is too small to give this limiting behavior for large N. For the SU(2) case (KRI 77), N=40 particles could be solved exactly since the number of basis states is proportional to N, and the agreement was very good. In our case, with an extra level populated, we would expect to need N>40 to be sure of a similar comparison. However, since the number of basis states is proportional to N 2 , we - 1::11 simply can't run N>15 for any length of time . From the change in the exact results from N=lD to N=15 , the stability parameter seems to work. That is. a pr~dic.:Uon of stability re:'mlts in exact soiutipns which approach the TDHF solution as N increases while a predicticm of instability seems to result in an approach to something different than the TDHF result. The real test, however, comes from increasing x so that the limit ing behavior of the exact solutions for large N can be reached at a smaller value of N (as for the static solutions described earlier) . Thus, we look for solutions near the minimum, which are always predicted to be stable and increase x until the exact solution agrees fairly well with the TDHF solution. Figures 23-26 plot 13- 1 and 13-2 as functions of time for the TDHF solution and two exact solutions N= 10 and N=15 . These show the progression from poor agreement for x=5 to very good agreement for a full period in both 13- 1 and1J 2 for x=lOO withN=15 . A plot of Tr(p 2 -p) in Figure 27 also shows that the exact state in this case is approaching an independent particle state with Tr=O. For an unstable example at x=100, we test the path from Figure 16. In Figures 28-30, plots of both 1J 1 and~ vs. time and a trajectory plot of 1J 1 vs. ~ 2 show that TDHF is a poor approximation to the exact solution. Figure 31 also shows that the Tr(p 2 -p) climbs quickly to a relatively large value, .6, and then levels off - a state far removed from a TDHF state. We also notice that all these results in Figures 28-31 don't vary much from N=10 to N=15 so that these calculations probably are in the "large N regime" where the exact solutions have reached their limiting behavior . - 132Section 9: Summary This TDRPA frequency can be used in a qualitative way to get fairly accurate ideas of the shape of the energy surface, and might, therefore, prove useful in more complex models . It also works well for predicting the reliability of TDHF calculations in extreme cases, with x= 100. However, because of the scaling of the frequency with x. this means the frequency has only proven reliable for I~ $::j100. The motion of the exact solutions as N increases shows that G>2 seems £ to be a good indicator of the reliability of TDHF calculations in general. However, even this relatively simple model can't be solved exactly for a large enough number of particles to permit quantitative testing of this criterion or the development of quantitative ideas about its meaning . The results also suggest another possible use for this frequency. Its magnitude might give a qualitative idea of the number of particles which can be considered large in a given situaUon. -133- APPENDICES - 134 Appendix Al: Derivation of the TDHF Equations L Variation.al Derivation (KE 76) We begin with the definition of the action, 1 =f d.t<vlih:t -H Iv>. (Al. l) and parametrize the many-body wavefunction, v, as a Slater determinant of single-particle wavefunctions,1/1;. (Al.2) where the a; refer to spatial. spin, and isospin nucleon coordinates and A is the total number of nucleons. These single-particle wavefunctions will be orthonormal, <'ljl, l't/t;> = o;,; . (Al.3) The calculation of <v Iv> checks that the many-body wavefunction is normalized. (Al.4) Each determinant has A! terms so the product has A! terms where the coordinates pair identical single-particle wavefunctions, (Al.5) The other terms include at least two mismatched wavefunctions, with <1/l;.11/11 ,.;> = 0 ; and thus, do not contribute. Therefore , <vlv> = 1 . (Al.6) In a similar way, we can rewrite <vlm:t Iv>. As above, this product has (A1) 2 terms. Of these, A! properly match the single-particle wavefunctions and the rest mismatch at least two of them. The derivative of each of the product terms can be expressed as a sum involving derivatives of single-particle wavefunctions, - 1358 TI n. B (Al .7) ;/a.fit . 1/1; = L:i <ihF'l/I; > _ 1/1, ' ' J Tb.us, any of the mismatched terms still have at least one mismatch, with <1/li 11/l;f0\>=0, and do not contribute . Each of the A! matched terms gives L:<li1 lih:t 11/l;>Il <1/li 11/1,> = L:<'l/I; lih:t 11/1;> (Al.8) <irl'ih:t lir>=L:<'l/I; liii.:t 11/1;> (Al.9) i Tb.us, our net result is , ..; j ; ih:t 1/1; · and the variation with respect to 1/1;.will give I We then arrive at the TDHF equations in the form ih.E_~,,. Bt "' = 6 <H > (Al.10) 61/1/ To see that Eq. (Al.10) results in a mean-field formulation, ih:t 1/1; =h'l/I; , with h a one-body operator, we calculate tS~H? for a two-body Hamiltonian. The 1/11 generalization for higher-body potentials is straightforward. A two-body Hamiltonian can be written in second-quantized notation as, * =L;ta.B a.: a.p + L: Vah6 a: a; a4a7 a.P a.fJ'r6 H (Al .11) The greek labels refer to the spatial, spin, and isospin nucleon coordinates; and the annihilation and creation operators satisfy the usual fermion anticommutation rules: fa..1" ,a/~ = laa ,apJ = 0 and ~a: ,a,~ = 6'1.p . (Al .12) Since H is hermitian, t a/1 = t fla = and Va.'16 Vfla67 With ir a Slater determinant as in Eq. (Al.2), = v;6a.B (Al.13) <H> = l; tt/$1/(a.)1/I; (p) *L: a/Ii + ~ij v~ 1/lt(a.)1/l/(f3)[1/l,(-r)1/1;(6)-1/li(o)1/1;(7)] (Al.14) - 136 The variation of <H> with respect to 'if;; gives o<H> o'if;;(r;) - ~tap6;k 6at1/l;(f3) + apj * ~ Va.8?6 (6ik 6at1/l/(f3) + 6;k 6pt1/l/(o..)l x a,8?6ij f1/li ( 1)1/11 ( 6)-1/li ( 6)1/lj b') l . = ~tcp1/lk(f3) +~Vc,8?6[~1/l/(f3)1/11(6)]1/lk(/) , fJ p-y6 ; =~hcp'l/lk ({3) ; (Al . 15) fJ where V are the antisymmetrized elements of V, Va.8? 6 = Va.8? 6- Va,867 : h =t + W : and the mean-field potential, W, a one-body operator, is given by Wap = ~Varp6[~1/l;lr)1/1;(6)] . (Al.16) ?6 ; For >Ir a Slater determinant, we can write the one-body density matrix, (A1.17) and rewrite Eq. (Al.16) as WaP = ~ Ya1P6 P67 · (Al. 18) "'16 We can also rewrite Eq. (Al.14) for <H> in terms of the mean field and density as <H > = ~tap Ppa + }2 ~ Ya,8?6 P"'1a P6p • aP afJ?'6 = 1;[tap + *Wap]Ppa · (Al. 19) a{J Similarly, for a general n-body potential, Va 1 ... anp 1 .. . fJn ; (Al.20) where Ya 1 . .. ~p 1 .. . /Jn is completely antisymmetrized with respect to inter- changes of the cx's and interchanges of the {J's . The variation gives (A1.21) - 137 Thus the mean field, (Al.22) and (Al.23) Il. Derivation by Truncation (KO 79a) This derivation begins with the definitions of the density matrices. Pop = «If Ia/aa I'11> , P~6 (Al .24) =<it Ia;a!aaafJ I'11> • (Al.25) etc. We once again use a two-body Hamiltonian, as in Eq. (Al.11),so * <H> =~tafJP~ + 2: Va.&?'6 pJ~~Jl . afJ afty6 The exact equation for the evolution of the one-body density is ih:t (Al.26) Pop= <'11 J-Halao.+a/aaH I'11> , = ~t715<'1t I[a/aa, a.;a6] lit> 76 * + ~ V7 6c('<'1tl[apaa, a;a/a('a J I'11> 1 (Al.27) 76c(' Using the anti-commutation rules from Eq. (Al.12), we get (Al .28) (Al.29) Therefore, ih:t Pop= ~(ta.6 P6TJ-Pa6 f411) + I * L;('i1o.6c(' Pt~b6-Pt~~6 Vc6fJ(') {Al .30) le(' This equation gives the exact evolution of Pap for a two-body potential in terms of a two-body dynamical quantity pZf.r6 . Of course , higher-body potentials would involve higher-body dynamical quantities . The evolution of p<2 > can be described by a similar equation involving p(3 ) (or higher for higher-body poten- - 138 - tials). The hierarchy terminates at p(A) and is equivalent to the exact Schroedinger equation. We make the independent-particle approximation by saying that the only two-body correlation is the necessary antisymmetry, (A1.31) and we ignore all the higher-body equations. Eq. (Al.31) is true for '11 a Slater determinant since a comparison of Eqs. (A1.26) and (A1.14) shows that PafJ = "E'lfj(f/)1/1; (a) and PA~~P = "E'lft(a)'!fj(p)[1/I, (7)1/i; (o}-'!f,(o)'!f;(/')]. However, it ~ j can also be shown that the assumption of Eq. (Al.31) leads to a Slater determinant wavefunction. Setting o = a and summing over a in Eq. (A1 .31) gives (Al.32) a a a Cll Since the number operator is N ="Eaiiaa, "EPaa =<N> =A the total number of a a particles and <a; Nap> = (A-l)<a;a 11 >. Therefore, p2 = p and p is a projector. This means it has A eigenvalues of l, the rest zero. Therefore, we can write a spectral expansion for Pa.fl , Po.p = t1"j(P)1/l;(a.) : i=l , where the functions 'If; span the space of unit eigenvalues Dap1/I; (p) = 'If; (a) and , are orthonormal "E'lfj(p)y;,(p) =oil'. The above expression for p is then consistent with a Slater determinant many-body wavefunction as in Eq. (Al.2) . Substituting this approximation, Eq. (Al.31), into our exact expression gives ih:t Pap= "L,(ta6P6P-Pa6tlfl) +l,;(Va.6e(P('6Pcfl-Pa.&P{6~c6fJ(} • I le{ =I;[ (t al+ Wa6)P6p-Pal(t lfJ+ W6p)] • I =[h , p ]o.p : (Al.33) where h and W are the same as the previous definitions in Eqs. (Al.15) through - 139 - (Al.18). Similarly for an n-body potential (Al. 34) we get <[a..;a4, v<n>]>= ~! [V4~ · · · °'n/J1 · .. fJ11 PJ:>. · · /J17a2 · .. an (Al.35) -p!:>. ··ll26a1 · · · °'n Val ··· ~7fJ2 ··· /Jn] where V is antisymmetric under permutations of o through the a ' s or equivalently -y through the P's while keeping the a's or P's in order. For Slater determinant wavefunctions p{n) is written as an antisymmetric combination of one-body p's, with n! terms . Thus there are n terms where I' is permuted through the a 's and the a's remain in the same order. Because of the symmetry of V these all give the same result. The (n-1)! permutations of the a's for each position of -y form the completely antisymmetric matrix elements V. Therefore, <[a..;a4, v<n >]>= (n .:.l)! ~ ( V6a2 · · · °'n/Ji · · · fJ11PfJ17P/Jt!'-2 · · ' P/Jn °'n 3 -P6a1 Val · ·· °'n7fJ2 · · · fJ PfJ~2 ' · ' P/Jn°'n) =I; ( W6fJ1PP1r-P6fJ1WfJ{r) (Al.36) ' • 111 with Win agreement with Eq. (Al.22). Therefore, (Al.37) ID. Hartree.Potential and J"ock Potential The traditional Hartree and Fock terms can be calculated by evaluating this mean field, W, for a local, Gaililean-invariant two-body potential, V, dependent only on the relative position of the particles. ) (r1+r2 rs+r4l16r (...... ( ... ... ) (...... ) ( .......... Yah6-.Vr 1 ,r 2 ,r 5 ,r 4 =6 1 -r 2 -r 5 -r 4 1 v r 1 -r 2 2 2 Ignoring the spins of the particles, Eq. (Al.18) gives a mean field (A1.3B) - 140 - (Al.39) In this first term, the two o functions require r 1 =r 3 and r 2 =r 4 . Thus, the first term, the direct term, gives Un(r1,-ra) = o(r1-ra)f a:r2p(r2.r2)v(r1-T2) . =6(r1-Ta)fdr2p(r2)v(r1-r2) , (Al .40) the local Hartree term. In the second term, the exchange term, the two 6 functions require r 1 =r4 (Al .41) the nonlocal Fock term. IV. Conservation Laws We want to show that for a one-body operator the TDHF equations conserve any expectation values conserved by the exact H, but that this is not true for many-body operators. The time evolution of <0> for the TDHF state is given by Ui~O>=~[ 6<0> o<H> - 6<0> o<H> at jo. 01/J; (a) o'ljlj(a..) 61/lf(a..) 61/1; (a) (Al. 42) Using the results from Part I, for a general n-body operator we get 6 ;n8t <O>= ~ ('l/l/(o)04-r1f7 4·1/J;(o')-1/lj(o)H6'70-r6'1/l;(o')) j76cS' =~cSd' [0 ,H ]cSctPcS'cS =<[D,H]> (Al.43) where 0 and H are one-body operators defined by Eq. (Al.22) . Of course, the exact evolution would be B ih'Bf<D>=<[ O,H]> (Al .44) Therefore, we must compare these two commutators from Eqs . (Al.43-44) . For a one-body operator 0 and an n-body H, the expectation value of the commutator is - 141 - <[ O,H]>= L0f'6<[a.;-a6,H]> (Al .45) f'6 Using Eq. (Al.36) from Part II this gives (Al.46) which can be rewritten {Al.47) Since O=O for a one-body operator, the TDHF evolution with the one-body h=H conserves all expectations of one-body operators which are conserved by the exact n-body Hamiltonian, in addition to the energy and normalization. However, it will not conserve expectations of many-body operators. For instance, with both 0 and H two-body operators and Ji the antisymmetrized elements of H the commutator is 1 X' .... [ 0 ,H] = 4 L.i Oo.1P17161 (H-r1P27e62a;1 a;1 a61 a;2a62a-Y2 (Al.48) -Jio.lifl20.162a:_2a;2a6iil1 a61 ~i where the sum is over all repeated indices. Collecting all a+•s to the left yields two-body and three-body pieces. The three-body pieces can be paired by defining an antisymrnetrized 0. The expectation values for a Slater determinant many-body wavefunction are antisymmetric combinations of one-body density matrices. For the two-body pieces, the two pairs of terms can be combined by changing 0 to 0. For each three-body piece, the terms form two groups : 4 terms .... .... have only one p with one index from each 0 and H, while 2 terms have all 3 p with one index from each. Within each group the terms are equal because of the symmetries of 0 and H. The four terms with only one mixed p give (Al.49) - 142 - Therefore, the difference in expectation values is (Al .50) - 143 - Appendix A2: <H> for the !:kyrme Potential I. lntr .xluction and Kinetic Energy Term We <H> = want to evaluate the energy functional. Hs (r) , defined by Jd.r H. (r), for a Skyrrne potential and a Slater determinant wavefunc- tion. For a Slater determinant wavefunction <H> =~ <i IT Ii> + 21 ~<ij I v<2>l ij-ji> \i (A2.1) i + .L ~<ijk I y<3 > lijk +iki+kiJ-jilc -kji-1lcJ > . 6 ijk where T is the one-body, kinetic energy operator, and v< 2> (or vC 3>) is the two(or three-) body potential. In coordinate space, T(r) =- ! V2 or (A2. 2) Similarly in coordinate space, the Skyrme potentials are : t VC 2>(r ,r') = t 0 (1 + Xo Pa) 6(r) o(r') - T[v; o(r) o(r') + v;. 6(r ') o(f )] (A2.3) + t 2 [Vr o(r)l[Vr·6(r ')] + iWa (a 1 + C12) ·[Vr o(r)] x (Vr·o(r ')] , or v<2>(r1.r{Z.r3,r4) =o(i? - R') v<2>(r,r') . .. .. .. ~ .. .. .. where r =r 1-r 2 , r' =r 3 -r 4 , R = r1+r2 2 .. , and R' = rs+r4 2 and (A2.4) or In these expressions: Pa=)?(l+a1 ·a2 ), the spin-exchange operator, and 0-1 . a2 are Pauli spin matrices . The rewriting of the expressions above is to make Gaililean invariance ( conservation of total momentum ) manifest ; that is, the potentials only depend on the relative coordinates and do not change the - 144- centers of mass . A general state, Ii>. is l9'1 (r)>. a four component spinor in spin, isospin space. (A2.5) where 9'i(r,a,,q,) is a wavefunction characterized by position r, spin a-;., and isospin qi. and Iai >. Iq, > are each two component spinors, in spin and isospin, respectively. Since each single-particle state is either a proton or a neutron. specifying Ii> also specifies its isospin q, (either p or n ) and the spinors Iq,: > are eigeruitates of the isospin operator, Qz . ( The more traditional 'T is being used for the kinetic energy density.) Therefore. the spinors are orthogonal. (A2.6) We consider only spin-saturated systems, so th.at for each Ji> as above there is a corresponding state j i '> with opposite spin. but with the same isospin. We also neglect the spin-orbit force. so the state Ji'> has the same spatial orbital as the state Ii>. (A2.7) where <~ I0-1~> =- <-~ I0-1-~> and ip,(r,-~.q,) =9'i(r.~.q,) . (A2.B) These two states. l~ > and 1-0'.\ >. form a complete set of spin states, so I;t (~)=I(~)+/(~) (A2.9) Ac Therefore, I;<a'.\ l<fla,> =0 I (A2.10) Ac I;la'.\><a,1=1, Ac and I;<a'.\ I~> =2 . Ac Since ip,(r,ai,q,) is independent of the spin ai, it can be rewritten 9'i(r ,q-;.) and - 145 - the sum over i can be replaced by a sum over i and a.i . The ability to sum over the spinors, lcxi>. without weighting factors from rpi(r,a.i,qi) is a crucial simplification in dealing with spin-saturated systems . The one-body piece, ~<i IT li>,can be evaluated as follows . i ~<i l T Ii> = ~fd.T1f dr2rpt(r2.qi )<ai I <qi I 6(r i-f2)(- 2-w. vi1)rpi (r1 .qi) Ia.i> Iqi> '~ ' m =ia, ~Jdr1rpt(r1.qi)(- : v; )rpi(r1.qi:)<ai la.i><qi lqi> m 1 =2~Jdr1rpt(r1,qi:)(- 21f.m v;1)r;i(r ,qi) i =2Jdr 2! 41v~i(r , q") 12 \ =f dr:m r(r) (A2 .11) where 'T'q(r)=2~ 1Vr;i(f,qi)l 2 and 'T'=Tp+'Tn . (A2 .12) i with q,=q D. Two-Body Terms Since the two-body potential has been written in terms of the relative coor- .. .. dinates, rand r', and the center of mass coordinates, R and R'; the expectation value is calculated by setting r 1 =R+~. r2 =R-W, r5 =R'+~·. and f d.r d.r dr dr becomes f drd.r'cJRdk . r 4 =R ' -~' . Therefore, the integral The integral J dR'6(R-R') can be done by replacing R' by R in the rest of the 1 2 5 4 expression. Since none of the terms in v<2 > involve any isospin operators. the inner product of the isospin spinors can be done very easily for any term. The direct terms <q;. Iq; ><q; Iqi> give <q, Iqi ><q; Iq; > = 1, and the exchange terms give =oq,qr The two-body contribution is evaluated term by term, neglecting the spin-orbit term. The t 0 term from Eq. (A2.3) gives - 146- V0 ='h"Efdrdr'ci.Rr;;t(R+1/il',q,)rp/(R-1/il'.qi )o(r)o(r') (A2 .13) ij "E to( <Cl.t Iai><a; Iai >+x <ai I<ai IP Ia;.> Iai> )rp;.(R +}2-r ,qi)rpi(R-}2T ,q1 ) [°"f.aJ 6q,q1 "E t 0 ( <ai I ai ><ail a, >+xo<ail <ai !Pa Ia;> l ai> )rpi (R-}2T .qi )rpi(R +}W,q.Jl 0 11 °"f.aJ The simplest way to evaluate the contribution from P 17 is simply to regard it as a (A2 .14) and Since "E <ai Ia;.><ai Ia;>= 2x2 = 4 and "E <ai: Iai ><a; Ia;.>= "E<ai Iat> = 2, °'i a.i tl.t. aJ tl.t. the contribution from the t 0 term is Vo='h"Ef d.Tdr'dRr;;t(R+1/il' ,q,)rp/(R-1/il' .q; ) fto( 4+2xo)rtJi (R +1/il ,qi )<Pi (R-1/il ,qi) 1 -6q,q/ o(2+4xo)<Pi (R-'hf ,q; )rp; (R +Y.f ,qi) ]o(r)o (r') ij ='hl:ifd.i?rpt(R .q;,)rp/(.R .q; )[to(4+2xo)'Pi (R .qi)<P;(.R .q;) . (A2 .15) ij -c5q,q/o(2+4xo)rp.JR,q; )rp;(R ,qi)] . With Pq (R)= "E rpt(R ,q, )<Pi (R .qi)=Z "E rpt(R .q;,)rpi (R .qi) i with q,=q (A2 .16) i wilh q~=q and p(R)=Pp (R)+Pn (R) . Eq. (A2 .15) can be written Vo='hf &{to(4+2xo)(}?p) 2 -to(2+4xo)[(}?pp ) 2 +(¥1Pn )2 ]] = , JdH~oJ(1+*xo)(pp+Pn) 2 -('h+xo)(pi+p~)] . =J aR~o~(l-x.o)(pj+p~)+(Z+xo)PpPn] . (A2 . l 7) -147 - The ti term from Eq. (A2.3) gives Vi= - ~ ~f ctrdr'dRSDt(R+W'.qi)Soj(R-W',q; ){[v;o(r)]o(r')+[v;.o(r')]o(r) (49',(R+}W,qi)Soj(.R-W .q;)-oq.q1 2rpi (.R--W.q;)i;P;(R+W.qi)) . )x (A2.1B) Since r and r' are independent variables, integrating by parts twice for each v; and v~. gives Vi=- i ~f ardr'ctR6(r)6(r')[9't(R+W'.qi)9'j(R--W',q;) ~v;[rpi(.R+W.qi)rp;(.R-W.q; )J -26,.q1v;[i;P, (.R-~ .q; )rp;(R+*7' ,qi)J)+v;.[rp/(R+*7''.qi)rp/(i?-~· .q; )Jx ~,(i?+}W,q,)9';(A-~ .q; )-20,,,1 1iPi(.R-~.q; )9';(A+~.qi) )] : (A2.19) where V2[i;Pi9';] is shorthand for V2[9',9';]=(V2rpi)SD;+2(VSD,)·(Vi;P;)+i;P,;(V2i;P;). and VrSDi (R±*7' ,q )=±*VSD, (R ±~.q ). Therefore the contribution from t 1 is t1 J .., 1 .. .. ... .... .. .. .. 2 .. Vi=-4ij d.R 4( 'i1[v2i;Pi (R ,qi)SD; (R .q; )-2(V9'i: (R ,q,))· (Vi;Pj (R,q; ))+Soi (R,q,)V i;P;(R .q; )] -269.q1 ( 1...r ~[Vlrp,(R .q; )rp; (R ,qi)-2(VSDi (R .q; ))· (Vrp;(R ,q,))+i;P, (R .qj )V2rpj(A .qi)]] ne;! • .... • .... ... • ... .. • .. • .. ).,.., • .. + 4l( v-rp, (R ,q,))rpj(R .qj )-2(Vrp, (R ,q,))-(Vrpj (R .q; ))+i;P, (R ,qi Y-9'j (R .q; )]x [ 49'i (R ,q.Jrpj(.R ,q1)-261, 111 1iPi(A .q; )rp; (R .q1)]) . (A2.20) The sums over i and j can be done, to express this result in terms of densities, by defining a current density 1"t9 =Im " £J ...... ~ .... j rpi(r ,q,)V9'-c (r ,q,) = i1 " £J (rp, v9',-rp, Vrpi (~ wUh ( with. ~~ ~~ and rearranging the terms in Eq. (A2.20) as follows: (A2 .21) -148- (A2.22) f t1 ~fr[ •( .. )"2 .. "2 • .. .. • .. .. Y1=-4 dRL...ij rpi\R.q, v-rp,(R ,q,)+v-rp, (R ,q,)rp, (R ,q,)]rpi(R .q; )rp; (R .q; )x2 i,i -2[rpt(.R ,q,)Vrp,(R ,q,)-rpj(R,qs )Vrp; c.R .qi )+Vrpt(.R ,q, )rp,(R ,q,)-Vrpj(R,qs )rp;(R .q; >1) -*6,. 111 ([rpt(.R ,q,)V2rp,(R ,qi )+V2rpt(.R,q,)rp,(R ,qi )]rpj(R .qi )rp;(R .qi)x2 -2[rpt(R ,q,)Vrp, (R ,qJ )-rpj(R .q1 )VrpJ(R .q")+rp, (R .q; )vrpt(.R ,q,)-rp; (R ,q,.)Vrpj(R .q;) 1]) where the extra factors of two come from terms with i and j switched. With p=2I;,t<r ,q,)rp,(r ,q,). f. ~p=2I;(vrpt(r,q,)rp,(r,q,)+rpt(r,q,)Vrp,(r,q,)) . and (A2.23) f. V2p=2I;[v2rpt(r ,q,)rp, (r .q,)+2vrp,,r .q,)-rp,(r ,q,)+rpt(r,q,)V2rp, (r ,q, )J Therefore, ' (A2.24) ~,t(r,q,)vrp,(r,q") = t<vp+2iJ) . The contribution from the t 1 term becomes v,=- ~ J ~[(J!V"p-,-)p-2( {6 (Vp+?fi.f)•+ 1~ (Vp-21:Jl21J ~[<m2P11-'T1)P11-2[ 116 (vp,+2ij~)2+ 1~ (Vp11-2il11>21J] . =-!!._fdR[(YiPV2p-pT--!--IVp 12 +;2) 4 4 (A2.25) ~<*1?11V2P11-P,,..,-tivpq 12+;l>J ., Doing an integration by parts on l~P 12 • 8Ild combining terms gives, 3 ~ +Pn 'flpp)+ppTn +Pn1"p-2Jp"Jn] JdR..[t1t-£-~V2Pn +1}t-~p V'lpp+Pn V2pn)+ppTp+PnTn -JJ-j~]l Yi= The contribution from the t 2 term is "t (A2.26) . - 149 - V2=1h~fdrdr'd.i?rp;."(R+1h-;t',q;,)t2Vr6(r)·Vr·6(r')x i.j [~<a;, la;,><a; la;>rp;.(R+1h-;t,q;,)rp;(R-1h.;t.q;) ll(.llj -6q(qi ~<a;. la:;><cx; la;.>rp;.(R-W .q;);o;.(R+W.q;,)] , ll(O.i =1h~f d.Td.T'ctRt26(r)6(r')Vr{Sot(R+W'.qdrp/(R-W'.q; )l '! [4Vr(<,O;.(R+1hT .q;,)<,O; (R-*1 .q; ))-26q,qi Vr (rpi (R-1h.;t .q; )rp; (R+W .q;,))] ' = t 2 f cti?~1h[Vrp;."(R ,q.;.)rp/(R .q; )-rp;."(R,q,)Vrp/(R .q; )l 2 ii ~ 4(Vrp, (R,q,)r,o;(R .q; )-rp;.(R .q;,)Vr,o; (R .q;) [ = ~ J dil[2pT-Jl(Vp +2if) (iip-2if) +*~ 2p, T. -J>(VP, + f?:if,) (VP, -2ij~) iJ f = ctR[it1h(pn V2Pp +pp V2pn)+2(pn 'Tp +pp'Tn -2,h Jn)] (A2 .27) l 3t2 ~272\ +1if{1h(pp V2pp +pn V2pn)+2(ppip +pnin -Jp]nJ] fil Three-Body Term. and Density-Dependent Two-Body Term The three-body contribution can be written quite easily since t 96(ri---r2 )6(r 2 -T 3 ) requires all three particles in the same place . This means that three nucleons must be two like nucleons, one spin up and one spin down, and an unlike nucleon in either spin state. The density in a particular spin state is }',p . Therefore, the total three-body contribution is fs[(p; )2Pp +(pi )2Pn ]= ~P +pn)PpPn The full calculation verifies this result. (A2 .2B) - 150 - Adding all these contributions gives ':n "r+*to[}« 1-xo)(p:+p~)+(Z+xo)PpPn]+ ~PnPp Hs(r)= 2 + (A2.29) 1 (t2-3t i)(pn V2pP +pp V2pn)+ ; (t2-t i)(pp V2Pn +pn V2pp) 16 2 ~ ~ 1 ) (pp'rp +pn Tn -Jp -:-2 72 + .!..J 4't1+t2 )(pn Tp+Pp'rn -2Jp "Jn)+ f}{t1+3t2 -Jn) The density-dependent two-body potential v2 =~1 +P,,)o(r-i'')o(r)p(R)o(R-R') (A2.30) is Gaililean-invariant and local with zero range. Its expectation value <v 2 > is (A2.31) .. ts .... •(R... ,q,)rp;•(R... .q; )(6rp, (R... ,qi )rp; (R... .q; )-6oq,q rp, (R... .q; )rp; (R... ,qJ) *t"fd.Rs--P(R)rp, 1 since (A2.32) Therefore, <v2>= !.1_f ctRp(R)[( ~2-2:(Pq(R) )2] 2 2 q 2 tsf ... =4 d.RPPnPp It is therefore equivalent to the three-body potential above. (A2.33) - 151 - Appendix A3: Spatial Discretization of. the Energy Functional I. Introduction and Zero-Range Terms We solve the TDHF equations by finite difierence methods in coordinate space. We obtain the discrete representation of the TDHF Hamiltonian by the variation of a discrete approximation to the energy functional with respect to the values of the single-particle wavefunctions at the mesh points . This procedure guarantees a hermitian approximation to the single-particle Hamiltonian and would yield exact conservation of the discrete approximation to the energy if we could use an evolution operator which is exact in time. We split the energy functional, H, into zero-range, Hz; kinetic, Hx ; nonlocal, HNL; Yukawa, Hy; and Coulomb, He . energies. (A3. l) H=Hz+Hx+HNL+Hy+Hc Our TDHF equations are "( ' . k) 6H '11.:t YJ>.q (i ,j ,k )=hv (i,j ,k) where ""'111,·1· =~,,/,.(' ' k) ·1 • (A3.2) u r>iq .,, We use a uniformly-spaced cartesian mesh within a rectangular box. Vanishing boundary conditions are imposed outside of this box and the separable approximation allows some simplification. We use a relatively large mesh spacing (typically 1 fm in all directions) because of computing time and storage consideralions; so we need care to obtain an accurate representation of the energy functional. We begin with the wavefunctions at each mesh point TJ>iq (i ,j ,k ) . Since 1/1>.q =O outside the box, we must be sure nothing gets too close to the edges. From these wavefunctions we calculate the proton and neutron densities at each mesh point p, (i,j ,k )= ~ 11>.q (i,j ,k) 12 >. (A3.3) - 152- Previous work has shown that the zero-range terms, Hz, involving t 0 and t 5 are accurately discretized by a simple trapezoidal approximation to the integrals. (A3.4) and f d.rp:(r)p 9 .(11)= ~p:(i .; ,k )p 9.(i .j ,k )(.6.x )s (A3.5) fJA; Therefore, the discretized expression for Hz is Hz=~~o[*(l-'X.o)(p:(i,j ,k)+p:(i ,; ,k ))+(2+'X.o)Pp(i,j ,k )Pn (i,j ,k )] ii" + ~p(i,j ,k )+pn(i,j ,k ))pp(i,j ,k )Pn(i,j ,k )] (A3.6) Il. Kinetic Energy Term The kinetic energy term. Hx. is more difficult because of the gradient operator. The treatment here is described in (FL 78) . Instead of first approximating -r= IVV; j 2 at each mesh point and then summing the contributions, it is more accurate to directly approximate the integral of T over each interval and sum these contributions . The difference in accuracy is shown by Taylor's expansions in the above reference. This method consists of approximating the integral ifi+(») * I d:T(.:r:>=*(-r++,-><.6.x) (A3.7) Ir(-(») where (A3.8) with b - 2 1 o- ""3 (Az) and b.t= ( -1 L (A3.9) 1 vs ±27 (Az) This approximation is the simplest of the most accurate group of formulas - 153- involving only the points at i,i±l. We can see that this is an application of Gauss' method, using the two points at i ± .Js with equal weightings = }2, from the Taylor expansion for -r. b~(i-l)+b 01f(i)+b:1:'¢'(i+l)=±'¢''(i)-(~:;/ ~(i)± (~) 'f/13+ ... 2 =±'!f'(i- (A3 .10) .Js) with error=order ((~) ) 3 Thus in 3D, Hx= ! ~*[ (A3 .11) -ri(i,j ,k )+-ri(i,j ,k )+-r:(i,j ,k )+-r;(i ,j .k )+-rt(i ,j ,k )+-r;(i ,j ,k )] where ill. Nonlocal Terms The nonlocal term HNL is the new feature, so its discretization is described in more detail. The fundamental quantity of interest is again the one- dimensional integral over an interval of±(~). Just as for the kinetic energy, we first discretize the integral for each component of -r and j in that direction and then perform a simple sum, like for the zero-range terms, in the directions not involving the differential operator. We look at the general case Jdx/g. The Taylor expansion for the integral is then * J =1u=<~>[1u+ <t;' ~+(~) S(-(~) . 2 (J 2g+21·g·+1g2)+. ·. Jl:r:=:r:, (A3 .13) Once again we want to use a formula involving only the 3 points i,i±1 . The values of r and g at these points can also be expanded. We will find that the best formula is similar to that for -r. We define I :1:R!.J (%;. ± (':s) )and g ±~g (%;. ± (~) ~ - 154- and approximate the integral (A3.14) If we try a simpler formula involving only one combination off's and g's, our error will be order ((Ax)2). (b-/ (i-l)+bof (i)+b +/ (i+ l))(a._g (i-1)+ao!1 (i )+a+g (i + 1)) (A3.15) =(b _+b 0 +b +)(a._+a. 0 +a+)/ (i )g (i) +(Az)[(b +-b _)(a._+a. 0 +a.+)/ '(i)g (i)+(a.+-a._)(b_+b 0 +b +)/ (i)g '(i )] + (~ > ( (b ++b _)(a._+a. 0 +a.+)/ 2 (i )g (i)+2(b +-b _)(a.+-a.-)/ '(i )g '(i) 2 +(a++a._)(b_+b 0 +b+)I (i)g 2 (i)+ · · · 'lb.us, setting the first coefficient to 1 and the next two coefficients to 0 requires b+-b-=O=a ... -a_ and we can't match the middle term of order (Ax) 2 . If we use two combinations, we can use symmetry to remove all odd powers of (Ax) just as in the expansion of the integral. We can then match all the terms of order (Ax ) 2 and be left with an error of order (Ax }4 . Because the odd powers have a factor (a.+-a_) we simply add two terms with a.+ and a._ switched. (Each with weighting *and b from the previous formula replaced by a..) Thus, we eliminate all odd powers of (Ax) and are left with simple equations for the coefficients a. = (ti_ +a 0 +a+ ) 2 1 (A3 .16) Choosing (a._+a 0 +a.+)=1 simplifies the other equations to become (a.++a_)= ~ (a.+-a-)=± (A3.17) J:s Choosing the minus sign above gives 2 -o-3 ft _ (A3.18) Now if we substitute these coefticients into 1•=a.+/ (i-l)+a.of (i)+a.~1 (i+l) (A3.19) and do a Taylor's expansion, we see that we have derived an approximation for - 155 - I (i+~) with error of order (fix ) 8 . fi=J (i}+(fix) /'(i)+ (fix)2 /2(i)+ ... ~ (A3.20) 2X3 1t1/(i+~) There is some ambiguity, however, to applying this rule to integrating functions of 1/l'from the values of V; at the mesh points. For instance, in calculating the kinetic en~rgy, we calculated V1/I at the points±(~) and then squared. We could also imagine using V1/I at the mesh points, squaring and then interpolating to the points ± (~). Of course, this would be foolish in this case since V1/J isn't known at the mesh points. Calculating it at the points ±{fix) from the three allowed points, a; .Z\ ±(fix) would be less accurate than calculating it inside the interval, for instance at Z\ ± <.e;> . This method would also require two succes- sive numerical approximations and would be more complex. Thus, for both T and ; we want to approximate the gradients directly at the relevant points (A3.21) and ;•="im.[(a+;i(i-l)+a01/l(i )+a:t:'ft(i+ 1)) •(b+t'(i-l)+b 01/l(i )+b :1:1/l(i+ 1))](A3.22) For calculating p+ there is no clear di!lerence between the accuracy of the methods p+(i )=a-.p(i-l)+ac,£>(i)+a:t:P(i + 1) (A3.23) p+(i)= Ia+l'(i-l)+a0.Y,(i)+a:1:1'(i +1) ! 2 (A3.24) However, since we are already calculating the densities at the mesh points (it is our primary result, in fact), and the resulting formulas for both the energy functional and the TDHF Hamiltonian will be much simpler (involving no cross terms - 156 - like 1/l•(i-1)1/l(i)), we interpolate the densities at the mesh points . (A3.25) Our 1-D approximation (A3 .26) becomes in 3-D (A3.27) where Rn':& =*l;(p;11;(i ,j ,k )-r;·z(i ,j ,k )+piz(i,j ,k )Tqo11;(i,j ,k) (A3 .28) iii: -J;(i,j ,k )fq\(i,j ,k)-J~(i,j ,k )fq-;z(i,j ,k)) (We have included the x in Pz also since this determines the direction of the points ±(6.x) .) We must also discretize they and z components to form R 11 and R1e . Then with (A3,29) The net result is (A3 .30) We can simplify the expression for R considerably, since we used a linear combination of p(i,j,k) to make our p•(i,j,k), by rearranging the terms in the sum. In 1-D l;[p+(i);+(i )+p-(i )T-(i)]= ~([ a+p(i-l)+aa0(i)+a_p(i + 1)]T+(i) ' (A3 .31) ' +[ a_p(i-l)+a.o.0(i)+a+P(i +1nr-(i)) : U(i)(a+;+(i + l)+aoT+(i)+ri_;+(i-1)+a_r-(i+ l)+a.01-(i)+a+ ,--(i-1)) ' In 3-D we rearrange the entire sum or pT terms where - 157 - r(i,j ,k )=-r:r (i,j ,k )+-r'll(i ,j ,k )+-r• (i,j ,k) (A3.33) and 'Ts (i,; ,k )=*(a_[ -r;(i+1,j ,k )+ii(i-1,j ,k )] (A3.34) +a 0 [ r;(i,; ,k )+1i(i,j ,k )]+a+[ -ri(i+1,j ,k )+-r;(i-1,j ,k )]) and -r11 and Tz are similarly defined with the index j or k, respectively, varying ±1. This me.ans we only need two extra arrays -r11 rather than p: and -r; in all three directions. As we shall see later, this -r1 is precisely what we need for the TDHF Hamiltonian. (Because 6'1/!{q(~.j ,k) p(i,j ,k )=•rfl-A. (i,j ,k ).) 1 IV. Coulomb and Yukawa Terms The Coulomb and Yukawa energies are simply summed using the densities and the Coulomb and Yukawa potentials. Ern·= ~Pq (i .; ,k) Uyq•(i ,j ,k )(~ )S (A3.35) ijk Vi, Hy= "t-<ElH> +Eim) + Vr.r E:ijm Hc=l/ZD>p(i.j ,k )Uq,(i,j ,k)(~) 8 ii" These potentials, U, are calculated by the solution of discrete Poisson or Helmholtz problems with the conjugate gradient method as described in (FL 78). We also choose the boundary conditions zero for the short range Yukawa force resulting in a repulsion from the edges of the box which destroys energy conservation when the density approaches too near the edge. For the Coulomb boundary conditions, we use a multipole expansion including multipole moments of the mass and charge density through order 3. These are important to calculate anyway since they will be used to divide the density into two spherical fragments and calculate their mass and position at each time step. (See Appendix A7) - 158 Appendix A4: Finite Ditference Form for the TD:HF Equations in Space I. Introduction and Zero-Range Terms The TDHF equations are partial differential equations in time and space . Since we need ordinary differential equations for numerical computations, we first need a spatial discretization of the equations, specifically for hq(i,j,k)= ·t. 61/1>.q i ,3 ,k )fl=hqz+hqx+hqNL+hqy+hqc (A4.1) Then we need to specify the time evolution algorithm. We obtain our discretization by the variation of the discretized form of the energy functional with respect to the single-particle wavefunctions at the mesh points. Thus, our discretized equations will properly embody the conservation laws satisfied by the continuous equations, particularly energy conservation and the orthogonality of the single-particle wavefunctions. We first do the variations in 1-D and generalize to 3-D for which all directions are equivalent. First, we calculate the zero-range piece of the TDHF Hamiltonian It _ • "VZ - 6 rr Z,J, • ( • • k ) ¥l Z u.,.>.cz 'L,J, (A4.2) by varying the densities in Eq. (3.6) d-i/I• t. E k) ( >.q .,, .1 • 1/1. t. E 6 >.q 'I, ,J, at• t. E >.q 1. .J. k >( k) ( P:·(i.',j'.k'))=2.p"(i,j,k)l'>.q(i,j,k)oqcz· i';'I;. (A4.3) Pp(i',j',k')Pn(i',j',k'))=pcz'""J(i,j,k)'l/l>.q(i,j,k) i'j'J:' p(i'.i ',k ')Pp(i •.;· .1, ')p"(i' .; ',k '))=(p 2(i .; ,k )-pi(i.; .k ))'1/1,.q (i .1 .k) i'j'lc' Therefore, (A4.4) h,,z(i ,j ,k )=~0 [(1-'Xo)Pq (i .j ,k )+ (2+xo)Pcz·,oq{i ,j ,k )]+ ~(i,j ,k )-p:(i ,j ,k )]]'1/1>.q (i ,j ,k) -159D. JCi.netic Energy Terms Varying our expression for the kinetic energy (A3. 7-9) gives a 5 point formula for the operator v2 as in (FL 78). ~( ") ~)~[ ,.-(i~+,.+(i')]=*l;[(b_6u·-1+bo6ii'+b +6u·+i)(b-1k (i'-1)+bo1k (i')+b +1k (i '+ 1) 6l'A: 'l. i' ,. +(b +Ou·- 1+bo6tt·+b _6«·+ 1)(b +1• (i'-l)+botk (i')+b-'fl• (i •+1))] = (~) 2 ( 112 ft;(i-2)-~k(i-1)+ ~A:(i)-~t;(i+ 1)+ 112 ¥'t;(i+2)) (A4.5) Thus, in 3-D where ~(i,; ,le)= 112 ,,.,(i-2,; ,le)-~>.q(i-1,; ,k) (A4.7) + ~>.q(i,; ,k)-~>.q(i+l,j ,k)+ 112 ¥1>.q(i+2,j ,le) and Kva, and Kq. are similarly defined withj or k, respectively, varying ±1. ID. Nonlocal Terms The nonlocal terms HNL• Eq. (A3.27-34), are the new contribution. Since this discretization also i.p.vol~ed only the points i,i±1 ( like the kinetic energy), it also gives a 5 point formula in each direction. We begin with ~A:~(t.") D(i')T(i')='T(i)'ft.t (i)+ l;p(i') 61/J ~(") or(i') ,. ,. A: 'I. (A4.8) From Eq. (A3.21), we get 6';~(i) 'l"(i')= 2{!)2 [1t(i-2)( ia )(~·+2+46ii'+l+c5 (A4.9) •.) +1.t (i-1)( tHc5a·+2-l3c5,i·+1-136a·+c5u·-1) +1'1 (i )( ia )(-56a·+2+446ii'+1+ 1026a·+446a·-1-50a·-2) +t.t(i+ 1)( ~(~·+ 1 -136a·-136a:·- 1 +6u·-e>+t.t (i+2)( is )(c5u·+40u·-1 +Ou·-2) - 160 Therefore, the variation of the p-r term gives IS'f~(i) t.J'(i ')T(i ')= (A4.10) Yl1c (i-2)A(i-1)+11c (i-l)B (i )+'f1: (i )C(i )+'J'A: (i + l)B (i + l)+'f1: (i +2)A(i + 1) where A(i)= ?2{!z)2 (p(i-1)+4p(i)+p(i+1)) (A4.11) B(i)= !B{!z)2 (p(i-2)-13p(i-1)-13p(i)+p(i+l)) C(i)=r(i)+ (!z) 2 (-5p(i-2)+4-4p(i-1)+102p(i)+44p(i+ 1)-5p(i +2)) 72 Now we must do the ; 2 term. From Eq. (A3.22), with reasoning like above, we get ~( ') ~(;+1 (i')+;J(i'))=-i[ 'f1:(i-2)D(i-1)+1'1: (i-l)E(i) t51J1c 'I. (A4.12) ,. -Yl1c (i + l)E(i +1)-1',_ (i +2)D(i+ 1)] where D(i)=- lZ(~) [;+{i)+j-(i)] (A4.13) E(i)=- (~) b·i+(i)+pj-(i)+pj+(i-1)+7;-(i-l)] 2 and 7=(-~~) (A4.14) ~=-{~ ~) Th.us, we see that the coet!'icients for the 5 points are complex functions of the densities and currents. For the 3-D result, we need tS (A4.15) 41.;,,(i .j ,k) (Rpp +R,.,J and (A4.16) The first, Eq. (A4.15), is straightforward since one term contributes ~(:.j.k) (p1 T 1 -I:) as previously discussed. and the other term is zero. We - 161 can generalize the previous discussion to give t. 61/1.Ml .,,,, ' k) (Rp,, + Rnn )=rqq (i ,j ,k) (A4 .17) ""Ip: (i,j ,k )=[.Ap;(i-1,j ,k }+i.Dp{i-1,j ,k )]1'>.q (i-2,j ,k) (A4.18) where +[Bp(i,j ,k )+iEp(i,j ,k )]'¢1,.,,(i-l,j ,k )+cqz (i,j ,k )'1/1>.q (i,j ,k) +[Bp:(i+1,j ,k )-i.Eqz(i+l,j ,k )]'l/;>.q (i+ l,j ,k) +[Ap; (i+ 1.; ,k )-iDp: (i+ 1,j ,k )]1'>.q (i +2,j ,k) The subscripts (x,y,z) in the coefficients give the index (i,j,k) respectively to vary and the component of; and j to use . For instance, using the definition of D(i) and A(i) from Eq. (A4.13 and 10), Dp(i,j ,k)=--&<l/z(i,j ,k )+lp(i,j ,k )) (A4.19) Ap(i.j,k)= 72(~)2 (pq(i-1,j ,k)+4pq (i,j ,k )+pq (i+ 1,j ,k)) For the second term, Eq. (A4.16), we adopt the convention q ·~q (A4.20) We can see that by adding the terms Rw and R,9 , we have arrived at the same formula as for R'l'l with one exception. The wavefunctions involved are precisely the same -;,.,, , but the coefficients use q' rather than q. We call t. ~·Aq 'L ,, •k) (R,m+R,.p)=rq',(iJ.k) (A4.21) ntqz=[l\r'.s (i-1.j ,k )+iD11 ·.s (i-1,j ,k )]'¥',.,, (i-2.j ,k) (A4.22) where +[B,·z(i,j ,k )+iEq•z(i,j ,k ))'l/;>.q (i-1,j ,k )+c,.z(i,j ,k )'l/;>.1z (i ,j ,k) +[B,·z(i+l.j ,k)-iE,·.s(i+ 1.; ,k )]-;,.,, (i+l.j ,k) +[J\·.s (i +l,j ,k )-iD,·.s (i+ 1,j ,k )]l'!.>.q (i+2,j ,k) - 162 - Th.us, the final result is We can see that since rq and rq' are linear in p:r,j we can evaluate hqNL more easily by combi.ni.ng densities: (A4.24) and likewise for -r and each component of j . We can use these densities in the formulas for the coefficients A, B, C, D, and E in the three directions to form A,B, C, D. and E, respectively, in the three directions. Applying these new coefficients to the formula for rq will give h,,NL directly while cutting the coefficients evaluated in half. IV. Coulomb and Yukawa Terms The variation for Hy and He is straightforward as in previous work (FL 78) . The discretized rorm. or the potentials u (i ,j ,k) must be equivalent to a sum U(iJ ,le)= I; O(ii: ,jf ,/ck' )p(i' ,j' ,le') where Or is a ('j'I:' discrete representation for exp1;!; ~a) (A4.25) e2 and Oc for Ir-T' I . Since the continuous functions 0 are symmetric in r and r', the discrete representation must be also symmetric if we switch primed and unprimed . Thus, from Eq. (A3.35) we get 61~(~.j ,k) E.Yf'q"=[Uy,(i,j ,le ).Sw+ Uy, (i,j ,k)o'l'l" }lpAq(i,j ,k) (A4.26) We can see that these terms are also completely local. Therefore, (A4.27) Similarly, (A4.28) - 163 - Appendix A5: Implementation of the Separable Approximation I. General Features The main saving from using this approximation is , of course, that the wavefunctions which are evolved are only two dimensional. The discretized form of the separable ansatz is YlNi (i ,j ,k .t )=rpN/ (i .j ,t )B a(k) (A5.1) Evolving only the 'P>.q saves an order of magnitude in both storage and computing time . However, the entire problem is not two dimensional, and we need to keep many three dimensional quantities . Using the separable ansatz, we can rewrite the particle, kinetic energy and current densities: (A5.2) a r:;i:(i ,j ,k )= "Er:=(i,j)B~ (k) a 'T:V(i ,j ,k )= "Er:ay(i ,j )B! (k) a 2 'T:S(i,j ,k )= "f;>qa(i,j)D; (k) Cl ;~(i.i ,k)=Dl=(i.J)B:(k) Cl i:V(i,j ,k )="Eiqa11 B: (k) Cl and il.(i,j ,k )=O The z current is zero since the 19'Ni(i,j)j 2 factors out and Ba(k) are real. The new densities are : P,a(i,;)= "E 19'>.q(i,j)J2 (A5.3) h&Jilh a(>.)=a r;= (i .j)= I; I b+9?>.q (i-1,J)+bo()O><q (i .j)+b ±'P>-q (i + 1,J) 12 (A5.4) Awith a(>.)=a il= (i ,j)= ~ Im[( a~i;oNi (i-1,j)+ao()O>.q (i ,j)+a±rp>.q (i + 1,j)r AwWt ca(>.)=a [b+cpNi (i-1,j)+bofP>.q (i ,j)+b ±!IO>.q (i+ l,j) J) with r:a11 and jq+a11 the same except using j±1 instead of i±1 and (A5.5) - 164- D! (k )=[b+Ba(k-l)+boBa(k )+b±(k + 1)]2 (A5.6) Some of the savings are very obvious. First, ;; is eliminated and, second, -r: becomes trivial since Pqa(i,;) is already calculated and the D! (k) need be calculated only one time along with the basis functions Ba(k ). The savings for the other densities are not as substantial since we need to store them separately for the di.fierent z quantum numbers, a. This saves some storage since we use fewer numbers a than points k, and saves more computing time since each density is summed only over wavefunctions with a particular z quantum number rather than over all wavefunctions . This fact also indicates that the most efficient grouping of wavefunctions is by the z quantum number. The other quantities we need for the evaluation of the energy functional and TDHF Hamiltonian are the Yukawa and Coulomb potentials. These are calculated from the three-dimensional proton and neutron densities without any changes due to the separable ansatz, so we must still calculate these densities Pq (i ,j ,k) . II. Incorporation into the Energy Ji"unctional Since we still calculate the proton and neutron densities, and the Coulomb and Yukawa potentials in three dimensions, the evaluation of the zero-range component, Hz, and the Coulomb and Yukawa components, He and Hy, is unafi'ected by the separable approximation. Hz~~~ 1-:i'ollP:(i J ,k) +p:(i J ,k)J+'l'o( I +1€1' )p,(i J ,k )p. (i J ,k) (A5. 7) 0 + ~(i,j ,lc)+Pn(i,j ,k)}pp(i,j ,Jc)pn(i,j ,k)l(L\%) Hc=*LPp(i,j ,k) Uq,(i,j ,k )(6%) 9 9 (A5.8) ~ "'I. Hy= f{ElPP+Eim)+ VvE]jm (A5.9) where LJ>cz (i,j ,k )+ Uycz' (i ,j ,Jc )(Ax )9 EJW = ~ (A5.10) -165- The kinetic energy term, HK, is simplified. In the x direction, (A5.ll) since I;B!(k)(~)=l, and similarly for y. Th.ere is a slight savings because It; fewer orbitals a are used than points k. In the z direction: I;'T:.r (i ,j ,k )=N,,a Tz(a) (A5.12) 'iii: where Tz(a)=(~)I;}!{D~(k)+D:"(k)] (A5.13) I: is the kinetic energy in the z direction for each basis state and (A5.14) is the number of particles with isospin q and z quantum number a. Now, the result is (A5.15) and the z component need be calculated only one time along with the basis functions.+ The nonlocal terms, HNL; were given by Eqs. (A3.28-34). We look first at the p; terms and substitute the separable approximation. We can do the sums over k to form BM(a,a')=~B: (k )B:· (k )(~) (A5 .16) I: and DM(a,a')=~Da·(k )B! (k )(6.-r) (A5.17) I& where +We don't actually calculate the total kinetic enew in this way, instead, we U8e the Bingle particle energies •A=<1Alhl1A>· ESPT=~e,.. and ETOT=<TJHjir>. Then >. E'l'OT=ESPT-EPOT-E'l'3 where EPOT is the potential energy and ET3 is the three-body term (triple counted in ESP'/'). Therefore EKIN=ETOT-EPOT. The calculation of the "A is str&ghtfo:nrard followin& the calculation of ht as d.iscUDed le.ter in Appendix A6. - 1662 2 Da.(k )=a.-[D;2 (k +l)+Dci (k-l)]+a. 0 [D;2 (k )+D; (k )] 2 (A5.18) 2 +a+[D: (k +l)+D; (k-1)] Iben with 'TfH (i,j)=a_[ 'Tf"u(i+ 1,;)+'T;u(i-1,j)] (A5.19) +ao[ Tia.z(i ,;)+'T;az(i ,;)] +a+[ 'T';az(i+ 1,j)+'T'qaz(i-1,j)] and a similar expression for 'T''latJI, we calculate T1111(i,;)= ~ ([ r 9 a.·z (i ,j)+r9 a·11 (i ,j )]BM(a,a')+p 90 ·(i,j)DM(a,o:')) (A5.20) Once again, the z direction is only calculated once at the initial time . This quantity is stored since it will be useful for calculating h"/I . We can then write ~q (i,j ,k )rI( (i,j ,k )=°E°EP11 a(i,j)Tq'a(i,j)(&r) 2 iJlr. 0. f.J (A5.2i) We rewrite the j terms in a similar way, combining the functions of k and doing the sum over k. We form 19 a(i,j)=°Ei11 a·(i,j)BM(a,a') (A5 .22) (&r )3 °EJq (i .j ,k )jrf (i,j ,k )=(&r )2 °E'Eiqo:(i,j )J9•o;(i ,j) (A5 .23) a' and then o.f.J ~ for each of the four currents (x or y) and (+or-). We store 11 a.(i.j) since it will be useful in calculating h"/I. However, we do not need iva.(i,j), since we can do the sums needed to calculate the energy while constructing the J1 a.(i ,j). First. calculate the currents for each wavefunction, i"Nz · 'lberefore, for each wavefunction we can form the partial sums for each of the four currents (A5.24) We can write E 1111 =(flz ) 2 °E °E°Ei~ (i,j)[j~ (i,j)BM(a(>.).a(>.))+2J;; 1 (i,j)] (A5.25) ~fJ.>.. so that we can do the sum for E1f!l without storing j~(i,j) . For q' ,tq that is Epn. - 167 - the sum is simpler. Since we do all proton wavefunctions first, we have stored the full Jpo.(i,j) . We can then sum (A5.26) The total is then (A5 .27) Ill. In.corporation into the TDHF Hamiltonian The discretized version of the TDHF equation is given in Appendix A4. The separable approximation gives: ih:t 1 'h.q (i,j )=h qu.(i ,j) (A5 .28) where h' 1 u.(i ,j)= ~B 0 (k )liq (i ,j ,k )(.6.:r) (A5.29) t Notice that since the operator h' depends on a., we have a second reason for grouping wavefunctions by their z quantum number. The zero-range terms and the Coulomb and Yukawa potentials are very straightforward. First, these are all diagonal. that is h (i,; ,k )=J (i.i ,k )'lfAtz (i,j ,k) (A5 .30) and, second, the function. J, is dependent on proton and neutron densities and Coulomb and Yukawa potentials, all of which are calculated in the full three dimensions, unaffected by the separable approximation We, therefore, combine these elements hqz(i,j ,k )+h,,y(i,j ,k )+hqc(i,j ,k )=/, (i,j ,k )ifA.q (i ,j ,k) where (A5 .31) - 168 - /q(i,j ,k)=~o(l-xo)p,(i,j ,k)+pq·.,.11 (i,j ,k) (A5.32) t lto(l +*Xo)+ f<2p,(i,j ,k )+p, ..,.11 (i,j ,k ))] +Ji Uy,(i,j ,k )+ Vr1Uy, ..,.11 (i,j ,k)+ Uq,(i,j ,k )oqp The efi'ect of the separable approximation is to replace these h by h ' 10 (i,j)=~B! (k )/q (i,j ,k )(.6:.t )~Ni (i,j) (A5.33) k resulting in a two-dimensional potential /', 0 (i,j) for each set of wavefunctions with z quantum number a. The kinetic energy part hvJ< gives a discretized approximation to the Laplacian, V2• Each direction involves ofi'-diagonal elements in that direction only, so hq.Kz(i.j,k) and ~{i,j,Jc) are functions of i and j multiplied by Bo.(k) . This means the projections are trivial, since (Ax)~B!{k)=l. We can write ~ h' (' ')- ~ 1 ,.Kz "'·' -2m <A:t>2 4 ( ' 1 ') wNi ('.,,_2 ·'') -~;>." .,,_ ·' [ 1 (A5.34) +~Nz(i,j)-~Aq(i+l,j)+ 112 ,Ni(i+2,j)] and similarly for y. In the z direction, h' will be diagonal in i,j: h'1 ,ur.{i.j)= ~ (~)e~.>.q(i.j)~B 11 (k) 1 4 5 4 {A5.35) 1 [ l2'8a(k-2)-3Ba(k-1)+ 2Ha(k )-3-Ba{k +1)+'128a(k +2)](Ax) and only h' 1 o.XJt is dependent on the quantum number a.. Note: Although this would be quite easy to calculate, we don't actually use it. For the time evolution, we use an operator (h - <h>) which would eliminate h 'qKz . For the single particle energies e ,., we add a z kinetic energy to the <h> above using the discretization for Jd.z I cUJ;z) 1 as described earlier rather than the less accu2 2 rate discretization of J dzBa(z) d. Bci: (z) Therefore, 2 . (A5.36) - 169 - independent of ex. The new, nonlocal part, h.vNL, contributes non-constant off-diagonal pieces. For the separable approximation, in the x and y directions all the wavefunctions are "/l>..q (i' .j' ,k ), so we may form A'qu (i .j)= ~B: (k )~ (i ,j ,k )(6.z) (A5 .37) le and similarly form B', C', D' and E' in these two directions. In the z direction, all the wavefunctions involve 1/1>.,,z(i,j ,k') so the terms are diagonal in the separable approximation. Of course since j 111 =O, .D. =E111 =O. Looking first at the off-diagonal terms, for a particular a we calculate Pq 4 (i,j)=(6.x)~B: (k )p11 (i,j ,k )=~PqaBM(a,a.') (A5 .38) a' le and then calculate the needed coefficients : (A5.39) and B'qaz (i,j)= lB(!x ) 2 [Pqa(i-2,j)-13Pqa(i-1,j) (A5 .40) -13Pqa(i ,j)+ P1 a(i + l,j)] with similar equations for the y direction. To calculate D' and E', we need the analogous sum 11 a(i,;)=Dqa•(i,j)BM(cx,a') (A5 .41) a' which was mentioned in the last chapter and stored. Thus with (A5.42) we get (A5.43) and - 170 - E'qaz(i,;)=- 2 (~) [·y'J +az(i,j)+p'J,;az(i,j) (A5 .44) 9 +(J'J,taz(i-1,j )+/'jqoz (i-1,j)] The diagonal part is given by G94 (i ,j)rp>.q (i ,j)=[ C'qaz (i ,j)+C''lOll (i ,j )+h 'qo.NLz (i,j)]rp>.q (i ,j) (A5 .45) where (A5.46) and PM(a.,a')=(l::.x +' (A5.47) )! ( )~Ba.(k l ) (B a' 2 (k -2)+4B a' 2 (k -l)+Baa'2 (k )]B (k -2) A: 72 6x 2 4 4 4 1 ( ) [B a' 2 (k-2)-13B 4 a' 2 (k -1)-13B 4 a'2 (k )+B 11 a' 2 (k + l)]B 11 (k -1) 18 l::.x 2 4 + 72 (~ )2 ( -5B a' (k -2) +44B a' (k -1) 2 11 11 2 +102B11 a' 2 (k )+44B4 o.'2 (k +1)-5B 11 a.' 2 (k +2)]Ba.(k) + 18 (~ ) [Ba.a' (k -1)-13B a' (k )-13B a' (k + l)+B a' (k +2)]B (k + 1) 2 4 2 + 72 (!z ) [B (k 2 11 a' 2 2 11 2 11 2 )+4B11 a' 2 (k +1)+Ba.a' 2 (k +2)]B 11 (k +2) 11 l This matrix PM need be calculated only once along with the basis B 11 • The first term can be combined with the .,. terms from C'z and r, 11 (i,;)=(!::.x)~B!(k)t(i,j ,k) c·,, to give (A5.48) A: :( t1+t2 4 . . )Tv''"9a.('l.,,7)+( f1+f2 4 .. )T,a.('L,:J) using the r,, 11 as defined in the last section and stored. Thus, we can write the entire diagonal part - 171 - ~,a(i,j)= T,a(i,j)+ ~Pq 11·(i ,j)PM(a.,a.') (A5.49) ll' + ?2(~) 2 [-5.P,,.(i-2.;)+44.P, 11 (i-1,j) +102P, ..(i,;)+44P,a(i+1,j)-5Pqa(i+2,j) Then, the final result for the nonlocal term is h 'qczNL (i,j )=G, 111 (i,j)r,D>.q (i,j) (A5.50) +[A' 1 cu(i-l,j)+iD'1 cu(i-1,j)]rp>.1 (i-2,j)+[A'1 ,.Jl(i,j-1)+iD'1 "11(i ,j-1)]9'>.q (i,j-2) 1 +[B'1 -(i ,j)+iE'1 u (i,j) ]r,D>.q (i-1,j)+[B'1 "11(i ,j)+iE'11111 (i ,j)]r,D).q (i,j -1) +(B'1 CIS(i+ 1,j )--iE'1 az (i+ 1,j)]rp>.q {i+ 1,j)+(B'1 "11(i ,j + 1)--iE'1 "11(i,j+1) ]r,o>.q (i ,j + 1) +[A'1 az (i+ l,i}--iD'1 iu: (i+ l,j) ]'PAI (i+2,j)+[A' 1 CIJI (i ,j + 1)-iD'q"ll (i ,j + 1) ]'P>.q (i,j +2) and the total discrete representation of h'lf is h' 1 ,.(i,j)= /'qa(i ,j )r,D>.q (i ,j)+h 'qK(i ,j)+h 'gaNL (i ,j) (neglecting the v; operator). (A5.51) - 172Appendix AB: Discretization of the TDHF' Equations in Time I. Introduction and Statement of the Problem The evolution or the TDHF wavefunctions, 1fllz, is given by (A6.1) where 1/IN/=1/lllz(z,y,z,t) in three-dimensional calculations (or 1/IN/(x,y,t) in the separable approximation) and his a hermitian operator, h=h(z,z' ,y,y' ,z,z',t) (or h(z ,z' ,y ,y' ,t) in the separable approximation). For the rest of this appendix, all dependence except time will be suppressed. This equation is formally solved by using time evolution operator, u,, as Yl'Nl = u, (t ,t 0 )1/IN/ (t 0) (A6.2) where u, is the unitary operator . ' U9 (t,t 0 )=Texp[( ~f dt'hy(t')] (A6.3) to with T the time ordering operator. Discretization de.fines a time mesh f m=(m)M where (m) is an integer and calculates 'I/IN/. at each time step from YJN/,(m+l)= u,(m+l,m)'l/l>.q(m) This discrete evolution operator can be written as -i u,(m.+1.m)=exp[( T16.t hy(m+)0] (A6.4) (A6.5) where h,z(m+~ must represent the average effect of h(t') over the interval t to t+llt. Thus, determination of the correct operator .hy(m+*> accomplishes the discretization in time . An obvious' possibility, hy(m.+*)=*(hv (m)+hy (m+ 1)) is very close to the operator we use. Part II of this appendix compares the Taylor expansions in 6t for exact and discrete solutions. This procedure determines an arbitrarily accurate approximation for hy (m+)0 as a function or hy, 'I/IN/., 1/l;q. and functional derivatives of hy with respect to 1 and y/; but, of course, the higher-order terms ·are very - 173 - cumbersome to calculate. Therefore, this expression for hv_ (m+*) is not used in the actual evolution code. Instead, it is used to estimate the error arising from other choices for hy(m+*) . These other choices use the values of h, 'ljl, and 'I/I• at various discrete time steps, rather than derivatives or functional derivatives. We use a two point formula involving the time steps (m) and (m+l). With only two points, this approximation can match only the first two terms in the exact expansion for hq(m+*). This results in an error of order At 2 for h.q(m+}0; or order At 3 for Uq(m+l,m) (since the first order term is 1). There is an ambiguity in this two point formula since changing terms of order At 2 or higher won't change the order of the error. This ambiguity is resolved by choosing the higher-order terms to conserve energy. Part III of this appendix calculates the energy difference E(m+1)-E(m) and expresses it in terms of a single-particle operator 01 as E(m.+ 1)-E(m)=I;[ 'l/J~ (m+ 1)01 1/IAq (m+ 1)-'l/l~(m)01 '1/1>.q (m)] (A6.6) Aq This means that M'=~'ji~(m)[ U1,+01 Uq-Oq]l!>.q(m) (A6. 7) "' (where U is the evolution operator as before) and energy is conserved if [U1 ,Dq}=O From Eq. (AS.4), choosing h.q(m.+~=Oq conserves energy. Part Ill also checks that the operator h.q(m+*' determined this way matches the expansion through order At as expected for the best two point formula. These calculations are based on those of (FL 78) . A problem with choosing ftv (m+*> to depend on quantities at time steps (m) and (m.+1) is that the quantities at step (m+l) are not yet known (they depend on hq (m+}0) . Part N of this appendix discusses the algorithm used to achieve this self-consistency and other details of the actual evolution code. This kind of - 174- self-consistency is necessary because the simplest guess hq (m+}0=hq (m) doesn't work. This operator would evolve 'l/J>.q in a potential always lagging behind its proper position. This gradually slows motion and results in a continuous decrease in energy. D. Comparison of Taylor Expansions for the Exact and Discrete Solutions The exact solution at t +flt can be expanded as (A6.8) where the exact solution obeys the TDHF equation and h=h(t) because h is a function or 1/'>.q (t ). The various derivatives of 'l/J>.q can be written exactly in terms of derivatives of h as a -i ( ft-11/l>.q =(-,;1hv i/l>.q <ft121/l>.q=[(-::)2Jtf+( -::)( (A6. 9) :t )hq]i/l>.q etc. Th.e time derivatives of h can be written exactly in terms of functional derivatives as (3P11,, =~.[[ 6:,.., r..J(/piV,..·+ [6:;., hv right) ( fi-1JV;..] =( (A6.10) ~~·[! ~,... r..]11,,-t,..·-[6:;.. hv ]11, -.PA.] (3P'11,, =( -; }~· [ 6:,... [(/tl1i,, lh,,·JV,.,,.- 6:... [( fi-111,, ]liq ·iVA.·] etc. The discrete TDHF solution is given by ?~ (m+ l)=exp[ ( -:: )flt h.q (m+}0]i/l>.q (m) (A6.11) In the expansion, h.q (m+J.0 is time dependent, involving h at future times (A6.12) - 175 - The net result is (A6 .13) 1Ni (m+1)=1Ji.>.q (m)+~t (";::)/ o~>.q(m)+.6t 2 [}2( ; )2/ 9+( ; ) / 1 ]~.>.q (m)+ · · · Matching this to the exact solution determines the operators f"; and gives hq(m+}0 accurate to any order. The first-order terms give /c=h.q(t) (A6.14) The s.econd-order terms give (A6.15) The third-order terms give /2=t< IP2,,_~ ~[liq( IP11v-<( IPhq)hq] (A6.16) t< 1Ps11.q -t< ~) [hv (~211.q-( ( IP21iq )liq] (A6 .17) The fourth-order terms give f s= etc. So now given hq and 'f/!.>.q• we can compute the operator hy(m+}0 and compute to any order 1Ji>.q(m+l). However, we want to derive an approximate formula based on hy(m) and l!v(m+l) . The easiest way is to look for conservation of energy. m. Derivation of the Discrete Evolution Operator from Energy Conservation Defining (A6.18) 6pq=pq(m.+1)-pq(m), c5Tq='Tq(m+1)-T9 (m). c5fq=fq(m+l)-f9 (m), and taking E(m.)=<H> from Eq. (2.7) or (A2.29), the energy difference between two time steps can be written - 176 (A6.19) E(m+l)-E(m)= J ctr[,! [o'Tp+OTn]+ T~L l6pp[Pp(m+l)+pp(m)]+6pn[pn(m+l)+pn(m)]~ To!'! + ~opp[pn (m+l)+pn(m)]+opn[pp (m+ l)+pp (m)]~ 2 1 1 2 + Ts~OPp [,W'P (m+ 1)Pn (m+ 1) + W'P (m+ 1)Pn (m)+ W'P (m)pn (m+ 1)+ S"PP (m)pn (m) + ~:(m.+1)+ ?n (m+ l)Pn (m)+ ?:(m)] 2 1 1 2 +6pn[ S"P (m+l)pn(m+l)+ ~.P {m+l)pn(m.)+ ?p(m)pn(m.+1)+ W'P (m)pn (m) + ~(m+ 1)+ ?P (m)pp (m+ 1)+ ?i(m) ]~ CL + "t-16pp[ 'Tl' (m.+l)+rp (m)]+o'Tp [pp(m+ 1)+pp (m)] +opn ['Tn (m+ l)+'Tn (m)]+orn lPn (m+ l)+Pn (m)] -6Jp [j~ (m+ l)+jp (m) ]J Cy + 2'°10Pp [ 'in(m+1)+'Tn(m)]+c5pn ['Tp (m.+l)+'Tp (m)] +o'T,, (pn (m+l)+Pn (m)]+61n[pp (m+l)+pp (m)] -26]p [h (m+l)+Jn(m)]-26Jn [JP (m+l)+Jp (m)]~ ~ + "°i'16Pp [ U1P (m+ 1)+ U.Jj1(m)]+6pn [ Uin (m+l)+ Uin (m)]~ Yu + -t-topP [Uin(m+ 1)+ Uin(m.)]+opn[ ulj1 (m.+1)+ u.Jj1 (m.)]~ +opp*[ Uc(m+l)+ Uc(m)J) We would like to write this difference in energy as ~(m.+ 1)0.'lfAiz (m+ l)-f'At(m.)09 .Y,Alz (m) to choose an evolution operator which conserves {A6.20) energy. Since 1'Aiz(m.+l)=Uq'f'Ni(m.), the energy ditference can then be expressed as (A6.21} ·Thus, for any u, which commutes with 09 energy is conserved exactly. Since we have Uq=exp(( -;)flt 12.q(m+*)), choosing hy(m+}0=01 conserves energy. (A6.22) Therefore, - 177 (A6 .23) Oq=l Similarly, for 6pqF(r) (A6.24) .. For 67' and 6jq the calculation requires integrating by parts as in Section 2. 1 Therefore, J d:r6'1'q(r)F(r) yields Oq=-V·FV (A6.25) where Voperates on everything to the right, and f d:r6]q(r)F(r0) yields Oq=FV~VF {A6.26) Applying these Eqs. (A6.22-25), we rewrite the energy difference as J AE= ttr2;[ li~(m+l)Oq'\VJ.q (m+ 1)-l'~(m)Oq'\lt>.q (m)] "q (A6.27) where with the exception of the three-body T 5 terms; Oq =*[hv (m+l)+h.q (m)]=h:z (A6 .28) with (A6.29) ii:, =h(p,, =*[p1 (m+l)+p9 (m)].rq=*['1"q(m+ 1)+'1"q (m)].Jq=~Jq (m+ l)+Jq (m)]) 'lb.e three-body terms are TaOf= T8 [ ~:·(m+l)+~,,.(m)pq•(m+l)+ }P:·(m) (A6.30) 2 1 1 2 + 3Pt (m.+ l)p,,.(m+l)+ s-Pq (m+l)p,,.(m.)+ 31'q(m)pq•(m+1)+ s-Pq (m)pq.(m)] - where rt ~q. Since we will use h'l, that is a calculation of the discretized version of h for the average densities of tbe two time steps, we must calculate the correction C1 , due to three-body forces O" =~ + C11 . Cv =Ts[ of-(2*[p" (m+ 1)+p11 (m)]~q.(m+ l)+p11 ·(m)]+[}Hp,·(m+ l)+p,.(m)]] 2) =~~ [(6p1 .)2+26p9 6p9 ·] (A6.31) Thus, choosing fly (m+*>=hq +C9 conserves energy exactly, but how accurately does this operator approximate the TDHF evolution? - 178 - To check this, we write 011 as Oq =hq (m+)0=*!:h,, (m+l)+h.q (m)]+C11 (A6.32) with a difierent correction term, where again C11 only involves three-body terms, Ts. (AS.33) Now we can expand dp dp q =flt~ ... dt (A6.34) h.q(m+l)=h.q(m)+6t~ '1~ d~';t 1 · • · (AS.35) and 2 giving q 2 (m)+M* dhq +6t2[ 1 d. 1tq - Ts [( d.pq )2+2 dpq dpq· ]] ., dt 4 rJJ;2 4 dt dt dt Q :n (AS .36) By comparing this expression with our previous expansion for hq (m+*), Eqs . (A6.14-17), we see that terms through M match, meaning that 1f;(t +dt) matches through terms of '1t 2 . With or without the correction term c, or c,. the error is order '1t 3 , but including the correction term conserves energy for this approximate TDHF evolution. thus providing us with a check on the numerical accuracy of the program. IV. Evolution .Algorithm We must calculate '¥'>.t(m+l)=exp[(-;:)M .h.q(m+)0}Yi>-q(m) for each orbital, with the Hamiltonian for the average densities and a correction term, as dis- cussed previously. There is a self-consistency problem because the operator hq(m+*) depends on the wavefunctions at the next time step (m+l) . A doublestepping procedure gives sufficiently accurate sell-consistency (as measured by energy conservation). First, the wavefunctions 1/J>.11 (m+1) are estimated usiilg - 179- lzv(m+}0=h.q(m). Then, these estimates for "A.q(m+l) are used to construct the lzv(m.+*) which performs the actual calculation of 1'1-q(m+l). Further iterations could be done if necessary. The exponential is expanded in a power series and the evolution operator is modified slightly to give better convergence . This modification consists of adding a phase k6t <hq >. This phase is physically irrelevant (since the particles are independent), but now the evolution operator becomes U11 =exp[ ( -;; ).6.t (hg-<hq >)] (A6.37) The exponent is now a "smaller" operator and the expansions of the exponential are more accurate . The wavefunctions , 1'A.i (m) are used to construct the densities p 11 ( m), ,-1 (m) , and J,(m). These densities are stored and used to calculate the total energies He. Hy, HNL and Hz = ETO and ET'3 at time step (m) (as discussed in Appendices A3 and A5 part II) and to calculate the operators hv_ (m) (as discussed in Appendices A4 and A5 part III) . This Hamiltonian is used to calculate "' -i '1A.i(m.+1)=( ~.6.t h.q(m)1',.q(m) (A6.38) and the single-particle energies 1h - (m+ 1)> eAll =<hy (m)>= A£<1'1'q (m.) l"Aq (A6.39) Th.e exponential of the modified evolution operator can be expanded yielding '1Ni(m+l)=~ ((~At [11.,, (m)-<hq (m)> ]).t =~1'fq(m+l) k! t;'>q(m) (A6 .40) I " where 'f',fg(m+l)=lxl'Aq(m) . Therefore, the quantities ~,.11 (m+l) and <h11 (m)> are used to calculate the first term, (A6 .41) - 180 - The energies e >..q are also used to calculate the total single-particle Inergy and total energy (as discussed in Appendix A5 part II) . Higher-order terms in the expansion are calculated from the previous term by (A6 .42) and a running sum is kept 1'Ni(m+l)= t 1/tf~(m.+1) . (A6.43) t'=O (This requires storage space for two wavefunctions - the most recent term, and the running sum.) For this first step, two terms are typically used in expansion of the exponential. From the estimate of 1/t>.q(m+l), we construct the average densities by returning to the same routine used for the densities in the first step, but not initializing the sums over the wavefunctions to zero. Instead, the final results are multiplied by one half. These 1/t>.q (m+ 1) are also used to calculate the correction term, C11 , from *6p 11 =pq-p11 (m) before storing the average densities where p11 (m) was previously stored. Now the same routines are used to calculate hv(m+*), except for the addition of the correction term. However, calculation of the energies is skipped. The same evolution routine is used to calculate the new wavefunctions, '1Ni (m+ 1)=exp[( ~)at (h,, (m.+}0-<J; (m+}0>)]1/tN/(m) (A6 .44) but for the actual time step four or five terms are used in the expansion of the exponential. Also, the single-particle energies are unchanged: they are not set to <hy(m+~> and the calculation of the total energies is skipped. - 181 - Appendix A7: Bi.nary Partition of. the Nuclear Density 'We choose a convenient r,eference frame centered at the total center of mass with axes parallel to the principle axes of inertia. These are easily found from the spatial distribution of nuclear matter . Due to the assumed refiection symmetry through the scattering plane, one axis OZ is perpendicular to the scattering plane. The remaining two axes in the scattering plane are denoted by OX and OY, OX being associated with the smallest moment of inertia (see Fig . 32a). For symmetric systems such as 160 + 111 0, the analysis of the nuclear density is rather simple . The center of mass, 0, is a center of inversion symmetry so that the Y,Z plane provides a natural division of the system. However, for asymmetric systems this is not an adequate procedure. Ai!. shown in Fig. 32, the total center of mass may be well within the heaviest fragment even when the two nuclei are physically separted and no longer interact through the short-range nuclear forces . Under such circumstances, it is desirable to divide the system into two pieces by a plane between the two fragments parallel to the Y, Z plane and in some sense as far as possible from each of them; for instance a plane through the dashed line in Fig . 32. 'We present two systematic procedures to automatically define this dividing line. The first one was used in the O+Ca calculations (see Section 7) and in the first few Kr + La calculations (see Section B). The second one has since been introduced into most of our Kr + La calculations. Th.ere are no significant difierences between the two methods , which are only useful tools to analyze TDHF solutions. However, the second one is better suited to very asymmetric events found in many IDHF studies which result in scattering for angular momenta below the fusion window. In both methods, we analyze the system in terms of two spheres with -182- uniform density, each centered on the X axis. These spheres are characterized by four parameters: the positions :r 1 and :r 2 of their centers along the X axis and their masses, A 1 and A 2 . Equivalently, their radii R 1 and R 2 are known if their density is given. The spheres are constrained to have the same total mass, A, and center of mass position ( =O) as those of the actual system: Ai+A2=A (A7.1) Aiz1+Ai!:r2=0 (A7.2) The two-sphere system must also reproduce the mean square extension along the X axis (A7.3) with <X2>= j dXdYdZp(X, Y,Z)x2 Finally, the sum of the mean square extensions in the Y and Z directions must be reproduced (A7.4) We now need only a relationship between Rand A to solve these equations. Jn the first method, a uniform time-independent density is assumed for the two sphere, namely (A7.5) where r 0 is ti.me independent and is determined by the initial condition through F.q. (A7.4), where A1 and A 2 are set to the masses of the colliding ions. Substitu- tion of (A7.5) into (A7.4) and a combining of (A7.3) and (A?.4) then results in a set of four equations with four unknowns (A7.6) (A7.7) - 103 ~ A 7. 8) Eqs . (A? .6) and (A7.7) are solved iteratively for A 1 and A 2 at er:>.c h t.lme s~sf , ~Ith the initial guess being the values of A 1 and A 2 from the previous time step . For convergence, the smaller A is substituted into Eq. (A7.7) and the new value of the larger A is substituted into Eq. (A7 .6) . Eqs. (A7.8) and {A7.9) are then solve d for.xf andxl : (A7.10) and These solutions do not specify the signs of .x 1 and x 2 . However, the two spheres must be on opposite sides of the total center of mass, i.e . .x 1x 2 < 0 . To completely determine the signs, we require that the expectation value ~>has the same sign in both the actual and two-sphere systems. - Having defined these two spheres, we divide the system along a line equidistant from the surface of each of the spheres, as illustrated in Fig. 33, i.e . at z+(R 1 -R2}/2 where z=(z 1+z 2)/2. (The location of this line is really the only quantity of interest in the procedure.) It should be noted that this procedure does not depend upon the orientation of the X axis, which is, of course, changing throughout the collision. In solving Eqs. (A7.6) and (A7.7), difficulties may arise if the coefficient P is 5 either too large or too small, since a solution is possible only for P between A 3 !.. and 2( ~ 5 . Since the definition of the axes requires <Y2~<X2> and <Z 2 > is frozen by the separable approximation, P will be large when a fused system is formed with an almost oblate shape: <X2>~<Y2>. This will then result in a very asymmetric solution with one mass, say A1, almost equal to A. and the other , A 2 , unphysically small, with the sphere corresponding to A2 unphysically distant ( <z~ > large ). - 164- To avoid this difficulty, whenever A2 (or A1) becomes smaller than a specified critical mass Amm +, a binary division is not attempted. Rather, the analysis is done as a single nucleus and the Coulomb boundary conditions are calculated assuming a single fused system. The maxi.mum value of P for a binary division is then (A7. l 1) The other extreme is when P is too small, as for a system which is highly elongated along the X axis . Whenever P becomes smaller than cs 3 Pmm=2(~ 2 during the collision, we have arbitrarily chosenA 1 and A2 to be (A7.12) _A Au-;t±l The above prescription in which r 0 is time independent, has proven to be adequate for scattering events of impact parameters larger than the maximum one for fusion. The accurate determination of an l-window for fusion does not require a refined division of the density, as we only need to know whether or not fusion is achieved. However, for scattering events at angular momentum smaller than l<, the outgoing fragments are highly excited and deformed and <X9> for the actual system is then quite different from its value for the twosphere system defined above . The ansatz o'f time-independent uniform density for the spheres is then not appropriate for those events. This defect is remedied in the second method, which we describe now. We shall use Eq. (A7.1) through (A7.5) but instead of having r 0 time independent, we require that the actual <.K-i> be reproduced by the two spheres. This furnishes a subsidiary condition which determines r 0 as a function of time: +Note that Am1n must be smaller than the lightest fnlgment to be emitted. Othenrise,when the fused 97Blem 9CiBsiOilll with a 'li&ht fragment mess smaller than AllllAo the analysis of the deDHity would be done Ill! i! the system had remained fused. - ~ 85 - Ai(xr+ ~rx1)+A2(x~+~~xz)=<X3> u.7 i J) The system of Eqs . (A7 . l) through (A7.5) (or equivalently (A7 .5-9) togethc v.r:.t:-, Eq. \A7 . l3)) can then b~ slated implicitly. E4s. ,A7.6,8,9} are combinc:d Li give: A 1 2 -Ai1± • -2' 1 v'1+QI2Ax 2 1 (A7.14) and Eqs. (A7.13, 14) give (A7.15) ,..2 _ _ _ P'__ •o- Ii. IL i.e. P'=Pr 2 and Eq. (A7.5) 0 (A7.16) (AP +A~) is unchanged. Note that x, R 1 and R 2 are the only quantities of int.erest. Having a first guess for the two masses A 1 and A 2 , r 0 is obtained from (A7.16), the two radii from (A7.17), and then x through (A7.1 5) . This procedure is then iterated, the new masses being given by (A7.14) . Although this algorithm may not be rigorously convergent for all values of Q, P, and <X3>, due to the physical nature of the problem these quantities are not strictly independent parameters and there is convergence in all practical situations. The convergence is rather slow but can be accelerated by averaging the x's obtained before and after each iteration. This alternate method otrers three advantages over the previous one. First, there are no limiting cases; if the system fuses , x goes to zero, the masses A 1 and A 2 become close to each other (s::i~ but Eq. (A7.15) is always soluble. Second. xis directly obtained, not j% 1 l and jx 2 j . Finally. less emphasis is put on the mean square extension perpendicular to the X axis and the asymmetry of the system is fully exploited. If the system is symmetric, <x5> goes to zero, but the denominator in Eq. (A7.15) has a finite limit, so that x is zero, as it should be. - 186Appendix Bl: Expectation Values Needed for TDRPA We calculate the various expectation values used in the derivation of our time-dependent generalization of the RPA theory. The normalization condition required evaluating several expectation values. ... • ... .. ...+ ............+ 1']=<rp Irp>+<c ·A>+<c ·A >+<c ·Ac ·A > (Bl.1) +*[ <(c··.4) 2>+<(c ·.4+)2 >J Since A+ creates a particle-bole pair and A annihilates a particle-hole pair, all terms with unequal numbers of A+ and A give zero. The only term left is (Bl .2) Therefore, fJ=l+c··c . (B1.3) Precisely the same reasoning is used in evaluating the ~ terms in the first step of <it Ii :t I it> which is to second order (Bl.4) _jJ_ ..!!..J .. • :t) + ... • \. =-n ;::,~ ,c ·c; c ·c dt This same reasoning is used in the orthogonality condition < O=<t lii.(1-}{C··c+c ·.A+ :t I '1l °'1'>=0. +*<c ...4+)2 ) I4>> (Bl.5) Since (Bl .6) can be written (Bl .7) ... *I .. + the only term with equal numbers of A+ and A is <h ·A c ·A >. Therefore, (Bl .8) - 187 - Now we need to calculate the expectation values used to get Eq. (3.24) in Part II . We first calculate the various expectation values for the full two-body Hamiltonian H . Then we are able to write down the expecatation values for the evolution operator h from the one-body terms in H if we simply remember that h has only hole-particle and particle-hole components. Throughout these derivations i, j, k, and l label general states whereas p or h refer specifically to empty or occupied states of the 'I'DHF solution. For the two-body Hin Eq. (3.7). (Bl.9) Since <"-'+ti..: >=c5 .. """' , '11io1a • (Bl.10) (meaning 6v plus the condition that both states must be normally occupied, therefore labelled by h) and (Bl.11) we get (B1.12) The next expectation (Bl.13) Since (Bl.14) and (Bl.15) <a.,tu,, ~+a,+a., ~ >=c5pi( c51i.,,,.. c5""' -c5jt,.,_ 6u..)-c5,,; (c5il...,_ 61m -011i;,..S 6Ur.) I we get ~H>=tph +l;'V.Ph'M' The next expectation value is " (Bl.16) - 188 - (Bl .17) The one-body terms are (Bl.18) Two-body terms are -1~~ t.J Y~A:t 4 "Jtl <a,l+Up ~+a.;+f1i ak a.p·+a.11.· > . (Bl.19) Using the symmetry properties of l1 to combine various terms such as opi. and -Op; gives a result (BUZO) Therefore, the net result is <ApnHc ·A+>= l;cp'h' [tpp· oM. -t11.·11.0pp· +opp· oM· })11."11." p'h' (Bl.21) h" + l1.J1h'hp' +oM' I; l11ft"p'h." -opp·~ ~·11."M" +opp· oM. ~ ~..11. Hh"h"'] It''. h." h."h"' The next expectation value is (Bl.22) The one-body terms give zero. Therefore, <Apn.Ap·11.· H>= 1 I; 'ViJA:L (-oyi Op; +op·10~(011.·i 6h.k -011:A: ohl) 4 "Jtl {Bl .23) I =V.wM· . Therefore, (Bl.24) (Bl .25) These expressions for expectation values of H can be simplified by using the equation {Bl.26) - 189 ~H>=~ • ~ l!C ·.A+>= L cp'h' ( 6M· hpp· -6w htth + Yp11.·11p·) + cp11. <H> • (Bl.27) (Bl.28) p'h' Now we can write the expectation values for ii. . From Eq. (Bl.12), (Bl.29) From Eq. (Bl.27), (Bl .30) From Eq. (Bl.28) , ~fur.j{+>=Lcyh' [6hh·if,.,p· -6w~·h]+cph <h>=O · (Bl .31) p'h' From Eq. (Bl.24), <A,inc .Ah>=O • .. - (Bl.32) w11h=Lcp·11.· (6hh·hpp· -6pp·hn·h + 'Ypn·11p·) (Bl.33) The evolution equation is thus p'h' + L c;,,,. vpp'hh' -~ . Jlh' - 190 - Appendix B2: Invariance of. the TDRPA Equations We prove that an arbitrary unitary transformation among particle states or among hole states does not change the eigenvalues, "'· of the time-dependent RPA equations. (Repeated indices are summed.) (B2.1) or (B2.2) where Cp11. =~e-i~ + ~e ttr.ii A=PAP (B2.3) B=PBP. ~;p'h' =6M' hp,,•-6pp· fln·11. + Vph'hp' and We transform either particle states (B2.4) or hole states ,.,,, =UM' rp,,,. (B2.5) with a unitary operator u+U=1=uu•. With t=detl~ 11J and '1r=det(XJi.L we get single-particle wavefunctions (B2.6) Therefore, (B2.7) Since the proofs for the two transformations are so similar, only the second, (B2.5), will be shown. Under this tranformation (B2.8) and - 191 (B2.9) Therefore, (B2.10) h 'l'P· =hpp· h ·,,,."' = u;.1r... ht.:·1r.·· u;.~ = h 'p11, U,,,,: hp,..· h ·,,,, =hn·p uh'h Although the position of the U's is arbitrary, they have been chosen for later convenience. The IS functions can be trivially rewritten, (B2.11) The transformation of the potentials; Yp11,·11p· =<ph' I v Ih:p' -p' h > (B2.12) gives (B2.13) and of j;HJ'M' =<pp' I VIM' -h'h> (B2.14) gives ~ "' rl v·irP'M' = uM.. Vw"'""''" un...11.· (B2.15) Therefore, the matrix transformations give A·llA..,,.,,,. =uNN' A,n·.. .,,,,.. u;.!n. =ulllt.... ~... .,,,,.. Uii-11.· • (B2.16) and B 'p11.~1r.· = UM" 13,,hn -,p'h."' u;.~. . (B2 .17) With h 'Jlhh ·,,,, =uM.h,,n· /in'p u,;.~ =k,,nhti;, . (B2.18) the projector transforms as (B2.19) - 192 We define (B2.20) and rewrite these equations as (B2.21) C'=UC A'=UAlf1 B'=UBU.,.1 P'=UPU- 1 p··= u• p• u'"1 The only change in the proof for the first transformation, (B2.4) is that U=U~·6M' · Therefore, A'=UAU-1 (B2.22) B'=UBU.,. 1 (B2.23) and so the first equation, (B2.1), is invariant. W'=A'C +l'J·c·· (B2.24) iUC= U(AC +BC) . (B2.25) becomes The eigenvalue matrix in the new coordinate system is ! e· J=lu oJ[ l' -'fi·· -.A·· o A. B J[u-' J o u• -l'J· -.A• o u ...1 (B2 .26) Since this is a similarity transformation the eigenvalues, e.>, are not changed. - 193 Appendix B3: Expectation Values for the SU(3) Model In this appendix, we calculate the expectation values of several operators for a general TDHF state 1-t'>=Ila.ti I0>. The number operator for each level i is Ga=~fln~Clm· " Remembering that cx+=A- 1a+, we can transform G from the n representation convenient for the levels, using a 's, to the representation convenient for our states, using a's. Thus (B3.1) and <Ga>=~Au-..\i ~6i' 1611 ff (B3.2) n =N.At1At1 To calculate the various expectation values used in the TDRPA analysis, we need to rewrite Hin the more convenient representation using a's. Thus, . ~ ~ . . H = """t,Asf At;· K;.·,. + i..J Vv Afi' Aij-A~· i\i·~i ~i· ftf (B3.3) i,ii'fl(r where Kt.J = ~ ~£X.n;. We could now calculate all the various expectation values. n'a' However, by finding the correct substitutes for t~ and i'iJA:l we are able to substitute directly into our previous general formulas for these various values (Appen- dix Bl). 'Ihis method requires fewer calculations but we have, of course, checked that the results of the complete calculation agree with the shorter transcription derived in this appendix. Thus, t'1 =<i It Ii >=<0 IO'.\t af I0> now becomes tj;t =<O Ia. t a.:!i l O> (B3.4) tor any particle n where k ,l =l,2,3. = tti ~ e,~ A;.~ ~o I~ a,!-K <Xn·i· ~ IO> wr n =I;e,Aa:Aa independent of the particle n . ' (B3.5) - 194- The potential term vijl:t becomes VijA:t=<OJan·1ani: Vari.t:a:.i IO> (B3.6) (the particles n, n' are distinguishable). This expression is true for any n, n' but for n =n' the two consecutive creation or annihilation operators give zero. Therefore, the result is evaluated for n,tn' and multiplied by (1-0nn·) . We have V1:imq= L ~i·mq' ~~~A;m·~X (B3. 7) L <01 an'l anA:a,!'2.l:·an2r a,tsm•ansq-a~a.,t- 11 I0> n2.n5 Using ~5q' to eliminate one of the creation operators to its right, using ~ 2 k. to eliminate one of the annihilation operators to its left, performing the sum over n2 and n3, and remembering that n' ffn gives I; <O>=oq'm [oKK Di·m· 6iq +olA:' 6krn' 611i·] (B3.8) n2.n5 +6w [ 61d:' Otm• 6ml' +61.t' Oi·m· 6krn] The two terms above with 6i·m· do not contribute since L ~ A;m• 6i•m• =LA;• A;i• =Oij r rm• (B3.9) and (BS.10) for all i,j. Thus, we are left with (BS.11) for anyn, n'. We substitute these values into our previous formulas. We also need 6M', which becomes Onn·, and Opp·, which becomes 6nn· 6uv· for a,u'=2.3. Thus, Wi.i=l:VihJ°h becomes, for anyn, " W.i=l,;VA:u1=l;Vi1Ap:A;A; 1At1 x22;(1-6nn·) n• ij n' = ~2Vv(N-l)A;.:Ati.A;1At1 f,f and hcz becomes (B3.12) - 195 (B3 .1 3) for any particle n . The elements hu 1 are used in the TDHF equations; all hkl are used in the TDRPA analysis . With these substitutions. the energy <H> becomes <H>= l:;(t i1+~W i1) (B3.14) n =N[l:;£iAuA/1 +l: Vv (N-l)(A;1) 2 (At1 ) 2) ( ii Substituting directly into the general TDRPA matrices. we get for Apn;p'h' a new matrix Ana;n'a'=Onn·hua·-Onn• oaa'h 11+ Yu11C7'( 1-0nn') ' (B3.15) and for B,m~'h' a new matrix Bna;n'a'= Vau'll ( 1-0nn•) (B3.16) We use these new matrices in our time-dependent RPA equations (B3 .1 7) - 196 - Appendix B4: Energy and Evolution Equation as Functions of the Four Rotation Angles In this appendix we would iike to rewrite some of lhe equations, which were written in terms of matrix elements of A, in terms of the four rotation angles 19 1 , ~2· 1/1 1 • and 1/12. We choose for our Hamiltonian, equal level spacing ( t 1=-t , e2 =0, t 3 =t) and equal interaction strength (Vij= f<1-oij ). We choose for A (B4. l) -cos~sin~ 1 ei°" 1 -sin~2sin~ 1 e '"'2 l+cos 2 ~(co~ 1 -1) sin~cos~2 (co~ 1 -l)e -i(1'i--"f2) C~1 t A= cos~sin~ 1 e ""'11 sin~sin:t'- 1 e -i'l'e sin~co~2 ( cos~ 1 -l)e i(ti1'e) 1+sin2~(cos~ 1 -1) The energy, from Appendix 3, is <H>=N[l:ti l~ 11 2 +l: Vv (N-l)(Aj 1 ) 2 (~"1 ) 2 ] i (B4. 2) ;,; Substituting the values of the matrix elements of A and x= V(N-l) gives t <H>=N t (sin~2 sin~ 1 -cos~ 1 -xsin2 ~ 1 [ cos~ 1 (cos~cos21f1 +sin~cos21f2 ) +cos~sin~sin~ 1 cos2(1/1 1 -1/12)]) (B4.3) Next we must derive the equations of motion for the rotation angles from the two complex equations of motion for ID5itrix elements of A . Because of the unitarity of A, Eq. ( 4. ) becomes (B4.4) for j =2,3. First we evaluate the left-hand side for j =2, il;...4; 1 .A,; 2 =ei°" 1 [i,j 1 co~-i~sirn'-2sin~ 1 It: 41i 1 cos~sirn'- 1 ( 1+cos~ 2 ( co~ 1 -1))-;?2 sin~2 co~2 sin~ 1 (co~ 1 -1)] Next we evaluate the left-hand side for j =3, (B4.5) - 197- i~A; 14: 3 =e •12[ i,j 1 sin"'2+i~cos1' 2 sin1' 1 (B4.6) " ~icos2~si..n1'2sin1' 1 (cos1' 1 -1)~ 2 sin~sin1' 1 ( 1 +sin21';?Sm)] Next we need to evaluate the right-hand side for j =2. 11 (cos"'2sin1' 1 [sin2"'2(cos1' 1-l)+cos1'i]-xsin1' 1cos1' 2 (B4. 7) 2 [ cos21'1 ( 1+cos~(cos1' 1 -l))e -2i't1+cos2 1'1sin21'2 (cos1' 1-1)e - 2 -cos2"'2sin21'1cos1' 1e 2i1'1 +cos 2"'2sin2"'2sin21'1( cos1' 1-1)e 2i C1'i -1'~ -sin2~sin21' 1 cos1' 1 e2i1e+sin2"'2sin2 1' 1 ( 1 +cos21'2(cos1' 1-1))e -2i(fli--1'z)]) h2 1=ee' i" Next we need to evaluate the right-hand side for j =3, h 31 =te ' 12sin"'2sin1' 1 (1 +sin~(cos1' 1 -l)+cos1' 1 -x (B4.8) [ cos2~cos21' 1 (cos1' 1 -l)e -2i1'1 +cos 21'1 (1 +sin~(cos1' 1 -l))e - il'a -cos21'2sin21'1cos1' 1e2i1' 1 +cos2 "'2sin21'1( 1 +sin21':a( cos1' 1-1))e 2i(1'i -ila) -sin21':asin21'1cos1' 1e ett2 +cos 21':asin21':asin21'1(cos1' 1-1)e -2i('ti-1'a>]) 2 If we divide by the common phase e •ti for the equation with j =2 and by e '"2 for the equation with j =3, the real and imaginary parts of the left-hand side are quite obvious. This results in four equations: . . h21 1'1cos1':a-1':asin"'2sin1' 1=Im( ~ e J (B4.9) , · . · . hs1 1'-1 si..n1'2+1'2cos1'2sm1'1 =Im( ~ e a 1 1 co~si..n1' 1 ( 1 +cos2"'2( cos1' 1 -1))+12 sin~cos1'2 sin1' 1 ( cos1' 1-1) =-Re(~~) e i. 11cos2"'2sin"'2sin1'1( cos1'1-1)+12sin1':asin1' 1(1 +sin21':a( cos1'1 -1)) =-Re(~~) 8 We can solve the first two equations for ,j 1 by multipl~ the first by cos1'2 , the second by sin"'2. and adding . · h21 . hs1 1' 1=cos"'2Im( ~+sm"'2lm( ~ a i e e (B4.10) Similarly, {B4.11) - 198- The second equation, (B4.11), looks like it might give trouble for ~ 1 =0 . However, both imaginary parts of h contain sin~ 1 also. The first equation, (B4.10) gives -d., =exsi..mS1 COS~1 (sin2'¥'1 COS~2+sin2;(i2sin2~2) (B4.12) The second equation. Eq. (B4.11), gives {B4.13) ~=excoS17-asi..mS2{sinfy1 (-cos~1)+sin2;(i2cos2~1+sin2(1/11-'¢'2)(-sin2~1)) Now having solved the first two equations for ~ 1 and ~. we must solve the last two equations for .,P 1 and .,P2 . We can simplify some of the mathematics by notic1.ng a common factor co~sin~ 1 in the third equation and sin1':!sin~ 1 in the fourth equation. After dividing by these factors. the determinant of the coefficients is co8'1) 1. Once again. these expressions look like they might cause trouble if cos~ 1 =O but the numerators also contain a co8'1)1. For .,j 1 we get "11 =-&(1 +x[ cos2'¥' 1(-cos ,j, 1 +cos~sin~ 1 )+cos21f2 (sin ~sin~ 1 ) 2 2 +cos2(1' 1 --1/t2)(-sin~sin~ 1 )]) For (B4.14) "12 we get ~=-e(2+xf cos2y 1 cos~sin~ 1 +cos2'¥'2(-cos2i7. 1+sin2~sin~ 1 ) +cos2('jt 1 -'f!2)(-co~sin~ 1 )]) (B4.15) -199- REFERENCES (AD 68) J·.R. Adam, N.R. Lindblad, and C.D. Hendricks, J. Appl. Phys. 39, 5173 (1968) 75) M. Beiner, H. Flocard, Nguyen Van Giai, and P. Quentin, Nucl. Phys. A238, 29 (1975) (BER 75) G.F. Bertsch and S.F. Tsai, Phys. Reports 18, 125 (1975) (BO 76a) P. Bonche, S. Koonin, and J.W. Negele, Phys. Rev. Cl3, 1226 (1976) (BO 76b) P. Bonche, J. de Physique 37 Collogue C5, 213 (1976) (BO 78) P. Bonche, B. Gramnaticos, and S.E. Koonin, Phys. Rev. Cl7, 1700 (1978) (BEI (BO 79) P. Bonche, K.T.R. Davies, B. Flanders, H. Flocard, B. Gram:na.ticos, s·.E. Koonin, S.J. Krieger, and M.S. Weiss, Phys. Rev. flQ, 641 (1979) (BR 60) D.A. Bromley, J.A. Kuehner, and Lett. i, 365 (1960) (BR 78) R.A. Broglia, C.H. Dasso, G. Pollarolo, and A. Winther, Phys. Rev. Lett. 40, 707 (1978) (DAS 79) C.H. Dasso, T. n,fssing, and H.C'. Pauli, Z. Phys. A289, 395 (1979) (DAV 74) K.T.R. Davies, R.J. McCarthy, J.W. Negele, and P.U. Sauer, Phys. Rev. ClO, 2607 (1974) (DAV 78a) K.T.R. Davies, V. Maruhn-Rezwani, S.E. Koorlin, and J.W. Negele, Phys. Rev. Lett. 41, 632 (1978) (DAV 78b) K.T.R. Davies, H.T. Feldmeier, H. Flocard, and M.S . Weiss, Phys. Rev. Cl8, 2631 (1978) (DAV 79) K.T.R. Davies, K.R. Sandhya Devi, and M.R. Strayer, Phys. Rev. C20, 1372 (1979) (DAV 80) K.T.R. Davies, K.R. Sandhya Devi, and M.R. Strayer, Phys. Rev. Lett. 44, 23 (1980) (DAV 81) K.T . R. Davies and S.E . Koonin, Phys. Rev. C23, 2042 (1981) (DE 78a) K.R. Sandhya Devi and M.R. Strayer, J . Phys. G4, 197 (1978) E. Almquist, Phys. Rev. -200- (DE 78b) K.R. Sandhya Devi and M.R. Strayer, Phys. Lett. 77B, 135 (1978) (DI 30) P.A.M. Dirac, Proc. Cam. Phil. Soc. 26, 376 (1930) (DY 80) P. Dyer, M.P. Webb, R.J. Puigh, R. Vandenbosch, T.D. Thomas, and M.S. Zieman, Phys. Rev. C22, 1509 (1980) (EN 75) Y.M. Engel, D.M. Brink, K. Goeke, S.J. Krieger, and D. Vautherin, Nucl. Phys. A249, 215 (1975) (FL 78) H. Flocard, S.E. Koonin, and M.S. Weiss, Phys. Rev. Cl7, 1682 (1978) (GA 75) B. Gatty, D. Guerreau, M. Lefort, X. Tarrago, J. Galin, B. Cauvin, J. Girard, and H. Nifenecker, Nucl. Phys. A253, 511 (1975) (GA 76) J. Galin, B. Gatty, D. Guerreau, M. Lefort, X. Tarrago, R. Babinet, B. Cauvin, J. Girard, and H. Nifenecker, Nucl. Phys. A278, 347 (1976) (GO 50) H. Goldstein, Classical Mechanics, (Addison-Wesley, Cambridge, 1950), p. 83 (HOL 74) G. Holzworth and T. Yukawa, Nucl. (HOO 77) P. Hoodboy and J.W. Negele, Nucl. Phys. A288, 23 (1977) (HU 76) J.R. Huizenga, J.R. Birkelund, W.U. Schroder, K.L. Wolf, and V.E. Viola, Phys. Rev. Lett. 37, 885 (1976) 76) A.K. Kerman and S.E. Koonin, Ann. Phys. 100, 332 (1976) (KE Phys.~' 125 (1974) tt (KO 76) S. Koonin, Phys. Lett. 61B, 227 (1976) (KO 77) S.E. Koonin, K.T.R. Davies, V. Maruhn-Rezwani, H. Feldmeier, S.J. Krieger, and J.W. Negele, Phys. Rev. Cl5, 1359 (1977) (KO 78) S.E. Koonin, B. Flanders, H. Flocard, and M.s. Weiss, Phys. Lett. 77B, 13 (1978) (KO 79a) S'.E. Koonin, Marmal Aid Preprint 5 (l'.979) (KO 79b) S.E. Kooniu• Marmal Aid Preprint 6 (1979) (KRA 77) (KRI 77) 11 II V. Kratz, H. Ahrens, W. Bogl, W. Bruchle, G. Fr~nz, M. Schade 1 , I. Warnecke, G. Wirth, G. Klein and M. Weiss, Phys. Rev. Lett. 39' 984 (1977) S.J. Krieger, Nucl. Phys. A276, 12 (1977) -201- (KR.I 78) S.J. Krieger and K.T.R. Davies, Phys. Rev. Cl8, 2567 (1978) (LE 73) M. Lefort, c. Ngo, J. Peter, and B. Tamain, Nucl. Phys. A216, 166 (1973) (LE 78) M. Lefort and c. Ngo, Ann. Phys. Series 15 V.3, 5 (1978) (LI 70) S.Y. Li, A. Klein, and R.M. Dreizler, J. Math Phys. 11, 975 (1970) {LIP 65) H.J. Lipkin, N. Meshkov, and A.J. Glick, Nucl. Phys. 62, 188 (1965) ~ (MA 70) J. Mathews and R.L'. Walker, Mathematical Methods of Physics, Second Edition, (Benjamin, New York, 1970), p. 107 (NA 79) J.B. Natowitz, G. Doukellis, B. Kolb, G. Rosner, and Th. Walcher, Nukleonika 24, 443 (1979) (NE 70) J.W. Negele, Phys. Rev. £!, 1260 (1970) (NE 72) J.W. Negele and D. Vautherin, Phys. Rev. C5, 1472 (1972) (NE 78) J.W. Negele, S.E. Koonin, P. MHller, J.R. Nix, and A.J. Sierk, Phys. Rev. Cl7, 1098 (1978) (NO 76) W. Norenberg and H.A. Weidenmuller, Introduction to the Theory of Heayy-Ion Collisions, Lecture Notes in Physics, 51, (Springer-Verlag, Heidelberg, 1976) (PI 66) D. Pines and P. Nozieres, The Theorv of Quantum Liquids Vol. 1, (Benjamin, N.Y. 1966), Chapter 1 (SC 77) W.U. SchrHder, J.R. Birkelund, J.R. Huizenga, K.L. Wolf, and V.E. Viola, Jr., Phys. Rev. Cl6, 623 (1977) (SI 76) K. Siwek-Wilczynski and J. Wilczynski, Nucl. Phys. A264, 115 (1976) (SK 56) T.H.R. Skyrme, Phil. Mag. !,, 1043 (1956) (TA 75) B. Tamain, C. Nga, J. Peter, and F. Hanappe, Nucl. Phys. A252, 187 (1975) (TH 61) D.J. Thouless, Nucl. Phys. 21, 225 (1961) (TI 74) C. Titin-Schnaider and P. Quentin, Phys. Lett. 49B, 397 (1974) (TS 74) C.F. Tsang, Physica Scripta lOA, 90 (1974) (VAN 78) R. Vandenbosch, M.P. Webb, P. Dyer, R.J. Puigh, R. Weisfield, T.D. Thomas, and M.S. Zisma.n, Phys. Rev . C17, 1672 (1978) II II -202 (VAU 72) D. Vautherin and D.M. Brink, Phys. Rev. CS, 626 (1972) (VI 79) S.E. Vigdor, D.G. Kovar, P. Sperr, J. Mahoney, A. MenchacaRocha, c. Olmer, and M.S. Zisman, Phys. 'Rev. C20, 2147 (1979) (WI 73) J. Wilczynski, Phys. Lett. 47B, 484 (1973) -203- FIGURE CAPTIONS Figure 1 . Schematic possibilities for the double differential cross section of heavy-i~n reactions. Shaded areas indicate peaks. Fig. la shows fusion-fission; whereas Fig. lb shows deepinelastic reactions. Figure 2. Schematic classical analysis of heavy-ion scattering; each initial angular momentum leads to a definite energy and scattering angle. These can ~e plotted as a trajectory for comparison with experiment as in Fig. 2a; or plotted separately as in Figs. 2b and 2c to determine various macroscopic parameters, such as the interaction time as in Fig. 2d. Figure 3. Rms. radius of the Kr+ La system as a function of time for four initial angular momenta. These show both scattering and fusion events, along with one ambiguous event. Figure 4. Comparison of our calculated trajectory of energy vs. scattering angle for Kr + La with the axially symnetric calculation and the -experimental contours of the double differential cross section. Figure 5. Energy, scattering angle and interaction times as functions of the initial angular momentum for Kr+ La. initial bombarding energy. Ei denotes the -204- Figure 6. Charge and mass transfers, charge and mass widths, and charge to mass ratio as functions of the initial angular momentum. Charge to mass ratios of the original Kr nucleus and the total system are also shown. Figure 7. Schematic illustration of (a) stable and (b) unstable TDHF solutions. Figure 8. Contours of the potential energy surface in the SU(3) model for X = 0 .5 • X= 2 .5 Figure 9. Energy surface for Figure 10. Energy surface for 7( = 5 Figure 11. Energy surface for Figure 12. Energy surface for Jf.= 20 Figure 13. Energy surface for Figure 14. Oscillations of the TDHF solution near the mininru.m of the SU(3) model with X= 10 X= 100 X= 10 and the periods derived from an RPA analysis. Figure 15. TDHF path and TDRPA frequencies for a stable path, in valley, with Figure 16. X= 100. Position is shown at time steps dt = 0 :01 . TDHF path and TDRPA frequencies for an unstable path, across ridges, with X= 100. dt = 0.01. -205 Figure 17. TDHF path and TDRPA frequencies for an intermediate path with 'X= 100. Figure 18. dt = 0.01. TDHF path and TDRPA frequencies for a stable path, in a valley, · with X = 5. dt = 0. l • Figure 19. Comparison of the HF and exact ground state energies. Figure 20. Deviation of the exact ground states from independent particle states, Tr C,-(") . Figure 21. Comparisons of the positions of the HF and exact ground states. Figure 22. Comparisons of stable and unstable TDHF paths with exact paths for X= 5. Figure 23. Comparison of TDHF and exact oscillations near the miniml.Uil for Figure 24. X= 5. Comparison of .TDHF and exact oscillations near the minimum for ){= 10. Figure 25. Comparison of TDHF and exact oscillations near the minimum for Figure 26. Comparison of TDHF and exact oscillations near the minimum for Figure 27. X= 20. "f...= 100. Deviation of the exact solutions from independent particle states for the paths in Figure 26. -206- Figure 28. Comparison of the TDHF and exact positions, 81, for an unstable path with X = 100. Figure 29. Comparison of the TDHF and exact positions, 62, for an unstable path with Figure 30. X= 100. C'o mparison of an unstable TDHF path with exact paths for X= loo. Figure 31. Deviation of the exact solutions from independent particle states for the paths in Figures 28 - 30. Figure 32. Slicing of the densities (projected onto the scattering plane) _for separated and unseparated fragments. Figure 33. Equivalent two-sphere systems for the densities in Figure 32. -207Energy - Mass Correlations in Heavy-Ion Scattering Energy Mass Figure la E.1. - Energy t Mass Figure lb ' bl) r:: E. 1 :>.. bl) "H r- ' - I I-< Q) bl) r:: Q) .w Q) .w .-1 t1I bl) Cl) < u i:: I ~ I r:: iil -I I-< I 00 •.-1 I-< Q) +.J .w ca u CJ) 00 r:: 0 H I 0 N I Fusion Figure 2a I 00 Scattering Angle ~ ca i Coulomb Scattering Figure 2c ~ ~ ~ 0 Ill •.-1 E. 1 Ill r:: >. 0 ~ "H .w ~ < ~ •u ,. ..,j Ill Cll u t1I :>.. I-< ... bO Q) .w Q) s r:: •.-1 Q) HH i:: iil ~ u I I Fusion I .R Figure 2b Coulomb Scattering I 0 Fusion J. Figure 2d Coulomb Scattering 86 TDHF Resu lts for 9.0 Separated - .. - - Two Fragments Kr + 139 La • • •• ' I I • I I I • I • • ' • •• ' ,,...._ 8 ~ 8.0 I I 0\ Ill 0 N ·~ I ::s \ '"O Ill c::: ~ ("°'I / J 7.0 10 30 20 -22 Time (10 sec) 40 ~ 1--1 ::! bO ~ ii. Comparison of TDHF Results with the Experimental Double Differential Cross Section E. 1 300 ~ •··. ~ I ,-... :> ~ '-" = Axially Synmetric Calculations 0 0 '\. 1 = Separable Approximation Calculations 250 ~ I 1-1 100 Q) c: ~ u .u .,.j I 0 ~ N I Q) c: ·" 200 . ~ ~ u -a Q) .u u Q) 1-1 1-1 0 150 u ,, ~ ~ , -:, " - -- - -. . L=80 ' 100 150 _) ______-/ 150 100 -------- 5 100 tr 2 I 20 40 60 C.M. Scattering Angle ( 0 ) 80 1 ..;t Q) 1-1 ;:I bO .,.j r.o. -211Interaction Time (lo- 22 sec) C.M. Total Kinetic Energy (MeV) N 0 0 N 0 0 N V'I 0 w 0 0 t".l ..... Scattering Angle (degrees) I w 0 0 w '°0 0 I-' 0 0 -., t-3 ., l"I Cll O" I-' I-' Cb "Cll > I-ti en ~ 'rj Cb "Cl ~ - i:: I-' \J1 ::i 0 ..... "Cl "Cl " Cll l"I 0 x 0 ~ i:: 0 '1 ~ ..... ::-:: ~ l"I "..... + HI ~ ......... N 0 0 0 ::i ......... Figure 5 I-' w .,i.-< l.O TDHF Results for 0 ·,-j .u ~ 86 Kr + 139 La (Separable Approximation) 0.42 tlJ tlJ ~I Total 0 .u I ~ 0.41 H tU ...c: u 6 I N ..-1 N 0 h I H <ll § z <ll -2 . ..-1 0 o,-j .u H tU p.. -4 I - ..., I I ·~ r ~--- ~ -t v 100 I 150 200 (units of ii) .R 250 4 tlJ ...c: r 'O ·,-j ::s: <ll H .J 2 2 o I ,........._ I 100 I I 150 I 1 I 200 I I 250 ;:l bO 4 , <r '" I......_ .u ~ ...... "' 0 ·,-j ~ 2 -213- t/) r-.. .._,, " ~ r-.. " -e,,....._ II) 0 .._,, (!) er ...... .._,, " Figure 7a ~ c:: ,,....._ ::I en 0 .._,, "().ea ,,....._ rt .._,, ,,....._ ~-e0 ------ .._,, Figure 7b " lb CT ...... (!) x= 0.5 Contours Hf Ne. -1.0 90 I ..;t r-1 N I - .8 - .4 0 .4 .8 x = 2.5 1.0 90 e~ .8 e~ (0) -1.0 & I ( o) I -1.0 -.8 I & I -.4 & I 0 I I .4 I & .8 I I 1 I .8 \ \\\\\.\.''' \. \ \ \. \. ' ' ~ .4 .6 00 Q) 45 .4 S-. ;:I bO •rl 0 45 o,..j µ.. .2 -1.2 -.2 0 0 45 e. (0) 90 Q) S-. ;:I bO µ.. 0 1-1.0 0\ 0 0 0 45 e, lo) 90 -215- 1.= 5 Contours H/NE. ~1.0 -1.2 -1.4 -1.4 -1.0 - .4 0 .4 .8 1 .8 80 .4 e.1. 0 ( ") - .4 -.8 40 -.6 -.2 0 .o -1.0 20 40 e. <·) Figure 10 60 80 -216- X= 10. Contours R/NE -1.0 -1.5 -2.0 -2.5 -2.5 -2.0 -1.0 0 1.0 80 0 min, • 40 -3 .43 -.5 -1.5 0 0 -2.5 0 40 Figure 11 80 -217- )C = 20. Contours HINt. -2 -1 -4 -4 -5 -5 -2 1 0 80 0 -2 82 (0) -4 min . • -6. 72 40 -4 -2 -2 0 0 5 0 40 Figure 12 80 -218- Contours H;NE. /( = 100 min. -1 -3 -6 -12 -18 -24 -24 = -33.34 -18 -12 -6 0 1 0 80 -6 e.l. (0) 0-33 40 12 -6 0 0 0 40 e (0) I Figure 13 80 -219TDHF Positions and RPA Periods X.= 10 50 e, (0) 49 48 41 8,._ (0) 40 39 38 0 .8 1.6 E.t Figure 14 -220- TillF Path with TDRPA Frequencie• ("'1£) 't. • 100 2 (104) 4t • .01 IO 60 el.. (0) 40 20 0 20 40 e. (0) Figure 15 60 80 -221- TDHF Path with TDRPA Frequencies (~/£)2 (104) x= 100 dt = .01 80 40 -1 . 0 -0.9 0 0 40 e. (·' Figure 16 80 -222- TDHF Path with TDRPA Frequencies (C.U/t) 2 (104) dt = .01 X= loo 80 -0.0 +o.o -0.3 e~ (0) 0 0 40 Figure 17 80 -223- TDHF Path with TDRPA Frequencies (~/~) x= s 2 dt "' .1 80 e.1. +11 (0) 40 -2 -211---~_-1~9---+-~--~l,..._______~~-1~3~ 0 0 40 Figure 18 e ( 0 I ) 80 Comparison of Ground State Energies 0 I • X"'2 .. x-s I "" N N -~ E/NE r . - HF Exac ----11.E Exact ___!!f - "X.=10 Exact I I a- .-4 <II ""00;:l ~ µ.. 1 HF --~ -8 ~ - E xact ~ )(•20 4 10 N 15 :?0 -225- Deviation of Exact Ground States from Independent Particle States 1.0 :>(=20 X=l5 ________,,_---- x=10 0.6 X=S 0 .2 10 N Figure 20 20 -226- Position of Ground States 1' HF Solution • Numbers Label N for Exact Solution 50 lo11i5 1of- 40 15 e~ ;t=20 20 ;ta:lO (0) 30 20 ):=5 20 40 50 Figure 21 60 -227- Comparison of Two TDHF Paths with Exact Paths X=5 • I - 40 ... t "' 0 - . 2 Exact N = 10 t "" .05, .1 TDHF Path • • I ' 1lt Exact N • 15 t "" .05, .1 JC - • - x - e:t (o) - 20 .. • x - • .....__ \ ~-HB w2--21 0 • • I 20 I e, l 0) Figure 22 I 40 - -228Comparison of TDHF and Exact Positions X=5 46 e, (o) 44 42 38 36 el. (0) 34 32 0 .4 ft Figure 23 .8 -229Comparison of TDHF and Exact Positions X"" 10 50 48 47 42 41 40 39 tt Figure 24 -230Comparison of TDHF and Exact Positions 54 • • - • N = 10 •-•·•N • 15 -TDHF e. (o) 52 \.: 50 44 • e.1 lo) 42 .4 ft Figure 25 .8 -231Comparison of TDHF and Exact Positions x- 100 55 54 • - • N = 10 53 • - - - • N = 15 - TDHF 46 e~ co) 45 44 43 0 .04 Et Figure 26 .08 -232Deviation from Independent Particle States of the Exact Paths in Figure 26 Figure 27 -233- Comparison of TDHF and Exact Positions z"" ioo N 0 --CD ~ 0 °' 0 0 0 00 \ •\ . , 00 0 \ ,• ~ ,. .' ) .... °' Figure 28 •• • •• I z• z• g • .... .... ' i-,'! I "rj V1 0 -234Comparison of TDHF and Exact Positions z., 100 ,,... l"T - 1\ • • J •I •I 1-3 z z I I-' °' t:I ::I: "rj II • I-' I-' Vl 0 •I\ ' • • ' .,• I • / • Figure 29 -235- Comparison of TDHF and Exact Paths Z-= loo ·-· 80 N = 10 IC•-• N = 15 ~ TDHF e~ (o) 40 0 0 40 Figure 30 80 -236Deviation from Independent P•rticle St•tes ef the Ex.ct Path• ln Figure 30 0 . . N C'I .00 0 .0 ~ ,,.. " .0 00 I~ 2! z I I "' 0 ........ Figure 31 -237- Slicing of the Density \ y x Figure 32a \ y \ x Figure 32b -238- Equivalent Two-Sphere System \ x \ \ \ \ \ Figure 33a x Figure 33b